Primary motor and sensory cortical areas communicate via spatiotemporally coordinated networks at multiple frequencies

Frequency-Specific Modulation of Spike–Spike Coherence (MS) During Task Performance.

Coherence between the spiking of MIo and SIo neurons (MS) at multiple frequencies was dynamically modulated; single neuronal pairs exhibited increases and decreases in coherence as well as firing-rate modulations relative to the onset of tongue-protrusive force (Fig. 2 B and C, and Fig. S1). Similar patterns were observed for the mean MS coherence across neuronal pairs with significant modulation of coherence (Fig. 2D, shuffle test, P < 0.01), consistent with the task modulation of theta coherence previously found within MIo and within SIo (16). Out of all of the possible combinations of paired signals (n = 44,152) from 10 datasets, the proportion of neuronal pairs that showed significant task modulation of MS coherence (i.e., “functional network”) was highest in theta and was significantly lower in the higher frequency bands (Fig. S1E, McNemar test, P < 0.01). The larger theta network (i.e., highest proportion of paired neurons with significant coherence) exhibited the strongest coherence (i.e., highest mean peak coherence), whereas the sparser networks in the alpha/beta/gamma bands exhibited weaker coherence. This was observed for each monkey (Fig. 2E) and for data pooled across monkeys [Kruskal–Wallis peak by frequency, ${\mathit{x}}_{\left(3\text{,}20332\right)}^2=10\text{,}683$, P = 0, post hoc, P < 0.001]. Differences in spectral power across frequencies cannot explain this result as the cross-spectrum is normalized by the autospectra. We also found frequency-specific differences in the time of peak MS coherence; pairs of MIo and SIo neurons exhibited a unimodal distribution of time of peak coherence in theta and a bimodal distribution in all other frequency bands [Fig. 2F, Kruskal–Wallis peak times by frequency, ${\mathit{x}}_{\left(3\text{,}20332\right)}^2=122$, P = 3e-26]. The mean time to peak theta coherence (0.03 s, SD = 0.07) occurred before and at force onset (note that coherence values were aligned to the right edge of the 0.5-s window). This suggests that theta coherence may play a role in generating the tongue-protrusive force. In contrast, the distribution of the times of peak gamma coherence was bimodal based on a two-component Gaussian mixture model having the lowest Akaike information criterion compared with models with one, three, or four components. The bimodal distribution of the peak times of gamma coherence was further confirmed by extending the time windows analyzed relative to force onset (Fig. 2F, Inset, and Fig. S2), indicating a suppression of gamma coherence at force onset. Peaks in gamma coherence [mean (M) = –0.35 s, 0.33 s (SD = 0.02) relative to force onset] occurred when monkeys were getting ready to apply the tongue-protrusive force or when monkeys were decreasing force production after being rewarded. On average, monkeys reached the force target at 0.15 s (SEM = 0.02) after force onset. Force offset typically occurred 0.5 s after force onset. The results of the coherence analyses in the theta band using a 0.5-s window were similar to the coherence results using a 1-s window (Fig. S3), and in the gamma band, the results using a 0.5-s window were similar to the results using a 0.3-s window.

In sum, task-modulated MS coherence differed between low- and high-frequency bands in terms of network size (i.e., the number of paired neurons with significant coherence in a specific band), peak time, and peak magnitude. Such dynamic modulation cannot be attributed to common task-related modulation nor to common inputs received by the neuronal pair because the shuffling procedure effectively estimates these coincident effects. A cross-validation approach (SI Results and Figs. S4 and S5) and performing the same analyses using only unique neuronal pairs further confirmed the dynamic modulation of coherence.

Changes in Spike–Spike Coherence (MS) with Learning.

Fig. 3 A and B illustrates the coherent activity of the population of paired MIo–SIo neurons for D1 and D5 of monkey Y. Proportions of neuronal pairs with significant modulation of MS coherence changed during training from day to day. We tested statistical significance between D1 and D5 for each band and found P < 0.05 (Fig. 3C, binomial test). We also found day-to-day changes in peak and time of peak coherence in the theta band [Fig. 3D and Fig. S6, Kruskal–Wallis, peak by days: ${\mathit{x}}_{\left(4\text{,}11017\right)}^2=53$, P = 7e-11; peak time by days: ${\mathit{x}}_{\left(4\text{,}11017\right)}^2=32$, P = 2e-6] but not in the other frequency bands (Kruskal–Wallis, P > 0.10, peak/time of peak coherence by days separately for alpha/beta/gamma). We did not find any systematic changes in peak and time of peak coherence in relation to changes in tongue protrusion force, success rates, and reaction and movement times (Fig. S7, Pearson’s correlation, P > 0.10). Furthermore, changes in firing rates cannot account for changes in MS coherence; no significant correlation was found between mean firing rates and mean MS coherence (Fig. S8, Pearson’s correlation, P > 0.10). These indicate that coherence does not directly relate to the encoding of specific behavioral parameters and suggests a role of coherence in the spatiotemporal coordination of different functional networks that emerge or are reshaped during learning.

Network-Specific Modulation of Spike–Field Coherence.

LFPs are considered to represent the aggregate subthreshold activity of neurons in a localized area near the recording electrode (21, 22), thereby providing information about the inputs to an area. A brief description of LFP properties in the orofacial sensorimotor cortex relevant to our analyses is presented in the SI Results and Figs. S9 and S10. Interareal spike–field coherence (MSf and SMf) also exhibited task modulation across multiple frequencies as illustrated for single pairs and for the population of paired signals (Fig. 4 A and B, and Fig. S11). The task modulation of MSf and SMf coherence was distinct from the modulation of LFP spectral power in MIo and SIo (Fig. S10). As was found in MS coherence, the spike–field coherent networks were large and strong in theta and were sparse and weak in alpha/beta/gamma [Fig. S12; McNemar test, P < 0.01; Kruskal–Wallis, peak by frequency, MSf: ${\mathit{x}}_{\left(3\text{,}17991\right)}^2=3558$, P = 0; SMf: ${\mathit{x}}_{\left(3\text{,}19832\right)}^2=2187$, P = 0, post hoc, P < 0.01]. Spike–field networks also exhibited preferred frequency bands; MSf coherence was stronger than SMf in theta, whereas SMf coherence was stronger than MSf in gamma [Kruskal–Wallis, peak by networks, theta: ${\mathit{x}}_{\left(3\text{,}46590\right)}^2=3491$, P = 0; gamma: ${\mathit{x}}_{\left(3\text{,}6170\right)}^2=70$, P = 5e-15; post hoc, P < 0.001]. Theoretic work has suggested that higher firing rates are correlated with stronger spike–field coherence (23). However, differences in firing rates of neurons cannot account for these results as no linear relations were found between firing rates of MIo or SIo neurons and the MSf/SMf coherence in either theta or gamma bands (Fig. S13). Temporal differences between theta and gamma also support the notion of frequency-specific roles of coherence in MIo and SIo. In theta MSf and SMf, the distributions of time of peak coherence were unimodal with peaks around force onset but the distributions were bimodal in gamma with peaks around ±0.35 s relative to force onset (Fig. 4C and Table S1).

Bimodal Distribution of Phase at Peak Coherence Reveals Subnetworks.

The coherence analyses also measure the phase difference, ϕ_xy(t), between two signals x and y. Two signals are in-phase when ϕ_xy(t) = 0°, or antiphase when ϕ_xy(t) = ±180°. The phase at peak coherence (Cϕ) may provide important information on the temporal organization of the coherent signals. For MS coherence, the distribution of Cϕ was bimodal in theta (Fig. 5A, blue, Rayleigh test for bimodal distribution, P < 0.00001; circular mean: –2° and 178°) but uniform in gamma (Fig. 5A, fuchsia, Rayleigh test, P > 0.10). The bimodal distribution of theta Cϕ could not be explained by paired neurons that modulated their firing rates in phase (i.e., both neurons increased or decreased their firing rate concurrently) or out of phase (i.e., one increased its firing rate while the other decreased its firing rate) with each other relative to force generation (SI Results, Fig. S14). A similar bimodal distribution of theta Cϕ was observed in MSf and SMf (Fig. 5B, blue; Rayleigh test for bimodal distribution, P < 0.00001; circular mean: MSf: –4° and 176°; SMf: –6° and 174° at 6 Hz), but the gamma Cϕ distribution in MSf and SMf was unimodal (Fig. 5B, fuchsia, Rayleigh test for unimodal distribution, P < 0.00001; circular mean: MSf: –31°; SMf: –26° at 40 Hz). The bimodal distribution of theta Cϕ was also found in the intraareal coherence MMf and SSf (Rayleigh test for bimodal distribution, P < 0.00001). These results suggest two subnetworks of coherent signals, i.e., in-phase and antiphase. Thus, we examined the theta Cϕ distribution of each neuron with all other neurons in MS to see whether there were distinct subpopulations of neurons. We found that, in MS, a larger number of neurons exhibited a significant bimodal theta Cϕ distribution (Rayleigh test for bimodal distribution, P < 0.05; mean and SEM across D1–D5; Y: 61 ± 13%; B: 40 ± 9%) than unimodal distribution (Rayleigh test for unimodal distribution, P < 0.05; Y: 7 ± 2%; B: 4 ± 1%). We found the reverse for MSf and SMf; over 50% of single neurons in MIo and SIo exhibited a single preferred Cϕ relation with all theta oscillations in MSf and SMf, respectively (Fig. 5C, top plots; Rayleigh test for unimodal distribution, P < 0.05), whereas a smaller number of single neurons showed two preferred Cϕ (Fig. 5C, bottom plots; Y: MSf, 17 ± 5%, SMf, 13 ± 3%; B: MSf, 21 ± 6%, SMf, 19 ± 5%). We verified that phase lags did not depend on specific LFP channels because Cϕ distributions of each LFP channel with all single units were also bimodal (0° and ± 180°). Last, we examined the temporal organization based on the Cϕ by comparing the time of peak coherence of paired signals whose Cϕ fell between –30° to 30° (i.e., the in-phase subnetwork) and between –150° to 150° (i.e., the antiphase subnetwork). In MS, time of peak coherence of the in-phase neurons was significantly earlier than that of the antiphase neurons (Fig. 5D, red; Mann–Whitney, P < 0.001). However, in both SMf and MSf, the antiphase subnetwork had significantly earlier time of peak coherence than the in-phase subnetwork (Fig. 5D, brown and black; P < 0.001). Time of peak coherence of both in-phase and antiphase subnetworks in SMf were significantly earlier than in MSf (Fig. 5D, green and orange; P < 0.001). These results indicate temporal differences based on Cϕ and the network: in-phase MS (n = 4,046) → antiphase SMf (n = 3,842) → antiphase MSf (n = 2,873) → antiphase MS (n = 2,895) and in-phase SMf (n = 2,229) → in-phase MSf (n = 3,217).

Coherent Activity Follows a Spatiotemporal Pattern.

Fig. 3 A and B shows that coherent neuronal pairs do not reach their peak activity all at the same time. Does the spatial location of coherent neurons influence the timing and the magnitude of their peak coherence? We tested this in MS, MSf, and SMf for both theta and gamma coherence but only found significant results in theta (SI Results, Table S2, and Fig. S15). Fig. 6A maps the mediolateral progression of the time of peak theta MSf coherence on the MIo array from D3 of monkey Y; the spatial gradient was apparent from 0.3 s before force onset when MIo neurons located medially exhibited the earliest peak MSf coherence. At 0.4 s after force onset, MIo neurons that exhibited peak coherence were mostly at the lateral border. The mean (across training days) spatial gradient for time of peak theta coherence for MIo neurons in MS and MSf was initiated close to the central sulcus and progressed in the rostrolateral direction (Fig. 6 B and C, MIo array). This may be related to our previous finding that MIo neurons whose spiking activity led the tongue-protrusive force were located closer to the central sulcus than neurons whose spiking activity lagged the force (16). For SIo neurons in MS and SMf, the mean spatial gradient for time of peak theta coherence was from medial to lateral (Fig. 6C, SIo array). Because opening of the lips preceded tongue protrusion, the spatial gradient of the time of peak coherence may relate to the order in which sensory information was received from the lips and the tongue. Indeed, mechanoreceptive fields (RFs) of the lips are located at the medial border of SIo and RFs of the tongue at the lateral border of SIo (24). For the LFP component of MSf and SMf, time of peak coherence progressed medially and toward the central sulcus for MSf and medially for SMf (Fig. S14A, dashed arrows). The results suggest a spatiotemporal organization of coherence based on the temporal relation of MIo neurons’ spiking to tongue force and on the spatial features of MIo and SIo neurons’ RFs.

Coherent Networks Involve Multiple Frequencies.

The simultaneous presence of cortical rhythms in multiple frequencies in the arm sensorimotor cortex has been reported (25, 26), but (to our knowledge) this is the first documentation of frequency-specific coherent activity in the orofacial sensorimotor cortex. Larger and stronger sensorimotor networks occurred in theta, whereas sparser and weaker networks occurred in alpha/beta/gamma, consistent with the general feature of cortical rhythms that slow oscillations engage larger networks whereas fast oscillations engage more localized networks (27). Synchronous activity of local- and large-scale networks may be organized through multiple frequency bands for different cortical processes to operate at multiple temporal and spatial scales. Specifically, different functional demands at different phases of the task may require interareal coherence at different frequencies. Previous studies in the arm sensorimotor cortex suggest that theta synchrony is organized according to movement phases (28, 29). Thus, the high theta coherence at force onset may relate to the generation of tongue-protrusive force. In contrast, the high gamma coherence before and after force onset may relate to other processes such as sensory gating (30, 31), attention (18), and memory encoding/retrieval in association with theta coherence (32).

Reciprocal Interaction Between MIo and SIo.

The observed interactions between MIo and SIo may be explained by a common source of modulatory activity, such as from the thalamus and not due to direct cortico-cortical communication. There are abundant projections from thalamus to MIo and SIo (11, 33), and thalamic neurons have been found to oscillate at 6, 10, and 40 Hz and thus have the potential to generate an oscillatory drive to the cortex in these frequencies (32, 34, 35). However, oscillations may be initiated in the cortex and propagated to the thalamus, which then sends oscillations back to the cortex, thus increasing the cortico-thalamo-cortical resonance (36). This is in agreement with findings in rats, which have spike-and-wave discharges originating from SIo, and then propagating to the thalamus (37). Alternatively, our results may instead represent a direct interaction between MIo and SIo. The dense anatomical connections between MIo and SIo provide a substrate for coherent firing of neurons that may underlie the formation of neuronal assemblies (2, 38, 39). Thus, the increased proportion of coherent neurons may represent the new coupling of a motor output to specific sensory inputs as learning unfolds. Our results also suggest that this process involves reciprocal interactions; MIo modulates sensory processing in SIo (40, 41) and SIo transmits afferent information to MIo critical for successful task performance (5, 6, 42, 43). Because both MIo and SIo have bilateral orofacial representations, cross-hemisphere interactions in the orofacial sensorimotor area may share some properties of the intercortical coherence observed here.

Coherence Phase Reveals Temporal Organization of Subnetworks.

To our knowledge, this study is the first to report a bimodal distribution of phase at peak theta coherence in the sensorimotor cortex, i.e., in-phase (zero lag) and antiphase (near 180° lag). Zero-lag coherence has been suggested to result from either common input, reciprocal interaction, or a combination of both (44), whereas non–zero-lag coherence in gamma has been suggested to reflect interareal conduction delays in unidirectional interactions (45–47). Antiphase synchronization (15–30 Hz) between spikes in the parietal reach region and LFPs in the dorsal premotor cortex has been implicated in down-modulation of communication to selectively prevent transmission of movement-related information (48). Likewise, bimodal distributions of relative phases (0° and 180°) between LFPs (8–25 Hz) in the prefrontal and posterior parietal cortical areas has been suggested to represent specific patterns of coupling among neurons in these areas (49). Given that conduction delays within an area and even across areas are much smaller than the half period of a theta cycle, the antiphase relations in theta coherence found here are likely not due to conduction delays but may represent different populations of MIo neurons whose activity either leads or lags the force generation (16, 49) or a down-modulation of communication between MIo and SIo to prevent transmission of irrelevant sensory or movement-related information (48). The results further refine the understanding of the temporal sequence of coherent activity between MIo and SIo and suggest that cortico-cortical coherence may use phase to organize the activity of neuronal assemblies at different timescales or for segregating feedforward and feedback influences so that orofacial sensorimotor behaviors using different muscles and motor patterns are appropriately coordinated during learning.

Subjects.

All experiments were performed in two adult male rhesus macaques (Macaca mulatta), B (10 kg) and Y (12 kg). All protocols were approved by the University of Chicago Animal Care and Use Committee and complied with the National Institutes of Health Guide for the Care and Use of Laboratory Animals (50).

Behavioral Task.

We used a long-term learning paradigm wherein the subjects were exposed to the same behavioral task parameters over days until subjects achieve a success rate >75% consistently for 3 d (16). Monkeys were trained to protrude the tongue onto a force transducer and apply isometric force at the level cued by target positions. Fig. 1A illustrates the sequence of events in a trial. Detailed description of the task can be found in SI Methods. The behavioral program was written using Spike2 software (Cambridge Electronic Design). Force transducer (Revere Transducers; mode 462-D3-2-10P1R) signals and the behavioral event logs and time stamps were recorded at 2 kHz and stored using a Power 1401 data acquisition system (Cambridge Electronic Design). User-designed pulse signals were generated to mark behavioral events and were sent to the neural data acquisition systems for off-line synchronization of time stamps across the different data acquisition systems.

Further details are available in SI Methods.

Description of the Behavioral Task.

The trial started with the appearance of a cursor that represented the amplitude of the tongue-protrusive force applied on the transducer. After a random period between 0.75 and 1.25 s from trial start, the base target window appeared to cue the monkey to keep the cursor within the base target window by generating force between 1 and 15 g for a random hold period between 0.5 and 1 s. Upon successful hold at the base target, the force target window appeared to signal the monkey to move the cursor into the force target window. The force target was set at 50 g (D1–D3) or 80 g (D4 to D5), and the size of the force target window corresponded to a range of ±30 g from the required force level. To achieve success, monkey had to generate a force of at least 20 g (for 50-g target) or 50 g (for 80-g target) within the allotted time (5 s). Because there was no requirement to hold the force within the force target range, a trial was still successful even when the force exceeded the range. To indicate success, the force target window changed color and the monkey immediately received a juice reward. We set an intertrial interval of 3 s.

Electrophysiology.

Under general anesthesia, each monkey was chronically implanted with two silicon-based arrays of 100 microelectrodes (BlackRock Microsystems), one in MIo and one in SIo of the left hemisphere. The microelectrodes on the array were separated from their immediate neighbors by 400 μm and their length was 1.0 mm for all implanted arrays except for one array that was 1.5 mm in length (MIo of monkey Y). Implantation sites were verified based on surface landmarks and exhibited evoked responses from the tongue and fingers following monopolar surface stimulation of MIo (50 Hz, 200-µs pulse duration, 2–5 mA) during the surgical procedure. Signals from both arrays were amplified with a gain of 5,000, simultaneously recorded digitally (14 bit) with a sampling rate of 30 kHz, and hardware-filtered between 0.3 and 7.5 kHz (Cerebus acquisition systems; BlackRock Microsystems). This was followed by digital bandpass filtering for spikes (0.25–7.5 kHz) and for LFPs (0.3–250 Hz). LFPs were then resampled at 1 kHz. Spike waveforms were stored and sorted off-line using Off-line Sorter (Plexon). Data from array channels with no signal or with large amounts of 60-Hz line noise were excluded.

Data Analysis.

We performed all of the analyses for each dataset as described below. Results were pooled across all datasets (five datasets per monkey) whenever appropriate. Only successful trials were used in the analyses. We included only task-modulated neurons and LFP channels with task-modulated neurons, i.e., neurons whose mean firing rate during the movement epoch (force onset to 0.5 s after) differed significantly from the mean firing rate during a hold period of 0.5 s before target onset until target onset (paired t test, P < 0.01). Of the total number of neurons recorded from both monkeys (MIo: n = 989; SIo: n = 906), 80% (n = 793) of MIo neurons and 72% (n = 650) of SIo neurons were task modulated. LFP channels included in the analyses were 576 for MIo and 528 for SIo (pooled across monkeys and sampled training days). Coherence was analyzed by using the multitapers method of the Chronux Toolbox (51, 52). Coherence, ${\text{C}}_{\text{xy}}$, is a frequency-domain representation of the cross-correlation between two signals, i.e., the raw (unsmoothed) spike trains of neurons x and y for spike–spike coherence, and the raw (unsmoothed) spike trains of neurons x and the LFPs of channel y for spike–field coherence. It is calculated as the cross-spectrum, $S_{xy}$, normalized by the geometric mean of their autospectra, $S_{xx},\hspace{0.17em}{S}_{yy}$, respectively, $C_{xy}=S_{xy}/\sqrt{S_{xx}{S}_{yy}}$, where $C_{xy}$ is a complex number whose modulus corresponds to the degree of coherence (0–1) and the phase as the relative phase difference between the two spike trains. We used a 0.5-s sliding window with 0.01-s steps and applied a time-bandwidth product of 3 and 5 orthogonal Slepian tapers to the first 100 trials. To determine the statistical significance of the coherence, we calculated the coherence between paired signals by shuffling the trials on one of the signals for the time window of the peak coherence. Trial shuffling was tested with 1,000 and 10,000 repetitions which yielded similar results. Since trial shuffling with 1,000 and 10,000 repetitions yielded similar results, we reported the results for trial shuffling with 1,000 repetitions. Paired signals with peak coherence that exceeded the highest 49th coherence amplitude obtained from shuffling were deemed significant (P < 0.05). To address the concern that coherence between pairs of neurons is not independent, we also performed the same analyses on unique pairs (i.e., no two neuronal pairs have a common neuron) from each dataset and found similar results. The unique pairs (i.e., no two neuronal pairs have a common neuron) were determined once by arranging MIo/SIo neurons in ascending order of the electrode number they were recorded from. We first paired, for example, unit MIo-1 to unit SIo-1. Then, we paired the next unit MIo-2 with another SIo neuron that had not been paired previously with a MIo neuron. To evaluate spatial patterns of coherence, we performed a multivariate linear regression with x and y coordinates of the electrode on which a MIo or SIo neuron was recorded as the predictors and the time of peak coherence (or peak coherence magnitude) as the response variable. Note that a single neuron contributes multiple values to the regression because it is coherent with many other neurons. Regression of the time of peak coherence was weighted by the peak coherence. We then derived the orientation of the significant linear gradient by computing the inverse tangent of the coefficients. Circular statistics was performed with CircStat (53). We used the nonparametric Kruskal–Wallis one-way analysis of variance and the Tukey–Kramer test for post hoc multiple paired comparison with significance level set at P < 0.05, unless otherwise noted. All other analyses were performed using built-in and user-defined functions in Matlab (Mathworks). The data will be made available upon request.

Cross-Validation of the Time of Peak Coherence.

We applied a cross-validation approach to verify the results of the shuffling test that determined the paired signals with significant modulation of coherence. The cross-validation involved determining the time of peak coherence on a training subset and then testing the significance of the peak value at that time on a test subset. Specifically, we took a subset (50%) of the paired signals that showed significant coherence and determined their mean time of peak coherence. We then tested the remaining dataset for significance of the coherence at this time of peak coherence using the shuffling method. This method addresses the issue of bias because the time of peak coherence to test for significance was predetermined on the training set. We found significant proportions of paired neurons from a test subset exhibited coherence above chance level at the time of peak coherence determined from a training subset (Figs. S4 and S5).

Local Field Potentials in MIo and SIo.

Fig. S9 A and B illustrates single-trial local field potentials (LFPs) (2–6 Hz) recorded simultaneously from electrodes in MIo and SIo. Spectral analysis of LFPs revealed dominant power in the theta (2–6 Hz) band, but weak power in the alpha (6–13 Hz), beta (15–30 Hz), and gamma (30–50 Hz) (Fig. S9C). Fig. S9D plots within-trial changes in the amplitude of the bandpass-filtered (zero-phase forward and reverse fourth-order Butterworth filter) and trial-averaged LFP signals recorded from a single channel in MIo and another in SIo.

For each LFP channel y and time point j over all trials, we quantified percent phase locking (PPL) of bandpass-filtered LFPs to the force onset (54, 55) as $\mathrm{PPL}\left(y,j\right)=100\left[1-H\left(\phi \left(y,j\right)\right)/H_{\max }\right]$, where $\phi \left(y,j\right)$ is the instantaneous phase and a(y,j) is the instantaneous amplitude of the Hilbert transform of the LFP signal. H is the Shannon entropy $H\left(y,j\right)=-\hspace{0.17em}{\displaystyle {\sum }_{k=1}^N{p}_k{\log }_2\left(p_k\right)}$, where N is the number of bins (n = 4) and p_k is the fraction of values of $\phi \left(y,j\right)$ of the kth bin. $H_{\max }={\log }_2\left(N\right)$. PPL is a measure of the trial-by-trial consistency of phase at a given time and is 0 when phases are randomly distributed and 100 if they are identical. The troughs (±180°) of the LFPs were found to be aligned to the force onset, suggesting phase locking of LFPs to behavioral events (55). Indeed, many single LFP channels in MIo and SIo exhibited phase locking across multiple frequencies (Fig. S9E). Typically in the theta band, SIo showed phase locking closer to force onset than MIo, but LFPs in most channels in MIo and SIo showed maximal phase locking after force onset (Fig. S9F). Such modulation of the amplitude and phase of the ongoing oscillations signifies evoked responses in the LFPs during the performance of the tongue protrusion task.

Bimodal Distribution of Phase at Peak Coherence (Cϕ).

One could argue that the bimodal distribution of MS Cϕ in the theta band might be due to some neuron pairs concurrently increasing or decreasing their firing (in-phase modulation) while other pairs modulate their firing rate out of phase (e.g., one neuron increases its firing rate while the other decreases its firing rate) . To test this, we examined the distribution of theta MS Cϕ obtained from shuffled trials. Shuffling the trials will break any coherence that is beyond the inherent rate modulation of neurons. Assuming that a neuron’s firing pattern (increasing or decreasing) was consistent across trials (because the movements are approximately repeated from trial to trial), the distribution of Cϕ from shuffled data would be bimodal if the coherence were to be explained by in-phase or out-of-phase firing-rate modulation. We found that this was not the case, i.e., the distribution of theta MS Cϕ in shuffled trials were not bimodal (Fig. S14).

Spatial Gradients of Peak Theta Coherence.

We also found significant spatial gradients with peak coherence (Table S2 and Fig. S15B); peak coherence of SIo neurons (SMf) and LFPs (MSf) progressed in the rostral direction and were highest near the central sulcus where SIo neuronal mechanoreceptive fields were smaller and less complex (56). Significant spatial gradients were observed for MIo neurons (MSf) and LFPs (SMf); however, the directions were not consistent across monkeys (Fig. S15B).