Theta Synchrony Is Increased near Neural Populations That Are Active When Initiating Instructed Movement

Abstract Theta oscillations (3–8 Hz) in the human brain have been linked to perception, cognitive control, and spatial memory, but their relation to the motor system is less clear. We tested the hypothesis that theta oscillations coordinate distributed behaviorally relevant neural representations during movement using intracranial electroencephalography (iEEG) recordings from nine patients (n = 490 electrodes) as they performed a simple instructed movement task. Using high frequency activity (HFA; 70–200 Hz) as a marker of local spiking activity, we identified electrodes that were positioned near neural populations that showed increased activity during instruction and movement. We found that theta synchrony was widespread throughout the brain but was increased near regions that showed movement-related increases in neural activity. These results support the view that theta oscillations represent a general property of brain activity that may also play a specific role in coordinating widespread neural activity when initiating voluntary movement.


Introduction
Theta oscillations (3)(4)(5)(6)(7)(8) in the human brain have been linked to a wide range of perceptual and cognitive functions, but their relation to the motor system is less clear (Buzsáki, 2006;Cavanagh and Frank, 2014;Jacobs, 2014;VanRullen, 2016). It has been hypothesized that theta oscillations provide a mechanism for temporally coordinating widespread sensory, goal-related, and motor neural populations that are behaviorally relevant for voluntary movement (Cavanagh et al., 2012).
In support of this view, a recent scalp electroencephalography (EEG) study observed that theta oscillations were phase locked to the initiation of voluntary movement and also related to performance on a visual discrimination task (Tomassini et al., 2017). These data suggest that theta oscillations are related to both sensory and motor behaviors, possibly playing a specific role in coordinating widespread sensory and motor neural activity during voluntary movement. However, because of the limited spatial resolution of scalp EEG, these data do not speak to the anatomic relationship between theta oscillations and brain regions containing neural populations that are behaviorally relevant for movement (e.g., sensory, goal-related, and motor). As an alternative explanation, it is possible that theta oscillations are anatomically widespread and independent from these behaviorally relevant neural populations, suggesting an indirect relationship to voluntary movement, rather than a direct role in coordinating movement-related neural populations We studied whether theta oscillations specifically synchronize regions that are behaviorally relevant for initiating voluntary movement by obtaining intracranial EEG (iEEG) recorded from patients with drug-refractory epilepsy as they performed a simple instructed movement task . high frequency activity (HFA; 70-200 Hz, often labeled "high gamma") from iEEG data provides a measure of local firing activity with high spatial and temporal resolution (Manning et al., 2009;Ray and Maunsell, 2011;Burke et al., 2015;Dubey and Ray, 2019) and can be used to identify electrodes near sensory and motor neural populations (Cogan et al., 2014;Flinker et al., 2015). Moreover, low frequency components of iEEG data, such as theta, can measure widespread network changes that can also be observed at the level of scalp EEG (Buzsáki, 2006;Burke et al., 2013;Solomon et al., 2017;Bickel et al., 2018). Our approach was to use task-related HFA increases to identify electrodes that were in proximity to task-related neural populations and then assess interactions between these regions by measuring theta synchrony between them (Solomon et al., 2017(Solomon et al., , 2019 We grouped electrodes based on their proximity to distributed behaviorally relevant neural populations based on patterns of cue-evoked HFA. These included a widely distributed "instruction-related" group of electrodes that showed HFA increases during instruction presentation and a perirolandic-localized "movement-related" group of electrodes that only showed HFA increases during movement, suggesting proximity to sensory/goal-related neural populations and movement-related neural populations, respectively. We found that theta synchrony occurred between widespread brain regions, including regions that did not show HFA increases, consistent with a general role in brain function. However, theta synchrony was further increased near regions that showed movement-related HFA increases. These results support the view that theta oscillations play a role in coordinating distributed neural activity when initiating voluntary movement.

Subjects
Patients with drug-refractory epilepsy underwent a surgical procedure in which grid, strip, and depth electrodes were implanted to localize epileptogenic regions (Table 1). Clinical circumstances alone determined the number and placement of implanted electrodes. Data were collected from our institution and was approved by the Institutional Review Board. Informed consent was obtained from all participants. In total, we recorded neural activity from nine subjects (ages 21-53; three females, six males). None of these patients had seizures originating from sensory or motor cortex. We did not specifically exclude electrodes based on epileptic activity (see Discussion, Limitations)

Instructed movement task
Subjects were asked to perform a simple instructed motor task at the bedside during the epilepsy monitoring period, ranging from 2 to 14 days after implantation. In each trial, subjects were asked to move either their right hand, left hand, or mouth and tongue (henceforth, "specific movements"). They were presented with written instructions on a laptop screen. Before the task, hand movements were demonstrated as opening and closing the hand and mouth and tongue movements were demonstrated as repetitive movement of the jaw and tongue without producing words or sound. Each trial consisted of three cues presented in sequence: a wait cue ("wait for instructions"), an instruct cue (e.g., "on GO! please move your right hand"), and a move cue ("GO!"; Fig. 1A) Each screen was displayed for 5 s with 24 trials per session. We could not perform behavioral analyses as reaction times were not recorded during the task.

Electrocorticographic recordings
We recorded iEEG from Ad-tech subdural (grids and strips, 4-mm contacts, spaced 10 mm apart) and intraparenchymal depth electrode (1.1-mm diameter, four contacts spaced 5 mm apart). Intraparenchymal depth electrodes were typically placed in medial temporal lobe structures but also used to target frontal lobe structures (subjects 2, 4, and 11). iEEG was recorded using a Nicolet or Natus EEG system. Based on the amplifier and the discretion of the clinical team, signals were sampled at either 250 or 512 Hz. Signals were converted to a bipolar montage by taking the difference of signals between each pair of immediately adjacent electrodes on grid, strip, or depth electrodes. The resulting bipolar signals were treated as new virtual electrodes (henceforth, "electrodes"), originating from the midpoint between each electrode pair (Burke et al., 2013). Analog pulses synchronized the electrophysiological recordings with stimulus presentation events. We excluded electrodes that recorded prominent 60-Hz electrical line noise, defined as electrodes that showed greater spectral power in 58-to 62-Hz range as compared with the 18-to 22-Hz range. Subject 5 underwent a montage change in between sessions resulting in different electrode labels. We only included data from the original montage.

Anatomical localization
Intracranial electrodes were manually identified on each postoperative CT scans. To map electrode coordinates from the CT scan onto the cortical surface, we registered each postoperative CT scan to each patient's preoperative MRI scan using a rigid-body 6 degrees of freedom affine transformation algorithm, and manually adjusted each transform such that electrodes were positioned as close to the cortical volume as possible. We co-registered each patient's preoperative MRI scan to the MNI152 brain to obtain anatomic labels (Jenkinson et al., 2012). Based on MNI152 labels, electrodes were manually assigned to one of several regions: prefrontal, perirolandic, parietal, temporal, medial temporal, or occipital (Fig. 1B). Electrodes that remained unlabeled based on the coregistration to the MNI152 volume were manually assigned to one of these locations.

Extracting spectral power
We extracted 3 s segments of iEEG data from 1000 ms before and 2000 ms after each cue presentation ("wait" cue, the "instruction" cue, and the "Go" cue). We extracted spectral power with 50 complex valued Morlet wavelets (wave number 7) with center frequencies logarithmically spaced from 2 to 200 Hz (Addison, 2002). We first squared and then log-transformed the wavelet convolutions, resulting in a continuous representation of logpower surrounding each cue presentation. We averaged these log-power traces in 500-ms epochs with 400-ms overlap surrounding the presentation of each task related cue (Fig. 1C). We z-transformed power at each frequency by the mean and standard deviation of power values obtained from randomly selected clips of iEEG data recorded from that session so as to not bias values toward any particular task-related event (Burke et al., 2013;Ramayya et al., 2015).

Identifying electrodes that showed task-responsive activity
We identified "task-responsive" electrodes as those that showed significant spectral power changes at any frequency in relation to the instruction or movement cue. For each electrode, we identified spectrally and temporally contiguous power differences between task conditions by performing a cluster-based permutation procedure that accounts for multiple comparisons (Maris and Oostenveld, 2007). We describe the statistical procedure for the comparison between movement and wait intervals ("move-wait"), but separately performed this procedure comparing instruction and wait intervals ("instruct-wait"). As suggested by Maris and Oostenveld, we began by performing an unpaired t test at each time interval comparing power distributions associated with all movement and wait trials performed by the subject. Using an uncorrected p = 0.05 as a threshold, we identified the largest cluster of adjacent time-frequency windows that showed positive t statistics (greater power following movement compared with wait trials), and the largest cluster of adjacent time-frequency windows that showed negative t statistics (greater power during wait trials compared with move trials). By taking the sum within each of these clusters, we computed positive and negative "cluster statistics," respectively. To assess the statistical significance of each cluster statistic, we generated a null distribution of cluster statistics based on 1000 iterations of shuffled data (on each iteration, "move" and "wait" labels were randomly assigned to power values recorded during the session). Based on where each cluster statistic fell relative to the null distribution, we generated a one-tailed p value for each effect, that we converted to a two-tailed p value. For instance, a clustered power increase with a onetailed p value of 0.025 was assigned a two-tailed p value of 0.05, corresponding to a 5% false-positive rate of identifying either a positive or negative cluster at that electrode. We considered an electrode to be task-responsive if we observed a cluster statistic with a p , 0.05 during either the move-wait comparison or the instruct-wait comparison (Fig. 1D).
Grouping electrodes based on cue-evoked HFA We grouped task-responsive electrodes based on HFA changes in relation to the instruction cue and the movement cue. For each electrode, we measured HFA as average z-scored, log-transformed power of wavelets ranging from 70 to 125 Hz (log transform before z-score, see above, Extracting spectral power). We defined a baseline interval as the 500 ms before the wait cue. We performed a paired t test between mean HFA values during the 0-to 1000-ms time interval after the instruction cue and this prewait baseline interval to measure instruction-evoked HFA. Similarly, we performed a paired t test

Wait
Instruction Move z-score z-score t-stat Figure 1. Methods. A, Instructed movement task. Each trial consisted of three cues: a "wait" cue, an "instruction" cue (where subjects were instructed to prepare of three movements: move their right hand, left hand or mouth and tongue), and a "go" cue movement (referred to as "movement-cue" throughout the text), where subjects were asked to carry out the instructed movement. B, We recorded from 490 bipolar electrode pairs (henceforth, "electrodes") throughout neocortex and medial temporal lobe that are shown here in standard MNI152 space (henceforth, "brain plot"). C, We illustrate our methods of extracting neural activity using a three second voltage recording from an example electrode in perirolandic cortex (indicated by a red circle in panel B). We used wavelet convolution of voltage time series to extract time-resolved measures of power and phase at a wide range of frequencies (3-200 Hz). We illustrate three example wavelet convolutions (3, 17, 41 Hz). We focus our analyses on two neural signals: (1) HFA, power in wideband 70-to 125-Hz frequencies, which is an established correlate of aperiodic local spiking activity (dashed red line indicates power estimates); and (2) theta oscillations, periodic membrane potential fluctuations in narrowband 3 to 8 Hz frequencies, that are hypothesized to facilitate interregional neuronal interactions (dashed blue line indicates phase estimates). D, We identified task-responsive electrodes using a non-parametric cluster-based statistical procedure that identified significant task-related power changes in any contiguous time-frequency interval. We illustrate our method to detect movement-related changes in power ("move"-"wait") at an example electrode (same as C). Time-frequency plots summarizing data across various time and frequency windows (3-200 Hz; centered from À500 to 1000 ms relative to cue onset). We show average cue-locked power (left, "wait"; center-left, "move") or cue-locked power differences (center-right, t statistics comparing "move" and "wait" conditions). Power values have been log transformed and normalized within frequency ("whitened"). White line indicates cue onset. We obtained non-parametric p values for each time-frequency cluster by comparing effect size to a null distribution of cluster statistics (right). Vertical blue line marks effect size of low frequency power decrease relative to null distribution, whereas vertical red line marks high frequency power increase. We considered an electrode to be task-responsive if it showed a significant instruction or movement-related power change in any time-frequency range.
between mean HFA values during the 0-to 1000-ms time interval after the movement cue and this prewait baseline interval to measure movement-evoked HFA. We grouped electrodes based on evoked HFA changes as follows.
First, we identified instruction-related electrodes as those that showed increased HFA following the inst ruction cue (t . 2.5, p , 0.05). Second, we identified movement-related electrodes as those that showed increased HFA following the movement cue (t . 2.5, p , 0.05), but that did not show instruction-related HFA increases. Third, we identified "HFA decrease" electrodes as those that showed decreased HFA either during the instruction cue or the movement cue (t , À2.5, p , 0.05). Finally, we labeled all remaining task-responsive electrodes as "HFA null" as they did not show significant cue-related changes in HFA during the task.
We selected these grouping criteria as a method to generally distinguish distinct patterns of local neural activity observed in this dataset, with the acknowledgment that specific boundary criteria are arbitrary. We used this datadriven approach to identify electrode groups rather than a region of interest analysis because some patterns of local neural activity might be widely distributed throughout the brain, and because a particular region might contain distinct and opposing patterns of neural activity . We grouped electrodes that showed both instruction-related and movement-related increases in HFA as part of the instruction-related group, rather than the movement-related group. This response pattern likely reflects a combination of visually responsive neural populations (that show increases in activity following both cues) but may also include preparatory motor populations that show increased activity during movement instruction and execution.

Identifying theta oscillations at each electrode
We identified oscillations at each electrode by assessing whether the power spectrum showed narrowband peaks above the 1/f background activity using a recently described parametric curve-fitting method (Donoghue et al., 2020). This approach has two main advantages over the simpler alternative of averaging power in narrowband frequency ranges. First, it avoids conflating any potential oscillatory components with the aperiodic background component of the power spectrum, which has been shown to reflect asynchronous neural spiking and noise (Manning et al., 2009;Ray and Maunsell, 2011;Dubey and Ray, 2019). Second, this approach does not require predefined frequency ranges to identify oscillatory spectral peaks in the power spectrum and can account for electrode-to-electrode variability in the center frequency of oscillations.
To apply this method, we concatenated iEEG data from each trial during the entire time interval (0-5000 ms following the wait, instruction, and move cues) into a single time series for the entire recording. We computed the power spectrum of this time series using Welch's method for frequencies ranging from 2 to 50 Hz. Briefly, we fit the aperiodic 1/f background component of the power spectrum using an exponential function (in log power vs linear frequency space) and then fit Gaussian peaks to the residual "flattened" power spectrum to assess for oscillatory peaks [python fitting oscillations and one over f (FOOOF) package Donoghue et al., 2020]. We labeled electrodes as recording theta oscillations if we observed a spectral peak with a center frequency in the 3 to 8 Hz range.

Measuring pairwise theta phase synchrony between electrodes
We studied pairwise connectivity between electrode pairs by measuring theta phase synchrony using methods similar to a recent study of medial temporal lobe theta (Solomon et al., 2019).
We used wavelets to extract instantaneous theta phase for each electrode throughout the task. We convolved the iEEG signal from each electrode with complex-valued Morlet wavelets from 3-8 Hz (wave number = 5, 6 wavelets spaced 1 Hz apart). Each wavelet was convolved with 6000 ms of data surrounding the instruction and movement cue (À1000 to 5000 ms surrounding each cue) and buffered with 1000 ms at the beginning and end of each segment (clipped after convolution). We averaged phase values across wavelets (circular mean) resulting in a single theta phase value for each time sample at each electrode.
For a given pair of electrodes, we measured within-trial theta phase coupling across time in 1000-ms time intervals spanning the trial epoch. This duration allows for at least three cycles of a theta oscillation. We focused on two time intervals that we hypothesized would be important for coordinating neural populations for voluntary movement: the 1000 ms surrounding the instruction cue and the move cue. We centered each of these time intervals from 250 ms before the cue to 750 ms following the cue to account for any temporal smearing into the prestimulus interval from the wavelet convolution. We also present data from surrounding time intervals using a sliding window analysis (ranging from À1000 ms to cue, to 1000 to 2000 ms postcue).
For a given time interval, we computed theta phase differences between the electrodes for each sample of time. The null hypothesis states that the distribution of phase differences across trials between two unrelated signals should be uniformly distributed on a unit circle. We assessed the non-uniformity ("tightness") of the distribution of phase differences for each pair at each time interval by computing the resultant vector length (RVL; python circstat package, Circular Statistics Toolbox; Berens and Velasco, 2009). RVL values can range from 0 (uniform circular distribution suggesting independent phases) to 1 (non-uniform distribution suggesting high phase coupling). We obtained a RVL value separately for each trial for the given time interval that we refer to as the "true" RVL distribution in the section below.
We used a non-parametric resampling procedure to estimate the null distribution of RVL values for a given pair of electrodes as follows. For each iteration, we randomly selected a 1000 ms clip of phase values from each electrode that had an intact autocorrelation structure but were mismatched in trial number and in temporal relation to cue presentation. This method ensured that the null distribution would not be influenced by event-related phase reset phenomena occurring at both electrodes. First, we randomly select a trial for each electrode, then circularly shifted each 6000-ms phase clip by a random value, and then selected a contiguous 1000 ms clip for each electrode. We computed a null RVL value by comparing phase differences between these random phase clips. We repeated this procedure 1000 times resulting in a null RVL distribution for a given electrode pair. To measure the extent to which phase coupling was greater than expected by chance during a given time interval, we performed an unpaired t test between true RVL distribution and the null RVL distribution. We refer to the resulting t statistics as "synchrony t statistics" when presenting results, and specifically refer to "instruction synchrony t statistics" and "movement synchrony t statistics" for the instruction-related and movement-related time intervals, respectively. We computed a t statistic to assess whether the distribution of RVL values across trials was greater than expected by chance, rather than the mean RVL value that might be heavily influenced by a handful of outlier trials.

Statistical tests
We performed ANOVA or Student's t tests to compare continuous distributions and x 2 test to compare categorical distributions. We performed false discovery rate (FDR) correction for multiple comparisons (Benjamini and Hochberg, 1995). We considered an FDR-corrected p value , 0.05 to be statistically significant. We also occasionally report uncorrected p values as noted. We performed analyses using MATLAB and Python using both publicly available packages (e.g., NumPy, for numerical computing; SciPy, for statistics and signal processing; MNE, for spectral analyses; pycircstat, for circular statistics; FOOOF, fitting oscillations and one over f, for power spectrum modeling; and statsmodels, for regression modeling) and custom code.
To assess whether task-related neural activity differed in relation to specific movements (left hand vs right hand vs mouth and tongue), we applied a one-way ANOVA on the distribution of mean power within the time-frequency range of the clustered power change at each electrode (uncorrected p , 0.01).

Data sharing
Behavioral and neural data obtained for this study and associated code will be made available on request

Code accessibility
The code described in this paper is freely available online at https://github.com/ashwinramayya/code_RamaEtal20_ sensorymotor.

Identifying task-related electrodes
We obtained iEEG recordings from nine patients as they performed a simple instructed movement task with distinct instruction and movement intervals (Fig. 1A). We recorded from 490 bipolar electrode pairs across widespread brain regions (Fig. 1B). We excluded 70 electrodes that showed prominent electrical line noise.
We found that 288 of the remaining electrodes were task-responsive, in that they showed a spectral power change in relation to the instruction or movement cue, which was more frequent than expected by chance (x 2 statistic = 336, p , 0.001; 23.4 expected based on the 5% false-positive rate). We identified these power changes without specifying a time interval or frequency range using a non-parametric statistical procedure clusterbased procedure (Maris and Oostenveld, 2007). Briefly, this method assessed whether an electrode showed consistent cue-related power changes in any contiguous time-frequency window that were greater in magnitude than expected by chance. We illustrate this method in Figure 1B-D using an example electrode from the perirolandic cortex that showed an increase in HFA and a decrease in wideband low-frequency power during the movement interval as compared with the wait interval (FDR-corrected p's , 0.001)

HFA identified distinct neural response functions throughout the brain
We identified electrodes that were positioned near behaviorally relevant neural activity using task-related HFA changes (70-to 125-Hz power), a known proxy for local neural firing rates (Dubey and Ray, 2019). At each task-responsive electrode, we studied how HFA changed following the instruction and move cues to estimate the response function of nearby neural populations. We illustrate the high spatial and temporal specificity of HFA by showing distinct response functions from four nearby electrodes from an example subject ( Fig. 2A). Two electrodes showed large time-locked HFA increases (blue electrode only during instruction and orange electrode only during movement). The other two electrodes did not show large HFA changes (yellow showed a small magnitude decrease following instruction and movement, whereas purple showed no reliable change).
We grouped electrodes based on cue-related HFA changes to generally distinguish distinct patterns of neural activity observed in this dataset (Fig. 2C-F). We quantified the HFA changes at each electrode by using t tests to compare mean HFA during the 1000 ms following each cue to a baseline interval (500 ms preceding the wait cue), resulting in an instruction-related t statistic and a movement-related t statistic. We grouped electrodes into one of four mutually exclusive groups as follows (HFA electrode groups). First, instruction-related electrodes were those that showed increased HFA following the instruction cue (n = 53 electrodes, 7 subjects) (Fig. 2C). Second, movement-related electrodes were those that showed only increased HFA following the movement cue (n = 34 electrodes, 6 subjects) (Fig. 2D). Third, HFA decrease electrodes were those that showed decreased HFA either following either cue (n = 54 electrodes, 8 subjects) (Fig.  2E). Fourth, HFA null electrodes were task-responsive electrodes that did not show reliable cue-related HFA changes (n = 147 electrodes, 9 subjects) (Fig. 2F).
HFA electrode groups showed distinct anatomic distributions (Fig. 3). Movement-related electrodes (Fig. 3B) were heavily clustered in perirolandic regions (n = 21/34, x 2 statistic = 19.8, p , 0.001), whereas instruction-related (Fig. 3A) and HFA null electrodes (Fig. 3D) were widely distributed throughout the brain (x 2 test, ps . 0.15). HFA decrease electrodes (Fig. 3C) were also widespread but were more frequently observed in parietal and temporal regions than expected (x 2 statistic = 12.3, p = 0.012). distinct HFA responses during the task from the same subject. Brain plot illustrates the location of each electrode (purple, blue, and orange are in perirolandic, whereas yellow is parietal, on the border between postcentral gyrus and superior temporal sulcus). The HFA responses in each of these electrodes correspond to the average response functions described below. B, Distribution of instruction and movement-related HFA across all electrodes (see main text for details). We identified four groups of electrodes as follows: (1) Instruction-related electrodes (blue, 53 electrodes, 7 subjects), (2) Movement-related electrodes (red, 34 electrodes, 6 subjects), (3) HFA decrease (yellow, 54 electrodes, 8 subjects), and (4) HFA null (purple, 147 electrodes, 9 subjects). Gray error ellipse marks 95% confidence interval of bivariate distribution. C-F, Average HFA (z-scores) at each electrode group and their anatomic distributions on brain plots.

Theta oscillations are widespread and independent of HFA
We found that cue-related HFA changes at each electrode were independent of narrowband theta power (3-8 Hz) after accounting of wideband frequency components of the power spectrum ("low-frequency activity," 2to 30-Hz power) using multivariate linear regression (Fig.  4). At each electrode, we measured instruction-related and movement-related measures of low-frequency activity and theta power using the same time intervals used to measure cue-related HFA (Materials and Methods). We found that movement-related changes in HFA were negatively related to low-frequency activity (t statistic of b coefficient = À3.68, FDR-corrected p = 0.002), but independent of theta (p . 0.5). Instruction-related HFA changes were independent of both low-frequency activity and theta (ps . 0.25).
In a separate analysis, we identified electrodes that showed theta oscillations using a curve fitting procedure to detect 3-to 8-Hz peaks in the power spectrum beyond the background 1/f component (wideband low-frequency changes in the analysis above; Donoghue et al., 2020). We observed good fits to the power spectra across taskresponsive electrodes (R 2 = 0.99 6 0.012; mean 6 SD; Fig. 5A) and observed theta oscillations in 134/288 taskresponsive electrodes (Fig. 5B). These oscillations were widespread across HFA electrode groups and brain regions (x 2 statistic, ps . 0.5; Fig. 5C,D). We focus the remainder of the analyses studying interactions among these electrodes that recorded theta oscillations.

Pairwise theta synchrony is widespread but dependent on movement-related HFA
We observed widespread pairwise phase synchrony between task-responsive electrodes that showed theta oscillations. For each pair of electrodes, we measured the extent to which theta oscillations were "coupled" over time by testing whether theta phase differences were more consistent than expected by chance. We focused on two time intervals that we hypothesized as potentially behaviorally relevant: 250 ms before, and 750 ms following, the instruction and movement cues. In each time interval, we quantified the extent to which phase coupling was greater than expected by chance, resulting in an instruction synchrony t statistic and movement synchrony t statistic. We illustrate our method for quantifying pairwise theta synchrony using three example electrodes from perirolandic cortex that were all task-responsive and showed theta oscillations (Fig. 6). We only included taskresponsive electrodes with simultaneous iEEG recordings that both showed theta oscillations (n = 1806 electrode pairs) We found that average pairwise theta synchrony was greater than expected by chance during both the instruction and movement intervals (paired t tests on distributions of average synchrony t statistics for each subject; instruction: t (8) = 4.42, p = 0.002; movement: t (8) = 2.5, p = 0.037). We observed significant (p , 0.05) pairwise instruction-interval synchrony in 213/1806 electrode pairs and movementinterval synchrony in 253/1806 electrode pairs (x 2 tests, ps , 0.001; n = 90 expected by chance during each interval). We also found a positive correlation between theta synchrony during instruction and movement (Pearson r across pairs, r = 0.65, p , 0.001; and across subjects, r = 0.71, p = 0.029), suggesting that theta synchrony was generally stable throughout each trial and not driven by event-related phenomena.
We found that theta synchrony was heavily influenced by the physical distance between electrodes in a pair and the similarity in the frequency of theta oscillations detected at each electrode. In both cases, we observed a supra-linear increase in pairwise synchrony when electrodes were closer to each other in physical space and theta frequency (Fig. 7). We observed linear relations between pairwise synchrony and log transforms of each of these measures (Pearson's r . 0.25, p , 0.001). Figure 5. Theta oscillations are widespread throughout the brain. A, Example power spectrum from a single electrode illustrating our method to identify oscillations. We applied a curve-fitting algorithm to model each power spectra as a combination of 1/f background (dashed blue line) and oscillatory peaks (black line shows original spectrum, red line shows model fit; see main text for details). Vertical green lines indicate the center frequency of detected spectral peaks (7 and 14 Hz) in this case. B, Distribution of lowest frequency spectral peaks across all electrodes. We identified electrodes recording u oscillations as those that showed a spectral peak with a center frequency between 3 and 8 Hz. Electrodes without a spectral peak are not shown. C, D, We show the distribution of these electrodes across HFA electrode groups (C) and various brain regions (D). We found that theta synchrony was positively related to movement-related HFA increases after accounting for the effects of physical distance and theta similarity using multivariate linear regression. We fit two separate linear models for the instruction and movement time intervals. During the instruction interval, we asked whether instruction synchrony t statistics were dependent on mean instruction-related HFA observed at the electrode pair, whereas during the movement interval, we studied the relation between movement synchrony t statistics and mean movement-related HFA. In both models, we included log-transformed physical distance and log-transformed frequency similarity as additional independent variables. . Theta oscillations are more synchronous than expected by chance. A, We illustrate our method for measuring theta phase synchrony using three example electrodes from the perirolandic cortex (color scheme matched to Fig. 2). Brain plot showing locations of each electrode and power spectra with green lines marking theta oscillation center frequencies (6.99 Hz in electrode 1, 7.07 Hz in electrode 2, and 4.84 Hz in electrode 3). B, Single trial data raw iEEG data (top, solid lines) and filtered 3-to 8-Hz iEEG (bottom, dashed lines) from electrode 1 (blue) and electrode 2 (red) during the instruction interval of interest (250 ms before and 750 ms after cue presentation). Polar plots show distribution of theta phase differences during this time interval showing a clustering around p. Orange line shows RVL (drawn on unit circle), which measures the non-uniformity or (tightness) of the phase distribution (RVL = 0.9 for this time interval). C, Null distribution of RVL values generated from random resampling of contiguous 1-s phase data from electrodes 1 and 2. Vertical red lines mark true RVL values from each trial. Blue vertical line marks the mean true RVL distribution. We quantified the extent to which phase synchrony across trials differed from the null distribution via unpaired t test resulting in a synchrony t statistic for the time interval. D, We show synchrony t statistics for the interaction between electrodes 1 and 2 surrounding several 1-s time intervals surrounding the instruction cue (À1000 ms before 2000 ms after using a sliding window analysis). Horizontal axis tick labels indicate the mean of each time window. Red line marks the mean of the instruction interval of interest (À250 to 750 ms). E, F, Same as C, D for interactions between electrode 1 and electrode 3. We observed significant synchrony between electrode 1 and electrode 2 (synchrony t statistic = 14.9, p , 0.001), but not between electrodes 1 and 3 during the instruction interval of interest (synchrony t statistic = 1.4, p . 0.15). G, H, Mean of synchrony t statistics over time across all electrode pairs (n = 1807) during the instruction (G) and movement intervals (H), with red lines marking the instruction and movement intervals of interest. We observed greater synchrony across electrode pairs than expected by chance (see main text for statistics). Width indicates SEM across electrode pairs. I, We observed positive correlation between instruction synchrony t statistics and movement-related synchrony t statistics across all electrode pairs.

A B C D
Theta synchrony showed an independent, significant positive relation to cue-related HFA changes during the movement interval (t statistic of b coefficient = 2.10, FDRcorrected p = 0.041), but not during the instruction interval (p . 0.15). During both time intervals, theta synchrony was heavily dependent on physical proximity (t statistic of b coefficients . 11.1, FDR-corrected ps , 0.001) and theta frequency similarity (t statistic of b coefficients . 7.3, FDR-corrected ps , 0.001). In a post hoc analysis, we found that instruction-related synchrony also showed a weak positive relation with movement-related HFA (t statistic of b coefficient = 1.97, uncorrected p = 0.049).
Consistent with these results, movement-related electrode group (HFA move) showed increased within-group theta synchrony (between two movement-related electrodes) and out-of-group theta synchrony (movement-related electrode vs non-movement-related electrode) relative to the remaining three HFA electrode groups during both the instruction (one-way ANOVA by group, F statistic . 4.9, ps , 0.003, post hoc t test movement-related electrodes vs all other groups ts . 2.5, ps , 0.012; Fig. 8) and movement intervals (one-way ANOVA by group, F statistic . 6.6, ps , 0.001, post hoc t test movement-related electrodes vs all other groups ts . 3.43, ps , 0.001, all ps FDR-corrected; Fig. 8).

Discussion
We studied iEEG recordings obtained as patients performed a simple instructed movement task. Our main goal was to test the hypothesis that theta oscillations Figure 7. Theta synchrony is heavily influenced by physical distance and similarity in oscillation frequency. A, D, Distributions of Euclidian distance in MNI coordinate space (mm, A) and differences between center frequency of theta oscillations (D) across electrode pairs. We show the relation between synchrony during the instruction interval and these measures but observed similar relations during the movement interval (see main text). B, E, Synchrony showed supra linear increases as a function of decreasing physical distance (B) and decreasing differences in theta center frequencies (E) that were fit by power law functions. C, F, Synchrony showed linear relations to log transforms of physical distance (C) and differences in theta frequency (F). Dashed red line shows line of best fit in each case. synchronize regions containing neural populations that are behaviorally relevant for voluntary movement. We used HFA (70-to 200-Hz power), a known proxy of local firing rates, to identify electrodes that were positioned near neural populations that showed temporally specific task-related neural activity (Manning et al., 2009;Ray and Maunsell, 2011;Dubey and Ray, 2019), and studied their anatomic relation to synchronous theta oscillations. We found that HFA identified heterogenous neural response functions throughout the brain, including distinct responses from nearby brain regions. These results suggest that HFA is local neural signal akin to measuring multiunit activity from neural populations near the electrode (Dubey and Ray, 2019). We grouped electrodes based on their response functions to study distributed neural populations that shared activity patterns. We observed two prominent patterns of HFA increases. First, we identified instruction-related group that showed temporally specific increases in activity following the instruction cue, and second, a movement-related group that showed temporally specific increases in activity during movement. Instruction-related electrodes were widely distributed across regions involved in the ventral visual stream (occipital, temporal) and goal-directed and movement-planning networks (prefrontal and perirolandic regions), suggesting a relation to heterogeneous functions that occur during instruction presentation (sensory processing, goal-related and movement planning. On the other hand, movement-related electrodes were heavily clustered in perirolandic regions, suggest a prominent sampling of movement-generating neural processes We most frequently observed neural activity that distinguished between specific movements (left hand vs right hand vs mouth/tongue) at movement-related electrodes, consistent with a role in selecting specific movements, consistent with prior studies (Crone et al., 1998(Crone et al., , 2006Miller et al., 2007;Schalk et al., 2008;Cogan et al., 2014;Korzeniewska et al., 2015). We rarely observed neural activity that distinguished between specific instructions, which is consistent with recent findings in non-human primates suggesting that population-level neural representations of specific movements only emerge during movement generation (Kaufman et al., 2014;Elsayed et al., 2016).
A challenge in measuring theta oscillations is that narrowband theta power can conflate periodic theta oscillations and aperiodic low frequency power changes that are thought to reflect a distinct underlying neural process Donoghue et al., 2020). Consistent with this view, we found, via multivariate regression, that movement-related HFA changes were related to wideband low frequency power decreases, but independent of theta power. These results are consistent with a power spectrum "tilt" that has been widely observed in iEEG studies (Miller et al., 2007;Burke et al., 2015;Solomon et al., 2017), and may be related to "b " desynchronization that has been observed in scalp and intracranial EEG studies (Murthy and Fetz, 1992;Crone et al., 1998Crone et al., , 2006Hari et al., 1998;Ohara, 2002;Jenkinson and Brown, 2011). However, instruction-related HFA changes were not related to wideband low frequency power, suggesting that HFA increases are not always accompanied by a power spectrum tilt. Further studies are needed to study the behavioral and structural underpinnings of spectral tilt in relation to HFA Gao et al., 2017;Herweg et al., 2020).
We observed theta oscillations at widespread electrodes throughout the brain, independent of local HFA changes following instruction or movement. We measured theta oscillations at each electrode by assessing whether the power spectrum contained narrowband peaks beyond the background 1/f shape using a recently described curve-fitting algorithm (Donoghue et al., 2020). This method allowed us to identify electrodes that showed periodic theta oscillations beyond asynchronous low frequency power changes. These results suggest that theta oscillations reflect a global neural signal, in contrast to HFA that measures local neural activity.
We found that these theta oscillations showed more phase synchrony over time than expected by chance. We measured phase synchrony over time in 1-s intervals throughout the trial rather than across trials in specific time windows to measure ongoing oscillations rather than event-related phenomena (Cohen and Donner, 2013; although see David et al., 2006). We found that theta synchrony was largely stable throughout each trial suggestive of ongoing oscillations, and in contrast to HFA that showed prominent within-trial dynamics.
We focused on phase-phase relations rather than other connectivity measures such as Granger causality or spectral coherence to mitigate the influence of asynchronous power correlations on our analyses (Lachaux et al., 1999;Herweg et al., 2020) and to relate our findings to previous studies of theta synchrony (Lega et al., 2012;Burke et al., 2013;Solomon et al., 2017Solomon et al., , 2019Donoghue et al., 2020). In keeping with recent literature, we define "synchrony" between two signals to imply any periodically coupled temporal relationship, but not to imply a phase difference of zero, as would be expected perfectly coupled oscillators (Mirollo and Strogatz, 1990). Instead, our definition allows for variable offset phase differences as might be expected from multiple uncoupled oscillators or traveling waves (Zhang et al., 2018).
Theta synchrony showed a linear increase as a function of logarithmic decreases in physical distance which is consistent with previous findings (Lachaux et al., 1999) and may suggest a small-world structure to network interactions (Buzsáki, 2006;Bassett et al., 2018). Taken together with the global distribution of theta oscillations, and slow fluctuations across trials, these data suggest that theta synchrony may reflect dynamic functional connectivity between brain regions (Bickel et al., 2018), rather than event related neural activity. We also observed increased theta synchrony between electrodes that showed similar theta oscillation frequencies, which is consistent with the view that multiple oscillations are multiplexed within the theta frequency range (Jacobs, 2014).
After accounting for the effects described above, we found increased theta synchrony involving movementrelated electrodes (that showed movement-related increases in HFA). This result supports the view that synchrony between theta oscillations plays a role in facilitating interactions between widespread behaviorally relevant neural populations during action selection (Cavanagh and Frank, 2014;Herweg et al., 2020). Our results provide an important complement to recent scalp EEG data showing a behavioral link between theta oscillations and sensory and motor functions (Tomassini et al., 2017). By showing that theta oscillations synchronize regions containing neural populations that are active when initiating instructed movement, our data provide an anatomic link between theta oscillations and movement-related neural populations in the human brain. Our data build on a large body of body of literature linking theta oscillations to human behavior in various domains, including perception (VanRullen, 2016), attention (Helfrich et al., 2018;VanRullen, 2018), spatial navigation (Jacobs, 2014), memory (Herweg et al., 2020), and decision-making (Cavanagh and Frank, 2014), and disease (Cavanagh et al., 2017;Singh et al., 2018).
In conclusion, we studied intracranial neural recordings patients with drug-refractory epilepsy performed a simple instructed motor task. We found that HFA measured distinct neural responses from nearby neural populations, suggesting a local signal. In contrast, theta oscillations were widespread and synchronous, suggesting a global neural signal. Theta synchrony was increased near neural populations that showed movement-related increases in local activity, suggesting that theta oscillations coordinate distributed neural representations during action selection.

Future directions
Future studies should assess how theta synchrony fluctuates over time in relation to bottom-up arousal systems (Joshi et al., 2016;Stitt et al., 2018) and top-down cognitive control signals (Cavanagh et al., 2012) during higher cognitive functions such as learning (Montague et al., 1996;Ramayya et al., 2015), decision (Gold and Shadlen, 2007) and memory (Ratcliff, 1978), and assess mechanistic interactions to local neural populations via phase amplitude coupling (Galifianakis et al., 2013;Helfrich et al., 2018). Additionally, the application of graph-theoretic methods may be useful in studying global changes in theta synchrony in relation to cognitive states (Bassett et al., 2018).

Limitations
We did not measure reaction times during the task and have limited ability to relate these neural signals to specific behavior. We had sparse electrode coverage across the nine subjects and have a limited ability to make claims about specific anatomic correlations of the observed neural responses beyond the region of interest analysis performed in the study. The opportunity to obtain human intracranial neurophysiology requires studying patient populations that may show systematic differences in neural structure and function relative to healthy individuals. In this case, increased functional connectivity from epilepsy (Bettus et al., 2008) may contribute to overestimating theta synchrony. However, it is unlikely to explain the HFA results, or the relation between theta synchrony and movement-related HFA.