Theta oscillations coordinate grid-like representations between ventromedial prefrontal and entorhinal cortex.

Grid cells and theta oscillations are fundamental constituents of the brain’s navigation system and have been described in the entorhinal cortex (EC). Recent fMRI studies reveal that the ventromedial prefrontal cortex (vmPFC) contains grid-like representations. However, the neural mechanisms underlying human vmPFC grid-like representations and their interactions with EC grid activity have remained unknown. We conducted intracranial electroencephalography (iEEG) recordings from epilepsy patients during a virtual spatial navigation task. Oscillatory theta power in the vmPFC exhibited a sixfold rotational symmetry that was coordinated with grid-like representations in the EC. We found that synchronous theta oscillations occurred between these regions that predicted navigational performance. Analysis of information transfer revealed a unidirectional signal from vmPFC to EC during memory retrieval. Together, this study provides insights into the previously unknown neural signature and functional role of grid-like representations outside the EC and their synchronization with the entorhinal grid during human spatial navigation.

INTRODUCTION

It has been proposed that mammals represent the world in the form of cognitive maps, which allow them to generate flexible and goal-directed behavior (). Grid cells in the rodent entorhinal cortex (EC) provide a periodic representation of self-location () and constitute a major cell type of the cognitive map. Grid cells have been directly recorded in the human brain (, ). In addition, several models have proposed how the activity of grid cells may translate into “grid-like representations,” i.e., network patterns with sixfold rotational symmetry that can be observed via functional magnetic resonance imaging (fMRI) in the human EC (, ). Recently, we () and others () also described grid-like representations in entorhinal theta oscillations recorded via intracranial electroencephalography (iEEG).

Human fMRI allows investigating the entire brain simultaneously and revealed grid-like representations outside the EC, particularly in the ventromedial prefrontal cortex (vmPFC), which may convey relevant information beyond physical space, including visual space, odor space, and abstract conceptual spaces (–). A recent primate study showed that multidimensional representations of options in a decision-making space in vmPFC (near Brodmann area 32) were based on a coding scheme resembling grid cells (). In addition, a growing body of evidence indicates that both vmPFC and EC are involved in the construction of cognitive maps during various tasks involving spatial navigation, episodic memory, and reasoning (–). Consistent with these findings, theoretical models put forward the ideas that grid cells may constitute a low-dimensional basis for cognitive maps () and that they are crucial for path integration and vector-based navigation (–). Despite growing evidence about vmPFC and EC grid-like representations, the electrophysiological basis of grid representations in vmPFC and their relationship to grid representations in EC have remained unknown.

Several studies have reported that theta oscillations in vmPFC support a range of cognitive functions including memory, imagination, and value processing (–). vmPFC theta oscillations have been observed with magnetoencephalography (MEG) in human navigation experiments (, ). We therefore hypothesized that theta oscillations in vmPFC may organize spatial maps in a grid-like manner. Furthermore, previous theoretical frameworks (, ) and experimental studies (, ) suggest that grid-like representations may be coherent across brain regions. Specifically, it has been proposed that vmPFC and medial temporal lobe (MTL) are cooperatively involved in the encoding and retrieval of contextual memory, whereby the vmPFC exerts a top-down control to select context-relevant information (–). In addition, vmPFC plays a critical role in the formation and updating of cognitive scenarios: By recruiting schematic knowledge from the neocortex and sending it back to MTL, vmPFC-MTL interactions mediate the schema-based retrieval and monitoring necessary for updating of mental scenarios (). Interregional oscillatory synchronization, particularly in the theta band, is considered as a basic mechanism of information transfer in the brain (–). We thus hypothesized that grid-like representations in vmPFC and EC may be coordinated by synchronized theta oscillations. Here, we addressed these questions by reanalyzing iEEG data recorded in epilepsy patients during a virtual spatial navigation task (). Our results indicate that a neural circuit involving vmPFC and EC organizes spatial representations into theta-based grid-like representations to support flexible navigation.

RESULTS

Behavioral results

We recorded from epilepsy patients (N = 15) with intracranial electrodes who conducted a virtual reality navigation task (Fig. 1, A and B, and table S1). The task required participants to repeatedly retrieve the location of eight target objects. During each trial, participants were presented with a target object and were asked to navigate to the location of the object based on memory. They then received feedback using a value signal with five levels and relearned the correct location of the object. Participants conducted the task for at least 40 min. For each subject, we divided the data into six sessions of equal length for later analysis (table S2). The average drop error of patients was 2362 ± 584 (mean ± SD) virtual units (vu), which was significantly above chance level (patient-wise permutation tests indicating random-level performance of 3116 ± 231 vu; paired t test: t14 = −5.39, P < 0.001; Fig. 1C). Drop errors were significantly lower during the last three sessions as compared to the first three sessions (paired t test: t14 = 7.09, P < 0.01; Fig. 1D). These results indicate that the patients could successfully learn the locations of target objects. Participants’ movement directions were evenly distributed in 360° space (Rayleigh’s tests: z < 0.05, P > 0.9; Fig. 1E) and in 60° space (Rayleigh’s tests: z < 0.93, P > 0.49), ruling out a possible confound of unbalanced movement directions during the analysis of grid-like representations (Fig. 1F).

Fig. 1.

Task, behavioral results, and hypothesis.

(A) Left: Object-location memory task. Right: Calculation of drop error in each trial as the Euclidean distance between the response location (green dot) and the correct object location (red dot). (B) Patients received feedback via one of five smileys that were color-coded according to the drop error. For example, light green smileys were shown for drop errors within a radius of 700 to 1500 vu. We divided the trials into three conditions: good (dark and light green), middle (yellow), and bad (dark and light red). (C) Histogram of drop errors across all trials and all patients. Red line, average chance performance across patients. (D) Change in mean drop error between the first three and last three sessions. Gray lines, data from individual patients. **P < 0.01. (E) Sampling per 30° bin of movement direction in 360° space. Error bars indicate SEM across participants. Gray dots indicate value from individual participants. (F) Left: Schematic depiction of firing fields for one grid cell. Right: Schematic depiction of a sixfold rotationally symmetric modulation of oscillatory power by movement direction.

Grid-like representations in prefrontal theta oscillations

To test for possible grid-like representations in the PFC, we considered all electrodes (459 contacts, mean number of electrodes per subject and SEM across subjects: mean ± SEM = 30.6±5.6; fig. S1A) that were located in orbitofrontal and medial prefrontal areas according to a published atlas (fig. S1B) based on multimodal human MRI data (). Using an automatic algorithm (MODAL) for the detection of narrowband oscillations (), we found notable oscillations in the high theta (6 to 9 Hz) frequency range in these regions (fig. S1, D to H). We thus filtered the iEEG signals in this band. We used the same procedure as in our previous iEEG study to quantify a sixfold modulation of iEEG power as a function of movement direction (Fig. 1F and Materials and Methods). Previous studies, including our studies (, , ), found that grid-like signals in EC were present only at fast speed level. Thus, we focused on these movement epochs [table S2; see fig. S6 (A and B) for grid-like representations at the middle and slow speed level]. For each electrode contact, one-half of the data were used as a training set to estimate the preferred grid orientation ϕ in different frequency bands; the other half served as the test set [using cos(6(θ-ϕ)) as a regressor] to examine whether the power in a given frequency band was higher when the direction of movement was aligned with ϕ than when it was misaligned with ϕ. The regression coefficient β quantifies the strength of the grid-like representation in each electrode, whose significance was then evaluated using a permutation test (transformed into Zβ; Materials and Methods).

To identify which subregions of the mPFC showed sixfold rotational symmetry, we applied a cluster-based method developed for iEEG for multiple comparisons on the β values obtained above (, ). We found a significant cluster in vmPFC (area 32 in the atlas; permutation test, P < 0.01; Fig. 2A and fig. S1C; Materials and Methods). This region was highly overlapping with results in previous fMRI studies on grid-like representations (, ). We then selected all electrode contacts in area 32 to illustrate the effect (88 contacts from 15 subjects; Fig. 2A). If a subject had multiple contacts, their β (Zβ) values were averaged first. At the group level, we tested for significant grid-like representations in the high theta band and three other control frequency bands: low theta (2 to 5 Hz), beta (12 to 30 Hz), and low gamma (30 to 60 Hz). We found that only the high theta band was significant (t test: t14 = 7.09, P = 0.028 after Bonferroni correction for multiple comparisons; Fig. 2, B and C). The sixfold modulation of theta power was apparent in 360° space, ruling out that it was only driven by higher power in one particular direction (Fig. 2D). Similar to previous studies, the preferred orientation was not clustered across subjects (Rayleigh’s tests: z < 0.14, P = 0.88; Fig. 2G), but it was clustered across electrodes within subjects (one-sample t test: t14 = −2.14, P = 0.02; fig. S2A). Example electrodes exhibiting substantial sixfold modulation of theta power are shown in fig. S2D. Next, we verified that the results were specific to a sixfold symmetry. We used the same analysis procedure for other rotational symmetries, and the results were not significant (t test: all t14 < 1.93, Pcorr > 0.35; Fig. 2E). We then checked whether the sixfold modulation was specific to the vmPFC. Using the electrodes from two neighboring subregions [dorsomedial PFC (dmPFC): areas 8m and 9 of the atlas; orbital frontal cortex (OFC): areas 11m and 14m of the atlas], we found that the results were not significant (t test: all t14 < 1.42, P > 0.18; Fig. 2F), consistent with a recent study (). Although we did not observe a significant correlation between sixfold modulation and drop error across subjects, we did find that the sixfold modulation in later sessions (with a lower drop error) was significantly higher than zero (t test: t14 = 2.27, P = 0.04), while this was not the case in earlier sessions (t test: t14 = 1.51, P = 0.15) (fig. S2C). See fig. S7 for distribution of sixfold modulation over the whole brain.

Fig. 2.

Medial prefrontal theta power exhibits sixfold modulation by movement direction.

(A) Top: Grid-like representations cluster obtained in the high theta band displayed in MNI standard space. The color bar denotes cluster-corrected t values (permutation test, P < 0.01). Bottom: Depiction of all electrode contacts in vmPFC and EC (black circles) relative to two control regions (dmPFC and OFC). A, anterior; P, posterior; V, ventral; D, dorsal. (B) Raw β scores for sixfold modulation of vmPFC oscillatory power in the high theta band. Black dots, individual subjects. *P < 0.05 after Bonferroni correction for four bands. (C) Mean Zβ scores across subjects with electrodes in vmPFC, separately computed for four different frequency bands: low theta (2 to 5 Hz), high theta (6 to 9 Hz), beta (12 to 30 Hz), and low gamma (30 to 60 Hz). *P < 0.05 after Bonferroni correction for four bands. (D) Theta power was higher during movements aligned with the grid axes as compared to misaligned movements. Purple, aligned; gray, misaligned. (E) Mean subject Zβ scores for different symmetry patterns. *P < 0.05 after Bonferroni correction for fivefold. (F) Mean subject Zβ in control regions. Error bars in (B) to (F) indicate SEM across participants. Gray dots in (B) to (F) indicate value from individual participants. (G) Distribution of mean grid orientations across subjects does not reveal a significant clustering (Rayleigh test: P = 0.88). n.s., nonsignificant.

Coordination of grid-like representations between vmPFC and EC

In our previous study, EC theta oscillations (at 4 to 8 Hz) showed a sixfold modulation of theta power as a function of movement direction. A natural question is thus whether grid-like representations in vmPFC and EC are coordinated. Hence, we focused on subjects who had contacts within both brain regions (N = 7; table S1). We first confirmed that these seven subjects’ EC showed sixfold modulation (one-sample t test of β values against zero: t6 = 3.57, P = 0.04 after Bonferroni correction for multiple comparisons by four bands; fig. S3). We then adapted analysis steps similar to previous fMRI studies to explore the relationship of these two types of grid-like representations (). For each subject, we selected all contacts in vmPFC and EC and paired them with each other (114 pairs in total). For each pair, we used the vmPFC contact to estimate the grid orientation and the EC contact to test the grid orientation, resulting in one regression coefficient βvmPFC->EC for each pair (Materials and Methods). βvmPFC->EC values were averaged across electrode pairs from the same subject. These averaged values were z-scored (relative to surrogate βvmPFC->EC values) for each subject and fed into group-level statistics across the seven patients. We found that ZβvmPFC->EC was significant (permutation test, P < 0.01; Fig. 3A). Again, this effect only existed for a sixfold symmetry, but not for other symmetries (permutation test, all P > 0.5). Next, to examine the specific pattern of vmPFC->EC sixfold modulation, we separated all aligned and misaligned trajectories based on the preferred directions in vmPFC and constructed 12 regressors (Materials and Methods). We found that the theta power of EC was higher in the aligned than the misaligned direction (Fig. 3B). To corroborate this effect, we conducted the same analysis in the reverse direction (i.e., estimating the grid orientation based on theta power from an EC contact and testing this grid orientation based on theta power from a vmPFC contact, resulting in the regression coefficient βEC->vmPFC). Similarly, we found that ZβEC->vmPFC was also significant only for sixfold modulation (permutation test, P < 0.01; fig. S4B). The β values of the two directions were significantly correlated, both across all electrode contact pairs (r = 0.70, P < 0.001) and across subjects (r = 0.91, P = 0.002) (Fig. 3C). Across electrode pairs, vmPFC and EC had similar grid orientations (circular correlation = 0.69; P < 0.001; Fig. 3D). We observed this effect also at the subject level (circular correlation = 0.72; permutation test, P = 0.03; Fig. 3D). These results indicate that the grid systems of the two brain regions share a similar preferred angle, consistent with results of previous fMRI studies (, ).

Fig. 3.

Coordinated grid-like representations between vmPFC and EC.

(A) Z-scored βvmPFC->EC scores across subjects show that theta oscillations in the EC are aligned to preferred grid orientations in the vmPFC. This effect was not observed for four-, five-, seven-, or eightfold symmetry. **P < 0.01. (B) Separate depiction for all aligned (blue) and misaligned (gray) trajectories. Error bars in (A) and (B) indicate SEM across participants. Gray dots in (A) and (B) indicate value from individual participants. (C) Pair-wise z-scored β values for all electrodes in vmPFC and EC in both directions are highly correlated. Inset: Interindividual correlation over subject-wise averages of grid cell–like representations in the two analysis directions. (D) vmPFC and EC show similar grid orientations both across pairs and subjects.

Coordinated grid-like representations are associated with interregional theta synchronization

Next, we tested whether the coordination of grid-like representations was associated with functional connectivity between vmPFC and EC and with behavioral performance. To calculate interregional functional connectivity, we selected continuous fast speed movements with a duration of more than 1 s from the retrieval phase and treated them as different epochs [Materials and Methods; Fig. 4A and fig. S4A; see fig. S6 (C and D) for interregional coherence at middle and slow speed]. For each pair of electrodes and each epoch, we calculated the normalized imaginary coherence (Z-Icoh; Materials and Methods) to eliminate possible confounding effects of volume conduction (). We then averaged the results over all epochs to obtain a coherence profile and compared it to zero across participants. This analysis revealed significant coherence at 4 to 9 Hz (t test: t6 = 3.52, P = 0.01; Fig. 4A; Materials and Methods). Moreover, we found that the amount of coherence in a given electrode pair was positively correlated with the interregional grid-like representations across this pair [linear mixed effect (LME) model, βvmPFC->EC as a dependent variable; normalized imaginary coherence as a fixed effect and subject IDs as random effects, t112 = 2.85, P = 0.005; t112 = 3.47, P = 0.0007 after controlling for basic differences in oscillatory power between channel pairs; Fig. 4B]. This was not the case for other symmetries (all t112 < 1.38, P > 0.17; fig. S4C).

Fig. 4.

Theta-based interactions between vmPFC and EC carry information about grid-like representations and predict spatial navigation performance.

(A) Top: Example raw trace during the retrieval phase illustrating strong theta band coherence between vmPFC and EC. Bottom: Z-scored imaginary coherence between vmPFC and EC signals. Green bar at the top denotes significant cluster after correction for multiple comparisons (**P < 0.01). Right: Z-scored imaginary coherence of vmPFC and EC at 4 to 9 Hz. Error bars indicate SEM across participants. Gray dots indicate value from individual participants. (B) Electrode-pair level (n = 114) z-scored βvmPFC->EC scores are correlated with z-scored 4 to 9 Hz imaginary coherence (LME, **P < 0.01). (C) Left: Trial-wise z-scored 4 to 9 Hz of imaginary coherence for one example vmPFC-EC pair correlated with drop error. Right: Slope of the best-fitting line of all pairs (LME, **P < 0.01). The negative slope indicates that imaginary coherence is higher in trials with better performance. (D) Mean z-scored 4 to 9 Hz imaginary coherence is significantly different between good and bad trials (LME, ***P < 0.001). Error bars in (C) and (D) indicate SEM across electrode pairs. Gray dots in (C) and (D) indicate value from each pair. (E) Angular histograms of mean vmPFC-EC phase lag (n = 114 pairs) in the theta band (Rayleigh test, ***P < 0.001). (F) Granger causality analyses show significant influence from vmPFC to EC (top) but not in the reverse direction (bottom). The gray lines denote 95% confidence intervals of the null distribution. a.u., arbitrary units.

Next, we examined the relationship between coherence and spatial memory performance (i.e., drop error). For each pair, we calculated the slope between trial-wise z-scored 4 to 9 Hz imaginary coherence and drop error. We modeled the slope as a fixed effect and subject IDs as random effects and found a significant negative correlation between trial-wise theta coherence and drop error (LME test: t112 = −4.17, P < 0.001; Fig. 4C), indicating that stronger vmPFC-EC coherence was associated with better memory retrieval. An exemplary effect in one subject is shown in Fig. 4C. In a complementary analysis, we averaged coherence in three conditions (i.e., good, middle, and bad trials; see Fig. 1B). Again, we found a significant difference between good and bad trials at the electrode-pair level, using LME analyses including trial condition as fixed effect and electrode pairs nested in subjects as random-effect variables (t112 = 3.57, P < 0.01; Fig. 4D). When applying the same analysis to the reencoding phase, we observed the opposite effects between coherence and drop error (fig. S5, A to C; Materials and Methods). This suggests that vmPFC-EC theta synchronization may play distinct roles for memory retrieval and reencoding. In sum, we found that theta coherence was related to the interregional grid-like representations and to the accuracy of spatial memory, indicating that network-level theta grid codes involving both EC and vmPFC support memory-based navigation.

Information transfer between vmPFC and EC

Last, we assessed the direction of functional interactions between EC and vmPFC. We first used phase lag analysis that quantifies the delay at which synchronization between two signals from different regions is maximal (Materials and Methods). Phase differences were calculated for each epoch and for each vmPFC-EC pair using their phase time series and were expressed as the phase locking value (PLV), indicating the consistency of phase differences. PLVs range between 0 and 1, with values approaching 1 if the phase differences between the electrodes vary little across time. This analysis revealed a significant PLV at 4 to 9 Hz (t test: t6 = 3.12, P = 0.02; fig. S5E; Materials and Methods). Phase lags were significantly clustered (Rayleigh’s tests: P < 0.001; Fig. 4E), indicating a consistent delay of around 20 ms of EC relative to vmPFC. We confirmed this finding using spectral Granger causality analysis, which quantifies the strength of directional information flow between two regions in the frequency domain by testing whether the signal from one region (for example, EC) can be better predicted by incorporating information from another region (for example, vmPFC) and vice versa. We found a significant Granger causal influence from the vmPFC to the EC at 4 to 9 Hz during memory retrieval (compared with shuffled data; permutation test, P < 0.05; Fig. 4F), but not in the reverse direction.

DISCUSSION

Grid-like representation in vmPFC

Single-unit recordings in epilepsy patients have revealed first evidence for grid-like cells in the human anterior cingulate cortex during a virtual navigation task (). fMRI navigation studies found that the human vmPFC organizes cognitive maps in a grid-like manner in both spatial and nonspatial space (, , ). Furthermore, a nonhuman primate fMRI study of two-dimensional (2D) decision-making space showed that the medial frontal cortex is involved in inferring novel choices and codes for multidimensional representations of options with a grid-like code (). In line with these studies, we showed grid-like representations in human vmPFC (around area 32), while we did not find grid-like representation in the adjacent OFC that was proposed to be predominantly involved in object-based value representations (). Thus, our iEEG results provide an electrophysiological basis for grid-like representations observed with fMRI. This helps bridge the gap between microscopic and macroscopic signals, as we recently elaborated in a review (). Together, these studies suggest that vmPFC grid representations may play vital roles in cognitive map and goal representations, advancing our understanding of neural computations and memory-based decision-making in human vmPFC (, ).

Grid-like representations across the brain

In this study, we observed grid-like representations in both the vmPFC and the EC. We found that, within subjects, grid orientations appeared to be consistent between the EC and the vmPFC, similar to what was found in previous fMRI studies (, ). A similar grid orientation may be beneficial to coordinate information processing between regions and to generate a coherent cognitive map for the encoding and retrieval of spatial memories (). Anatomical connectivity studies revealed bidirectional projections between these two regions (, ), which may provide a basis for the similarity of grid orientations. Future single-unit recording studies could directly record grid cell activity from these two regions simultaneously and study the phase, orientation, and scale relationship between them, as well as their phase locking to theta oscillations. Recording grid cell activity in different tasks may furthermore help understand whether grid cells in different brain regions support different functions.

Theta synchronization between vmPFC and MTL

We observed pronounced theta oscillations in the human vmPFC in our task. A previous study had shown that single-trial EEG theta power predicted simultaneous functional MRI signals in the anterior cingulate cortex including areas 32 and 24 (), which provides a possible explanation for how our results are linked to grid-like representations observed in prior fMRI studies. Furthermore, by measuring synchronization between vmPFC and EC, we found that theta coherence was related to the accuracy of spatial memory retrieval and to the interregional consistency of grid-like representations. Rodent studies demonstrated phase locking between single-unit activity in the PFC and local field potentials (LFPs) in the MTL as well as synchronous theta oscillations between mPFC and the MTL during active exploration and spatial memory tasks, which was found to correlate with behavioral performance (, ). Human MEG studies additionally found increased theta synchronization between MTL and vmPFC during memory-guided navigation and episodic memory retrieval (, ), as well as memory integration () and scene construction (). Together, these studies indicate network-level interactions between the MTL and the vmPFC that may subserve cognitive functions such as memory and decision-making.

Directional modulation from vmPFC to EC

Electrophysiological investigations in rodents have suggested that the vmPFC drives the hippocampus during location memory retrieval (), proposing that the mPFC influences memory through a top-down control of medial temporal areas to select context-relevant information (). In our task, participants were required to remember eight different objects. When the participants conducted a given trial, they needed not only to extract relevant goal information but also to suppress irrelevant nontarget information. This putatively requires vmPFC control over the retrieval of context-dependent memory representations. Notably, grid cell firing in the EC is also associated with theta oscillatory states (–), which may parallel with grid-like processing in vmPFC. Together, our study found directional modulated theta synchronization between vmPFC and EC, which complemented a previous hypothesis that both coupled theta rhythmicity and spatial cell firing might play important roles in spatial exploration ().

We previously reported that theta oscillations in human EC exhibited grid-like representations during virtual spatial navigation, which we interpreted as a mesoscopic representation of grid cells (, ). Reanalyzing iEEG data during this task, we here reported that theta oscillations in the human vmPFC also exhibited a sixfold rotationally symmetric modulation as a function of movement direction. We found that grid-like representations in the vmPFC exhibit a similar grid orientation with grid-like representations in EC and that theta oscillations in the vmPFC are coherent with EC theta oscillations. Our study thus provides evidence for a network of grid-like representations based on theta oscillations across the brain that supports flexible goal-directed behavior, deepening our understanding of the neural mechanisms underlying human spatial navigation.

MATERIALS AND METHODS

Most of the analysis methods in this study are consistent with those in a previous report (). Here, we introduce them in brief.

Participants

iEEG was recorded from patients who underwent stereotactic electrode implantation for localizing epileptic seizure to guide surgical treatment. This study used 15 subjects (4 female; all 15 right-handed, mean age ± SD: 24 ± 7.6 years; see table S1) based on electrode location.

Navigation task

Participants performed an object-location memory navigation task in a circular virtual arena adapted from (). The environment comprised a grassy plane (diameter of 9500 vu) bounded by a cylindrical cliff. During an initial learning phase, patients were asked to navigate toward and memorize the locations of eight different objects. Afterward, patients completed variable numbers of trials. The locations of the eight objects remained stable throughout each session (but differed between subjects). Each trial consisted of a cue, a retrieval, a feedback, and a reencoding phase (Fig. 1, A and B). During the cue phase, participants were presented with one of eight objects (duration = 2 s). During the retrieval phase, they navigated to the associated object location using the arrow keys on the laptop keyboard. The duration of this phase was self-paced. When they had reached the location that they considered correct, they dropped the object using the space bar. Depending on drop error, patients then received feedback via one of five color-coded smiley faces (duration = 1.5 s). The object then appeared at its correct location, and patients navigated to that location, allowing for relearning. Continuously pressing the forward button automatically accelerated movement speed until maximum speed was reached. Behavioral events including the subject’s location within the virtual environment were written to a log file with a temporal resolution of 10 ms. Triggers were assessed using an independent custom MATLAB (2017b, The MathWorks Inc., MA) script that sent triggers both to the paradigm and to the EEG recording software with randomly jittered intervals between 0.5 and 5 s.

Behavioral analysis

For each time point, we extracted the participant’s location in the virtual environment, which enabled us to derive instantaneous movement direction and speed. We checked whether movement directions were uniformly sampled across 60° space and 360° space using Rayleigh’s tests. This was done to rule out the possibility that our finding of hexadirectional modulation of theta power was due to nonuniform sampling of movement directions. To study the relationship of memory and coherence of vmPFC-EC interactions (Fig. 4), we split trials into three conditions according to the feedback: good trials (dark and light green smileys; drop errors below 1500 vu), medium trials (yellow smileys; drop errors between 1500 and 2500 vu), and bad trials (dark and light red smileys; drop errors above 2500 vu) (Fig. 1B).

iEEG recordings and artifact removal

iEEG recordings were performed at the Yuquan Hospital and the First Affiliated Hospital of People’s Liberation Army (PLA) General Hospital, Beijing, China, as well as at the Department of Epileptology, University of Freiburg, Freiburg, Germany. Our research protocol was approved by the appropriate institutional review boards at each of the three hospital sites. Written informed consent was obtained from all patients. At the recording site in Freiburg, iEEG data were acquired using a Compumedics system (Compumedics, Abbotsford, Victoria, Australia) at a sampling rate of 2000 Hz. At the recording sites in Beijing, iEEG data were acquired using a Nihon Kohden system (Yuquan Hospital) and a Blackrock NeuroPort system (First Affiliated Hospital of PLA General Hospital) at a sampling rate of 2000 Hz as well. Electrodes were provided by HKHS Beijing Health (HKHS Beijing Health Co. Ltd., Beijing, China) at Beijing or Ad-Tech (Ad-Tech, Racine, WI, USA) at Freiburg. Recordings were referenced to Cz (Freiburg) or to one electrode contact located in white matter (Beijing). Regarding the latter, candidate reference electrode contacts located in white matter were chosen by visual inspection of the post-implantation computed tomography (CT) images co-registered onto the pre-implantation MR images (see below). Then, iEEG traces from each candidate reference electrode were visually inspected, and contacts with little or no apparent EEG activity were chosen as the reference for all subsequent recordings. See table S1 for details of electrode contacts from the four different regions of interests.

Recordings of EEG were first inspected for interictal spikes (IISs) and other artifacts by an automatic removal procedure. Epileptic discharges were identified and excluded when one of the following two criteria were met: (i) The envelope of the unfiltered signal was 4 SDs above the baseline, which is the mean value of the entire signal, and (ii) the envelope of the filtered signal (bandpass-filtered in the 25- to 80-Hz range) was 5 SDs above the baseline, following previous studies. All iEEG traces were also visually inspected for epileptic activity, both before and after the automatic removal procedure. The timing of all IISs and other artifacts on each channel were stored for subsequent analysis.

Electrode localization

For patients from Beijing, post-implantation CT images were co-registered onto pre-implantation MR images using FreeSurfer (v6.0.0; surfer.nmr.mgh.harvard.edu/). Intracranial electrodes were identified using clustering-based segmentation and classified according to anatomical landmarks in native space (). For visualization of all subjects’ electrodes on the brain surface, each individual subject’s MRI was normalized to the Montreal Neurological Institute (MNI) standard space. Location of electrode contacts was ascertained by visual inspection of post-implantation CT and/or MRI scans by a consultant neurophysiologist and MRI expert. For patients from Freiburg, for whom one pre-implantation and one post-implantation MRI were available, electrode localization was performed using FSL (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSL) and PyLocator (http://pylocator.thorstenkranz.de/): The post-implantation MR image was co-registered to the pre-implantation MR image. Next, the pre-implantation MR image was skull-stripped and normalized to MNI space. The normalization matrix was then applied to the post-implantation MR image.

Epochs of interest

Analysis of EEG recordings focused on fast movement epochs. As in previous studies, fast movement epochs were defined as the fastest tercile of all movements, separately determined for each participant (table S2). For additional analyses, stationary periods during which patients did not move were also extracted.

Time-frequency analysis

Raw data of the entire experiment were notch-filtered at 50 Hz as well as their harmonics. Next, all epochs that included IISs or other artifacts were excluded from further analyses.

To study the hexadirectional modulation in each predefined frequency band (Fig. 2), the entire data were bandpass-filtered in different frequency bands using eegfilt.m from EEGLAB (low theta, 2 to 5 Hz; high theta, 6 to 9 Hz; beta, 12 to 30 Hz; and low gamma, 30 to 60 Hz). Subsequently, oscillatory power in each band at each electrode channel was extracted using a Hilbert transform. The power values were log-transformed and then z-transformed according to the mean power in that electrode over all fast movement epochs (defined as movements above the cutoff speed; see table S2).

To obtain a power spectrum at a high-frequency resolution, we convolved the raw signal with six-cycle Morlet wavelets at 40 logarithmically spaced frequencies ranging from 2 to 60 Hz. For the analysis of navigation-related oscillations, the power values were log-transformed and then z-transformed according to the mean power in that electrode over all stationary epochs.

Analysis of sixfold modulation of vmPFC

To analyze a potential sixfold modulation of theta power by movement direction, we closely followed the procedure of estimating grid-like representations in previous studies. Recent studies have shown that a sixfold symmetric grid-like code can be extracted from fMRI BOLD (blood oxygen level–dependent) signals, iEEG, and MEG signals. While the cellular mechanisms underlying this effect are still unclear, one possibility is that grid cells show reduced activity during navigation in aligned versus misaligned directions because of more pronounced repetition suppression effects due to narrower spacing of grid cells in these directions; another possibility is that grid-like representations rely on conjunctive grid by head direction cells ().

The initial learning phase was excluded from all analyses. We then averaged power values in each frequency band within consecutive nonoverlapping 10-ms time windows to align the iEEG time courses with the sampling rate of the behavioral data. After extracting periods of fast movement, movement directions (αt) over time were split into six sessions according to the total experiment duration of each participant (irrespective of trial phase). Using one-half of the data (sessions 2, 4, and 6), we identified the “preferred” movement direction ϕ (grid orientation) related to the strongest increases of iEEG power in each frequency band. To do so, we modeled the time series of the power in that band using a general linear model with two regressors, cos(6αt) and sin(6αt). This resulted in two weights of the two regressors, βcos and βsin, which were then used to calculate the grid orientation ϕ = [arctan(βsin/βcos)]/6 using the atan2 function in MATLAB (here termed GLM1). In the second half of the data (sessions 1, 3, and 5), we then used a second general linear model with a single predictor cos(6(αt − ϕ)) to examine whether power increased when the patient moved more aligned with the grid orientation ϕ (here termed GLM2). GLM2 provided a weight β value quantifying the amount of sixfold rotationally symmetric modulation of theta power by movement direction. The higher was the β value, the higher was the band power when the subject was moving closer to the preferred direction ϕ (and the lower when the subject was moving closer to the nonpreferred directions in between). The factor 6 accounted for sixfold rotational symmetry (Fig. 2B).

We also computed a z-scored hexadirectional modulation (Zβ) (Fig. 2C). Each electrode’s Zβ value was computed according to a surrogate distribution of β values that were calculated from shuffled data. Each of the 500 values of the surrogate β values were calculated by circularly shifting the vector of neural data of each fast movement time point relative to the orientation values by a random number of time points. This random circular shift created a random mismatch between the neural and behavioral data while preserving the temporal structures of both time series and offering a conservative control.

To identify which subregions of the vmPFC showed significant hexadirectional modulation, we adapted a cluster-based procedure that is widely used in fMRI studies and has recently been used in the iEEG field as well (, ). All electrodes in a PFC region of interest (ROI) (fig. S1, A and B) were registered to a standard brain atlas in MNI space, allowing for group analyses. For each subject, we selected every vertex on the surface that fell within 15 mm of an electrode’s coordinates (we also considered radii of 9 and 12.5 mm and found similar results) and tagged each of these vertices with the electrode’s Zβ values. Values at each vertex were averaged across electrodes within each patient to create patient-specific Zβ maps. Across patients, we then compared the subject-level Zβ values to zero using a one-sample t test, separately for each vertex. Brain coordinates with less than five subjects were excluded from analyses (, ). Contiguous 3D clusters of significant t values at P < 0.05 were identified. Statistical evaluation of the clusters was performed via cluster-based permutation testing as follows: We created a null distribution of t statistics by randomly sign-flipping each subject’s value and recomputing the group-level t statistics and performed this procedure 1000 times. In each surrogate, the cluster with the highest summed t value was kept. Empirical cluster summed t values were then ranked within the 1000 surrogate summed t values (using four predefined frequency bands and considering a rank of >990 significant, corresponding to P < 0.01). The identified cluster (Fig. 2A) then served as a mask for the selection of PFC in the analyses shown in Fig. 2.

In every given ROI, we used a one-sample t test across subjects (Fig. 2). If a patient was implanted with more than one electrode channel in the vmPFC, β or Zβ values were first averaged across electrode channels. Afterward, the averaged β or Zβ values were fed into second-level statistics across patients. One-sample t tests across patients were performed to test whether hexadirectional modulation of theta power was present at any frequency band.

To illustrate that theta power was modulated by movement direction across the entire 360° range, movement periods during the second half of the data were distributed into 12 different bins (each of 30° size) relative to the patient-specific preferred direction estimated using the first half of the data. Theta power was determined for each of the 12 bins and averaged across patients afterward (Fig. 2D).

Analysis of vmPFC grid orientation

To test whether preferred directions estimated using the first half of the data clustered at a specific angle across patients, we performed a Rayleigh test. If a patient was implanted with more than one electrode channel, the circular mean of the channel-wise preferred directions was fed into the Rayleigh test (Fig. 2G).

To test for grid orientations at the within-subject level, we tested the hypothesis that each subject has a consistent grid orientation across their own electrodes. (i) For every channel from each patient, we calculated the putative grid orientation using one-half of the data; (ii) we then compared this orientation to the orientation in a different channel using the other half of the data. This comparison was done by calculating the absolute angular difference between the preferred orientations. Absolute angular differences ranged from 0° to 30°, because preferred orientations comprised circular values between 0° and 60°; (iii) the results were averaged across channel pairs for each patient; (iv) across patients, we performed one-sample t test to examine whether the mean absolute angular differences were below 15°, which represents chance level (fig. S2A).

Control analyses of sixfold modulation

We performed control analyses to validate the presence of hexadirectional modulation of theta power with regard to sixfold rotational symmetry and to the vmPFC. We checked whether theta power during fast speed movements was modulated by movement direction with regard to other types of rotational symmetry (i.e., four-, five-, seven-, and eightfold rotational symmetry). Again, the same overall analysis procedure as outlined above was performed, now using the factors 4, 5, 7, or 8, respectively, when setting up the regressors for the first and second GLM2 [e.g., for eightfold rotational symmetry, cos(8αt) and sin(8αt) for GLM1, and cos(8(αt − ϕ)) for GLM2]. Second-level t tests were performed across patients afterward. As a region-level control analysis, we performed the analysis of hexadirectional modulation of iEEG power by movement direction for electrode channels located in dmPFC and OFC (Fig. 2, E and F).

Relationship between grid orientations in EC and vmPFC

To analyze the relationship of the grid-like representations in these two brain regions, we focused on seven subjects who had both contacts in EC and vmPFC. For each subject, we selected all electrodes separately for each region and then paired them (for example, if a subject had two EC electrodes and three vmPFC electrodes, this results in 2 × 3 = 6 pairs). This yielded 114 EC-vmPFC pairs in total.

We divided the data into six sessions, where three sessions were drawn for training and the other three sessions were drawn for testing. For each vmPFC-EC electrode pair, (i) we computed the grid orientation in the EC using the first half of the data and tested this grid orientation on the second half of the data from the vmPFC to obtain the regression coefficient β1; (ii) we computed the grid orientation in the EC using the second half of the data and tested this grid orientation on the first half of the data from the vmPFC to obtain β2; (iii) we then averaged β1 and β2 to obtain the grid-like modulation index βEC->vmPFC; (iv) βEC->vmPFC values were averaged across electrode pairs from the same subject. The averaged βEC->vmPFC values were then fed into group-level statistics across the seven patients. To assess the statistical significance of each of the grid-like modulation indices (βEC->vmPFC, averaged across subjects), we performed a permutation test by creating a null distribution of βnull EC->vmPFC and tested whether the mean value βEC->vmPFC across patients was significantly higher than the 99th percentile of the null distribution. In detail, we obtained 500 surrogate group mean β values across patients. Specifically, for each electrode pair of each patient, (i) we circularly shifted the vector of movement direction data relative to the neural data by a random number of time point, and (ii) we recomputed βEC->vmPFC according to the procedure described above to get surrogate values βnull-EC->vmPFC. This random circular shift created a random mismatch between the neural activity and the behavioral data, while preserving the temporal structure of both time series. Note that, in each loop, we also keep the random time points the same for each electrode pair from the same patient to preserve the temporal structures of electrode pairs. This also ensured that the consistency of grid-like orientations cannot only be explained by power correlations between the two regions. To further corroborate our results, we also performed the analysis in the reverse direction, where we used vmPFC data to estimate the preferred grid orientations and the EC data to test the consistency of the preferred orientations across both regions, resulting in βvmPFC->EC values. Again, we created null distributions of surrogate values (βnull vmPFC->EC). We z-scored the raw βvmPFC->EC values to obtain ZβvmPFC->EC to create Fig. 3A and, similarly, converted the raw βEC->vmPFC values to obtain ZβEC->vmPFC to create fig. S4B. As a control analysis, we used a similar procedure to obtain the Zβ values for other folds (Fig. 3A and fig. S4C).

Next, to examine the specific pattern of vmPFC->EC sixfold modulation (Fig. 3B), we used the previous method () to separate all aligned and misaligned trajectories and construct 12 regressors. In detail, we aligned the trajectory orientations to the grid orientation and split these aligned trajectories into 12 equal bins of 30°. We thus created 12 separate regressors for trajectories that belonged to one of these 12 bins. The 12 resulting regression coefficients were then averaged and expressed in Fig. 3B.

Analysis of hexadirectional modulation across the whole brain

To illustrate the distribution of hexadirectional modulation across the whole brain, we analyzed all channels in our dataset (2745 in total) and tested for statistical significance at the electrode level. To assess the statistical significance of the hexadirectional modulation from individual electrodes, we measured the rank of each electrode’s raw β value relative to the null distribution of surrogate β values. We designated an electrode as showing a significant hexadirectional modulation if its rank exceeded the 95th percentile of the null distribution. To assess the significance of the electrode counts for each brain vertex on the surface, we (i) counted the number of significant electrodes and total numbers of electrodes that fell within 15 mm of this vertex’s coordinates, and (ii) we used a binomial test to compare the observed value to the 5% expected false-positive rate at a chance rate of 0.05. In line with the results shown, we found that the percentage of electrodes showing grid-like representations was high in vmPFC (surface space; fig. S7A) and EC (volume space; fig. S7B). In addition, we found posterior cingulate cortex/retrosplenial cortex (PCC/RSC) and superior temporal sulcus (STS) (fig. S7A) showing a significant hexadirectional modulation. However, because of the dependence of LFP signals in neighboring electrodes for each subject, this analysis may overestimate the percentage of significant electrodes. Thus, this result shall mainly illustrate which other brain regions also show hexadirectional modulation of theta power by movement direction in our dataset.

Analysis of EC-vmPFC connectivity

To measure functional connectivity between EC and vmPFC, we extracted continuous fast-speed movement segments with a duration of more than 1 s and treated them as different epochs. For each pair of electrodes, the imaginary coherence of each epoch was calculated. Imaginary coherence () was applied to estimate spectral coherence without spurious connectivity because of volume conduction. The selection of the white matter reference electrode in vmPFC and EC electrodes was different, which further reduced the introduction of common signals. Coherence between two time series x(t) and y(t) was then estimated as a function of frequencywhere S(f) and S(f) are the power spectral densities of electrodes x and y, respectively, and S(f) is the cross-spectral density between two electrodes. Imaginary coherence Icoh = Imag(C(f)) is the imaginary part of coherence. The average imaginary coherence values across epochs were used as indicators of the strength of connectivity of the electrode pair. The statistical significance was then estimated for each electrode pair using a permutation test, in which a null distribution was created by shuffling epoch labels. EC time series were randomly exchanged among epochs and then the coherence was recomputed, and this procedure was repeated 500 times. The observed data were then z-transformed according to this null distribution to estimate a Z-Icoh value for each frequency. The Z-Icoh values within the 4- to 9-Hz range were averaged across electrode pairs of each subject to obtain subject-wise Z-Icoh4-9Hz values. Afterward, the averaged Z-Icoh4-9Hz values were entered into a second-level statistical test across patients (Fig. 4A). To study the relationship between functional connectivity and grid-like representations, we calculated correlations between Z-Icoh4-9Hz and ZβvmPFC->EC across channel pairs using a linear mixed model to disentangle within-subject and between-subject effects (Fig. 4B).

To study the relationship between functional connectivity and spatial memory performance, we examined the correlation between coherence and drop error. To do this, we first averaged Z-Icoh4-9Hz values of fast-speed movement epochs in the retrieval phase of each trial. Second, we extracted trial-by-trial values of Z-Icoh4-9Hz for each vmPFC-EC electrode pair. Third, we fit the data using a linear regression to obtain slopes of each pair (Fig. 4C). In a complementary analysis, we split trials into three groups (i.e., good, medium, and bad performance) depending on the feedback. We found a significant effect between good and bad performance trials across pairs of electrodes using a linear mixed model to disentangle within-subject and between-subject effects (Fig. 4D).

Analysis of directional information transfer between EC and vmPFC

To investigate the direction of information flow between the two brain regions, we used two methods: an indirect method using phase lag analysis and a direct method using spectral Granger causality.

Phase lag analysis

The PLV quantifies the consistency of phase differences between a pair of electrodes. Phase differences were calculated for each epoch and electrode pair (x, y) using their phase time series. The phase difference between the two time courses ∅[n] indicates the synchronization between two signals and is expressed as the PLVwhere N denotes total time points. PLVs range between 0 and 1, with values approaching 1 if the phase differences between the electrodes vary little across time. Then, the statistical significance for each electrode pair was estimated using a permutation test to get ZPLV4-9Hz (fig. S5E).

For each pair of electrodes, we extracted the time series of complex numbers at 7 Hz (center overlapping frequency between vmPFC and EC) by convolving both raw time series with a complex Morlet wavelet (wavelet number = 6) centered at 7 Hz. We measured phase lags between EC and vmPFC by calculating the mean phase lag between the phase time series of each pair of electrodes

The average phase lag values of different epochs were used as indicators of the lag between them (Fig. 4E). Notably, for oscillatory signals, these lags do not provide information about the direction of information flow.

Spectral Granger causality

Spectral Granger causality quantifies the strength of directional influences between iEEG signals in the frequency domain by testing whether iEEG signals from one structure can be better predicted by incorporating information from the respective other structure than by only taking preceding information from the same structure into account. For each pair, the continuous fast-speed movement epochs of iEEG raw data were first down-sampled to 250 Hz before fitting to an autoregressive model and computing spectral Granger causality. These data were detrended and normalized over time and across epochs to increase stationarity. Model order was determined by the Bayesian information criterion using the Multivariate Granger Causality Toolbox (), providing a compromise between spectral resolution and overparameterization. Model orders were estimated for each patient and converged to numbers varying from 10 to 16, corresponding to maximal lags between 40 and 64 ms. We then averaged all pairs from one subject to obtain the spectral Granger causality for this subject (Fig. 4F). Next, we tested whether this directionality was regionally specific. For each pair, a null distribution of spectral Granger causality values was created by shuffling the epoch labels 500 times before calculating the Granger prediction values. The threshold was determined by the 99th percentile value of this distribution under the null hypothesis. We averaged the threshold of all pairs from one subject and then averaged across subject to get spectral Granger causality at the shuffle level (gray lines in Fig. 4F).

Statistical analysis

All analyses were performed using custom MATLAB scripts or R software. Unless specifically mentioned, statistical tests were conducted across subjects. Otherwise, they were performed at the electrode (or electrode pair) level using a linear mixed model, rather than a simple t test. Given the fact that different participants contributed different numbers of electrodes, this approach maximizes statistical sensitivity while separating within-subject and between-subject variance [see (–) for a similar approach]. LME analyses were carried out using the LME4 package implemented in R. Power or coherence values were fitted with a random intercept model formulated as follows: Y ~ X1 + X2 + (1|Patient), where X1 and X2 are fixed effects, and the terms in parentheses are the random effect of “Patient.”