Introduction

The cognitive map, initially introduced by Tolman (Tolman, 1948) and later supported by the discovery of place cells in the hippocampus (HPC) (O’Keefe and Dostrovsky, 1971), has served as a framework for spatial navigation. Grid cells in the upstream entorhinal cortex (EC), known for their hexagonal firing pattern to provide the metric for space (Hafting et al., 2005), have also been shown to represent conceptual space beyond the physical reference frame (Constantinescu et al., 2016; Bao et al., 2019; Raithel et al., 2023). These findings suggest that the EC-HPC circuit fundamentally organizes spatial and non-spatial knowledge (Epstein et al., 2017; Behrens et al., 2018; Bottini and Doeller, 2020; Park et al., 2020), and to guide flexible behaviors such as retrieving knowledge from the past and making decisions for the future (Addis et al., 2007; Hassabis and Maguire, 2007; Schacter et al., 2012). However, a key question remains unresolved, how the EC and HPC collaboratively construct a conceptual map is still largely unknown.

To flexibly navigate, it is necessary to maintain internal representations of current location, goal location, and locations along the planned trajectory (Nyberg et al., 2022). This process must involve the activation of a population of spatial cells. Ormond and O’Keefe (2022) identified convergence sinks of place cells (ConSinks). These ConSinks, characterized by the vector fields of single neurons, showed dynamic activity distributions of place cells across space as a function of goal locations. These vector fields provided activity gradients that guided movement from current locations, thereby highlighting the vector-based representation in the HPC that supports flexible navigation. Muhle-Karbe et al. (2023) investigated distortion in the geometric representation of cognitive maps modulated by context-based goals. They reported the “goal-based spatial compression”, demonstrating that goal-dependent cognitive maps were encoded together in the HPC as if participants were simultaneously imagining themselves at those locations. These findings, as consistently observed in the brains of bats (Sarel et al., 2017), suggest that the HPC integrates discrete information, including both current and perspective locations along a spatial trajectory, into a vectorial representation. They also raise the question of what role grid cells in the EC play in this process.

Anatomically, the entorhinal cortex (EC) serves as one of the major sources of input to the downstream HPC (Witter and Amaral, 1991; van Groen et al., 2003; Garcia and Buffalo, 2020). Functionally, grid cells in the medial EC exhibited multiplexed and heterogeneous responses corresponding to multiple navigation-relevant variables, such as position, direction, and speed (Hardcastle et al., 2017). These dynamically mixed codes in the grid system may reflect updates of self-motion-based information before it is anchored into allocentric maps in the HPC (Sargolini et al., 2006). Computationally, grid cells have been proposed to provide the foundation for hippocampal place field formation by integrating multi-scale grid codes across different grid modules (Solstad et al., 2006; De Almeida et al., 2009; Bush et al., 2014; Bush et al., 2015; Bicanski and Burgess, 2019). Physiologically, EC lesions have been shown to disrupt the precision and stability of place fields (Hales et al., 2014), leading to reduced discharge rates and field sizes (Van Cauter et al., 2008). Collectively, these findings support the idea that the EC serves as a geometric framework for discrete spatial information, which is integrated in the HPC.

The periodic coding of grid cell population is known to exhibit nearly invariant orientation and spacing (Hafting et al., 2005; Sargolini et al., 2006; Krupic et al., 2012; Gardner et al., 2022). We hypothesize that the three primary axes of the hexagonal grid cell codes may manifest as vector fields observed in the HPC. If this conjecture holds true, the HPC should exhibit a periodic activity pattern as a function of spatial directions, phase-locked with that of the EC. However, simultaneously population-level recording EC and HPC activity poses significant challenges for electrophysiology, two-photon calcium imaging, electroencephalography, and Magnetoencephalography. Therefore, we employed functional MRI, which has successfully identified 6-fold periodicity in the EC during navigation in both spatial and conceptual spaces (Doeller et al., 2010; Constantinescu et al., 2016; Bao et al., 2019; Wagner et al., 2023; Raithel et al., 2023), to investigate the structured neural periodicity in the HPC.

A novel 3D object, named Greeble (Fig. 1a) (Gauthier and Tarr, 1997), was used to create a conceptual Greeble space. Within this space, locations were represented by Greeble variants characterized by two features, namely “Loogit” and “Vacso”. The feature length defined the two dimensions of the space. The Greeble positioned at the center of the space served as the prototype and acted as the target. Participants were instructed to morph Greeble variants to match this target prototype (Fig. 1b). This process generated a sequence of Greebles that behaved like movements along a navigational trajectory in a two-dimensional space (Fig. 1c). To ensure a comprehensive exploration of the Greeble space for detecting HPC periodicity, Greeble variants were pseudo- randomly sampled at the periphery of space. This provided a high-resolution vector of conceptual directions ranging from 0° to 360° (i.e., the orange locations in Fig. 1c), while controlling for distance. As the result, we observed a 3-fold periodicity in HPC activity within the direction domain, which was cross-validated using sinusoidal and spectral analyses, despite the participant’s unawareness of the existence of the Greeble space. The spatial phase of the HPC was coupled with a 6-fold periodicity in the EC, as well as a 3- fold periodicity in participants’ behavioral performance. Additionally, the EC-HPC PhaseSync model, developed to simulate EC projections into the HPC, validated the emergence of 3-fold activity periodicity under randomized goal locations. These empirical findings support our hypothesis that the formation of hippocampal vector fields relies on EC projections.

Fig. 1.

Results

Behavioral results

Participants adjusted Greeble features using two response boxes inside the MRI scanner. Each response box controlled one feature, with two buttons used to stretch and shrink the feature, respectively. To prevent the use of an “horizontal-vertical movement” strategy, where movements follow the cardinal directions (e.g., repeatedly adjusting one feature towards Northward and then Eastward), participants were encouraged to adjust both features simultaneously, promoting variance in movement directions for detecting neural periodicity. This desired strategy was referred to as “Radial Adjustment” involved making a one-unit adjustment to one feature with the left hand followed by a one-unit adjustment to the other feature with the right hand.

To examine participants’ strategy, we calculated the number of directions across every three consecutive steps within their original movement trajectories (Fig. S1). Participants showed an average of 3.8 directions per trajectory, which was significantly higher than the baseline of two directions predicted by the “horizontal-vertical movement” strategy (t(32) = 15.76, p < 0.001, two-tailed; Cohen’s d = 2.78), confirming stable movements across participants along intercardinal directions in the Greeble space.

To eliminate the impact of learning effect on potential BOLD signal periodicity, all participants completed a training task one day prior to the MRI scanning to familiarize themselves with Greebles morphing (see methods for details). Behavioral performance was defined as the sum of trajectory length and error size, with superior performance characterized by smaller scores. Learning effects were observed during the practice experiment on day 1 (Fig. S2; t(32) = −2.46, p = 0.019, two-tailed; Cohen’s d = 0.44), but were no longer present during the MRI experiment on day 2 (t(32) = −0.74, p = 0.462, two-tailed; Cohen’s d = 0.13). These results reflected a stable task accuracy across participants (Fig. 1d; trajectories of individual participants are shown in Fig. S3).

6-fold periodicity in the EC

Using sinusoidal analysis, the grid orientation, reflecting the allocentric direction of the grid axes, was calculated for each voxel within the bilateral EC using half of the dataset (Experimental session 1, 3, 5, 7). Participants’ movement directions were determined by the 2D position of Greeble variants relative to the ending locations, with 0° arbitrarily defined as movement from the East (Fig. 2a). Significant deviations from uniformly distributed grid orientations were observed in 30 out of 33 participants (Fig. S4, p < 0.05; Rayleigh test of uniformity, Bonferroni corrected), confirming that grid orientations across voxels in the EC tend to be consistent in the majority of participants.

Fig. 2.

Hexagonal 6-fold activity was examined using the other half of the dataset (Experimental session 2, 4, 6, 8), with movement directions calibrated according to participant- dependent grid orientations. This analysis identified a significant cluster in the right EC (Fig. 1b; initial threshold: p = 0.05, two-tailed; cluster-based SVC correction for multiple comparisons: p < 0.05; Cohen’s d = 0.63; Peak MNI coordinate: 32, −6, −30). Four control analysis were performed using rotationally symmetric parameters (3-, 4-, 5-, and 7-fold) to validate the robustness of this cluster in representing the 6-fold periodicity. These analyses only revealed significant periodic modulation at 6-fold (Fig. 2c; t(32) = 3.56, p = 0.006, two-tailed, corrected for multiple comparisons; Cohen’s d = 0.62). In contrast, no significant results were detected for 3-fold periodicity (120° periodicity; t(32) = 1.10, p = 0.28, Cohen’s d = 0.19), 4-fold periodicity (90° periodicity; t(32) = 0.31, p = 0.76, Cohen’s d = 0.05), 5-fold periodicity (72° periodicity; t(32) = −0.21, p = 0.83, Cohen’s d = −0.04), or 7-fold periodicity (51.4° periodicity; t(32) = 0.88, p = 0.39, Cohen’s d = 0.16).

We further inspected the hexagonal activity pattern in 1D directional domain reported in previous papers (Doeller et al., 2010; Constantinescu et al., 2016; Bao et al., 2019; Wagner et al., 2023; Raithel et al., 2023), where EC activity changes as a function of movement directions. The Bold signals reconstructed from sinusoidal analysis were averaged across voxels within a hand-drawn bilateral EC mask. A stable 6-fold periodicity was observed (Fig. 2d, left; p < 0.05, permutation corrected for multiple comparisons), characterized by stronger activity when movement directions were aligned with the calibrated grid axes and weaker activity when they were misaligned, confirming a 60° periodicity of EC activity in the Greeble space (Fig. 2d, right).

3-fold periodicity in the HPC

Potential neural periodicity in the HPC was examined using spectral analysis. This method does not require calibrations of spatial phase, as is necessary for the sinusoidal analysis of 6-fold periodicity. Therefore, the BOLD signals from all eight experimental sessions were modelled using a GLM to estimate direction-dependent activity. Specifically, participants’ movement directions were down-sampled into 10° bins, resulting in a total of 36 directional bins (e.g., 0°, 10°, and 20°, etc.). This approach enables the detection of activity periodicity ranging from 0- to 18-fold (Fig. 3a). These directional bins were entered into the GLM as binary regressors. Fast Discrete Fourier Transform (FFT) was performed on the direction-dependent parametric maps sorted in ascending order of spatial directions. Spectral magnitude maps were then extracted for each of the 0- to 18-fold periodicities (see Methods for details).

Fig. 3.

A significant cluster of 3-fold periodicity was observed in the bilateral HPC (Fig. 3b; initial threshold: P = 0.05, two-tailed; Cluster-based SVC corrected for multiple comparisons: p < 0.05; Cohens’ d = 1.06; Peak MNI coordinate: −24, −20, −18). Moreover, a significant cluster of 6-fold periodicity was identified in the right EC (Fig. 3c; initial threshold: p = 0.05, two-tailed; Cluster-based SVC-corrected for multiple comparisons: p < 0.05; Cohen’s d = 1.27; Peak MNI coordinate: −22, −14, −30), thereby replicating the finding from the sinusoidal analysis.

The hippocampal 3-fold periodicity has not been previously reported in the spatial domain by either neurophysiological or fMRI studies. To validate its reliability, we conducted additional analyses. First, sinusoidal analysis confirmed significant 3-fold clusters in the bilateral HPC (Fig. 3d; initial threshold: p = 0.05, two-tailed; Cluster-based SVC-corrected for multiple comparisons: p < 0.05; Cohen’s d = 0.68; Peak MNI coordinate: −24, −18, −12). ROI analysis indicated that HPC activity selectively fluctuated at 3-fold periodicity (t(32) = 3.94, p= 0.002, corrected for multiple comparisons; Cohen’s d = 0.70), whereas no significant activity periodicity was observed at other spatial fold (4-fold periodicity: t(32) = 0.09, p = 0.93, Cohen’s d = 0.02; 5-fold periodicity: t(32) = 0.31, p = 0.76, Cohen’s d = 0.05; 6-fold periodicity: t(32) = 2.36, p = 0.12, Cohen’s d = 0.42; 7-fold periodicity t(32) = 1.21, p = 0.24, Cohen’s d = 0.21). Reconstructed BOLD signals showed a significant 3-fold periodic pattern across voxels within an anatomical mask of the bilateral HPC (Fig. 3f, left; p < 0.05, corrected for multiple comparisons). In sum, both the spectral and sinusoidal analyses confirmed a 3-fold periodic activity in the HPC across the conceptual directions within the 2D Greeble space (Fig. 3f, right), which was distinct from the 6-fold periodicity observed in the EC.

Second, the 3-fold HPC periodicity was not affected by the precision of directional sampling, as consistent findings were observed when 20° bins were applied (Fig. S5; initial threshold: p = 0.05, two-tailed; Cluster-based SVC corrected for multiple comparisons: p < 0.05; Cohen’s d = 1.18; Peak MNI coordinate: 22, −24, −14). Third, we independently analyzed HPC and EC periodicity using spectral analysis for each of the three site-dependent experimental groups (See methods for details). These analyses showed reliable activity periodicities in both brain areas across experimental groups (Fig. S6; initial threshold: p = 0.05, two-tailed. Cluster-based SVC corrected for multiple comparisons: p < 0.05). Forth, we further examined the whole brain representation of both 3- and 6-fold activity periodicity, respectively (Fig. S7). The 3-fold periodicity revealed significant involvement of the medial prefrontal cortex (mPFC), precuneus (PCu), and parietal cortex (PC)(initial threshold: p = 0.05, two-tailed. Whole brain correction for multiple comparisons: p < 0.05), suggesting the engagement of the default mode network. In contrast, the 6-fold periodicity highlighted the Salience network, including the anterior cingulate cortex (ACC) and insular cortex (INS), in addition to the EC (initial threshold: p = 0.05, two-tailed; Whole brain cluster-based correction for multiple comparisons: p < 0.05). These results demonstrated distinct functional networks involved in integrating spatial metric (6-fold periodicity) and in judging self-motion within conceptual space (3-fold periodicity).

To investigate the relationship between the observed activity periodicity in the HPC and EC, we first questioned whether phase-coupling could be found between the two areas? The spatial orientations of both areas’ activity, estimated through rotationally symmetric parameters in sinusoidal analysis, reflected the spatial phase of brain activity. We therefore examined the cross-participant correlations between the orientations of both areas using Pearson correlation. The result suggested significant phase coupling between the EC and HPC (Fig. 4a; r = 0.41, p = 0.019).

Fig. 4.

We further examined amplitude-phase coupling using a cross-frequency amplitude-phase modulation analysis (Canolty et al., 2006). This method measured the composite signal z for each sample of EC and HPC activities, where z = A_ERC e^iφHPC, with A_ERC representing the amplitude envelope of the 6-fold activity extracted from the the EC, and φ_HPC representing the phase of the 3-fold activity from the HPC. The activities of both areas were reconstructed using a sinusoidal model based on participants’ original movement direction, instead of the orientation-calibrated directions. The modulation index M, defined as the absolute value of the mean of z values, was calculated for each phase of the HPC activity in the direction domain. This analysis confirmed significantly higher coupling strength between the amplitude of the EC and the phase of the HPC compared to a surrogate dataset created by spatially offsetting the EC amplitude relative to the HPC phase by spatial lags (t(32) = 7.24, p < 0.001; Cohen’s d = 1.14). This analysis further confirmed that the EC amplitude phase-locked to the HPC phase.

Next, given our hypothesis that the vectorial representations of the HPC depends on projections of grid cell activity in the EC, particularly along the three primary axes of hexagonal codes, we asked whether the periodic activities of the EC and HPC are peak-overlapped? To test this, the phases of the HPC were classified into four bins . The modulation index M was then calculated for each bin. This analysis revealed the peak offset of brain activity by identifying at which HPC phase the EC showed the strongest amplitude. If the activities of both areas were peak- overlapped (i.e., no peak offset), a symmetric pattern of M centered at phase 0 of the HPC would be expected, given the cosine function used during activity reconstruction (as shown by the one-cycle-based schema in Fig. 4c). Indeed, phase and phase [0 to ] showed significantly higher coupling strength than phase (Fig. 4d; F(1, 130) = 218.99, p < 0.001, Cohen’s d = 1.36). Moreover, no significant difference in M was found between phase and phase (t(32) = 0.05, p = 0.96; Cohen’s d = 0.01), nor between phase (t(32) = 1.11, p = 0.27; Cohen’s d = 0.19). This inverted U-shaped pattern of the modulation index M, characterized by higher M values near phase 0 and lower M values near phase −π and π, suggested that no peak offset was found between EC and HPC activity. Notably, these results differ from previous studies that examined neural activity coherence in the temporal domain (Buzsaki, 2006; Maris et al., 2011; Friese et al., 2013). Instead, our findings indicate selective alignments of activity peaks between the 3-fold HPC periodicity and the 6-fold EC periodicity in the conceptual direction domain.

3-fold periodicity on behavior performance

Considering the circular connections of the medial temporal lobe (MTL), where the EC and HPC reside, with the visuomotor system, HPC-based conceptual representations might be integrated with visually guided motor planning, particularly through parietal and retrosplenial pathways. It is plausible that the activity periodicity of the HPC and EC propagates to egocentric cortical areas, such as the parietal cortex (PC), thereby influencing visuospatial perception performance (Fig. 5a). To test this hypothesis, we examined participants’ behavioral performance using two metrics introduced above: (1) the trajectory length of each trajectory, and (2) the error size between the actual ending location and the goal location, given that participants seldom stopped exactly at the goal. This procedure generated a vector of task performance in the directional domain for each participant (Fig. 5b), with superior behavioral performance indexed by a composite score of shorter travel lengths and smaller error size.

Fig. 5

Using spectral analysis, we identified a significant periodicity in participants’ behavioral performance, with spectral power peaking at a 3-fold symmetry (Fig. 5c; p < 0.05, permutation corrected for multiple comparisons), demonstrating that participants’ visuospatial perception performance fluctuated as a function of movement directions. In contrast, no significant spectral power was observed at other spatial folds, including 6- fold behavioral periodicity (p > 0.05). To further investigate the causal relationship in brain-behavior coupling, we calculated the phase-lag index between the 3-fold behavioral periodicity and HPC activity for each participant (Stam et al., 2007). This analysis revealed a significantly higher phase-locking value (PLV) compared to surrogate dataset (Fig. 5d; t(32) = 8.10, p < 0.001; Cohen’s d = 1.14), confirming the phase coupling between the 3-fold HPC activity and the 3-fold behavioral performance. Contrary to the commonly held assumption of isotropic task performance across movement directions, these results demonstrated a biased behavioral pattern influenced by HPC and EC activity periodicity, with superior navigational performance (i.e., shorter trajectory length and smaller error size) observed when movement directions aligned with grid axes compared to when they were misaligned.

The “EC-HPC PhaseSync” Model

The above empirical results elucidated that the vectorial representations of the HPC, characterized by 3-fold activity periodicity, was spatially coherent with EC population activity over the direction domain. To further investigate this phenomenon, we developed a cognitive model, termed the “EC-HPC PhaseSync” model (Fig. 6), which demonstrates the possible computational mechanism of the EC-HPC circuit, through simulating how the periodic representations emerged in the HPC influences spatial movements toward the center of Greeble space, as well as qualifying the reliability of the 3-fold periodic performance, particularly when goal locations are randomly distributed across the conceptual space.

Fig. 6.

The EC grid code G was simulated using a cosine grating model (O’keefe and Burgess, 2005; Burgess et al., 2007; Blair et al., 2007; Bush and Schmidt-Hieber, 2018; see methods for details). A population of 45-by-45 grid cells were totally generated to tessellate the space (Fig. 6a, left), with spatial orientations kept constant while spatial phases varied across locations. The population activity of the EC, elicited by movements from a start location, was termed the trajectory code V, simulated by summing the grid codes across locations along the trajectory (Fig. 6a, right). Notably, the trajectory code V is exactly identical for movement direction φ and φ + 180°, naturally reflect the spatial orientation Ψ ∈ [0, π) across the space.

The downstream HPC activity was simulated using a two-dimensional δ vector to represent spatial orientations. Each δ value, encoded by a single HPC neuron, was obtained by linearly summing EC activity across locations on the trajectory code V (Fig. 6b, black dots). The spatial periodicity and spatial phase of δ depend on the number of primary axes of grid codes and grid orientation, respectively, while the magnitude of δ indicates the degree of alignment between the orientation Ψ and the grid axes. This model hypothesizes that spatial directions are embedded within the 3-fold structure of δ originally projected from the EC grid cells, enabling allocentric directions to be distinguished by the amplitude and phase of δ vector. The vectorial representation C of the HPC was simulated by integrating δ vector with a Gaussian-based distance vector centered at the goal location, thereby driving goal-directed movements (Fig. 6c).

During simulations, 100 goal locations were randomly generated within the space, each paired with 120 starting locations, arranged in a ring-shaped pattern near the boundary. At each simulation step, a vector of HPC activity representing the locations surrounding the current location was extracted from C, and movement was driven toward the next location with the strongest activity, following winner-take-all dynamics. The step size, initialized as 1 pixel, was dynamically adjusted using a stepwise approach, with a maximum limit of 5 pixels; for example, the step size was increased by 1 pixel if no stronger activity was found among surrounding locations. Model performance was defined by the trajectory length, calculated as the sum of Euclidean distances between consecutive locations, with superior performance corresponding to shorter trajectory length. The model performance was sorted by movement directions in ascending order prior to spectral analysis. Ultimately, the simulation confirmed a significant 3-fold periodicity in trajectory length across 100 navigations tasks with random goal locations (Fig. 6c; Right).

Discussion

Spatial navigation, the cognitive process of retrieving spatial relationships from one moment to the next in both conceptual (Epstein et al., 2017; Behrens et al., 2018; Bottini and Doeller, 2020; Park et al., 2020) and physical spaces (Tolman, 1948; Sargolini et al., 2006; Gardner et al., 2022; Ormond and O’Keefe, 2022), relies on the EC-HPC circuit. These two areas are known to have distinct functions. The EC encodes continuous space with a metric representation (Hafting et al., 2005), while the HPC represents localized spatial locations (O’Keefe and Dostrovsky, 1971). A longstanding puzzle is how these brain areas collaborate to drive navigation. In the present study, we modulated the EC and HPC activity using conceptual directions, and identified 6-fold periodicity from EC activity and 3-fold periodicity from HPC activity. The periodic activities of both areas were synchronized and peaked-overlapped. Additionally, 3-fold periodicity was also observed in participants’ behavioral performance, which was synchronized with the HPC activity, suggesting a direct influence of periodic brain activity on visuospatial perception. These empirical findings provide a potential explanation to the formation of vectorial representation observed in the HPC (Sarel et al., 2017; Ormond and O’Keefe, 2022; Muhle-Karbe et al.,2023).

The default mode network and the salience network, revealed by spectral analysis showing 3- and 6-fold periodicity, reflect a widespread distribution of conceptual- dependent periodic neural activity across the brain beyond the EC-HPC circuit. Although the brain areas within these networks were engaged through the neural periodicity of the EC and HPC, this does not imply that the origin is the hexagonal representations of grid cells. Instead, these periodic activities are more likely relying on the hierarchical organization of upstream cortical areas, such as the ganglion cells first reported by Hubel and Wiesel (Hubel and Wiesel, 1959; Hubel and Wiesel, 1962), which selectively respond to visual locations and provide external visual cues in order to support egocentric navigation, especially when an allocentric map is not accessible. Moreover, rhythmic spatial attention further supports the widespread presence of periodic representation in the brain. Fiebelkorn et al. (2018) reported that the PC activity fluctuated at 3-8 Hz. This periodic activity predicted the behavioral performance of spatial attention sampling in a discrimination task. Together with the 3-fold behavioral periodicity observed in the present study, these findings suggest that periodic representations might be a generalized mechanism by which the brain structures information beyond time and space.

Supplementing previous navigational models (Bicanski and Burgess, 2019; Gao et al., 2018; Whittington et al., 2020; Edvardsen et al., 2020), the EC-HPC PhaseSync model demonstrates how the EC and HPC are functionally connected in periodic manner. The model was inspired by two factors. First, the EC is known to serve as one of the major projection sources of the HPC (Witter and Amaral, 1991; van Groen et al., 2003; Garcia and Buffalo, 2020). Together with nearly invariant spatial orientation and spacing observed in grid cell population (Hafting et al., 2005; Gardner et al., 2022). The EC provides a foundation for the formation of periodic representations in the HPC. Second, mental planning, characterized by “forward replay” (Dragoi and Tonegawa, 2011; Pfeiffer, 2020), has the capacity to activate grid cell populations that represent sequential experiences even without actual physical movement (Nyberg et al., 2022). The hypothesized trajectory code in our model, which integrates the spatial representations of grid populations from sequential experiences, aligns well with several previous studies, including the oscillatory interference theory (Solstad et al., 2006; De Almeida et al., 2009; Welday et al., 2011; Bush et al., 2014; Bush et al., 2015), the activity pattern of band-like cells observed in the dorsal medial EC (Krupic et al., 2012), and the irregular activity fields of trace cells observed in the HPC (Poulter et al., 2021). In addition, the orientation representation in the HPC hypothesized by the model may correspond to the reorientation problems observed from both rodents and young children (Hermer and Spelke, 1994; Julian et al., 2015; Gallistel, 2017; Julian et al., 2018), where subjects tend to search for targets equally at two direction-opposite locations, potentially reflecting a successful identification of spatial orientation alongside a failure to form an allocentric direction representation based on geometrical cues.

In sum, the present study reported a hippocampal representation exhibiting 3-fold conceptual periodicity, which may play a key role in structuring information into cognitive maps for flexible navigation. A limitation of our cognitive model is that it solely accounts for afferent projections from the EC to the HPC, while the reverse projections from the HPC to the EC were not considered (Bonnevie et al., 2013; Dordek et al., 2016). Therefore, future biologically plausible models are needed to explain the recurrent mechanism of the EC-HPC circuit. Additionally, we encourage future neuroscience studies to validate the robustness of the 3-fold representation in the HPC, as well as to further investigate the hierarchical representation of activity periodicity.

Methods

Participants

Thirty-three right-handed university students with normal or corrected-to-normal vision participated in the experiment. Participants were recruited in three separate batches: 10 from Tsinghua University (Group 1: mean age = 23.9 [SD = 3.51], 5 females, 5 males), 10 from Peking University (Group 2: mean age = 22.10 [SD = 2.23], 5 females, 5 males), and 13 from Tsinghua University (Group 3: mean age = 20.85 [SD = 3.08], 2 females, 11 males). All participants had no history of psychiatric or neurological disorders and gave their written informed consent prior to the experiment, which received approval from the Research Ethics Committee of both universities.

Greeble space

Artificial Greeble space. The Greebles (Gauthier and Tarr, 1997) were generated based on the modified codes of the TarrLab stimulus datasets (http://www.tarrlab.org/) using Autodesk 3ds Max (version 2022, http://www.autodesk.com). A Greeble’s identity was defined by the lengths of two distinct features, represented by two nonsense words, “Loogit” and “Vacso” (Fig. 1a). These names were arbitrarily assigned and unrelated to the object-matching task. The generated Greebles were arranged in a 45-by-45 matrix, forming a two-dimensional space, with each Greeble representing a “location”. The feature ‘Loogit’ corresponded to the y-axis, while the feature ‘Vacso’ represented the x- axis (Fig. 1c). The prototype Greeble, positioned at the center of the Greeble space (Fig. 1c; blue dot), served as the target (i.e., “goal location”), with Loogit and Vacso having equal lengths. The length ratios of Greebles’ feature relative to the Greeble prototype ranged from 0.12 to 1.88, with a step size of 0.04, providing a smooth morphing experience as participants adjusted the lengths of the Greeble features.

Experimental design

The object matching task was programmed using Pygame (version 2.0, https://www.pygame.org/; Python version 3.9). Each trial began with a 2.0-s fixation screen, followed by a pseudo-randomly selected Greeble variant presented at the center of the screen (Fig. 1b). These variants were sampled from those near the boundary of the feature space (Fig. 1c; orange dots). Participants were instructed to adjust the features of the Greeble variant by pressing buttons to match the Greeble prototype as quickly as possible. Two response boxes enabled participants to stretch or shrink the two features. These boxes were counterbalanced between the left and right hands of participants. To mitigate the impact of learning effects on BOLD signals during MRI scanning, participants completed a training task one day prior to scanning. On the training day, participants familiarized themselves with task through observation and completed at least 10 practice trials before the main experiment. The training task was self-paced, allowing participants to adjust Greeble features without time constraints and press the return key to end each trial. Feedback was provided during the first 7 sessions (e.g., “Loogit: 5 percent longer than the target”), while no feedback was given in the final session. Day 2’s MRI experiment procedure mirrored day 1, except that 1) participants were required to complete a trial within 10 seconds. A timer was presented on top of screen (Fig. 1b), and 2) no feedback was provided for all experimental sessions. Based on the present design, the “movements” from a Greeble variant towards the prototype can be interpreted as forming a conceptual direction. To ensure an even distribution of conceptual directions, we arbitrarily defined the “movement” from the East as 0°. For each 30° bin, 24 Greeble variants were pseudo-randomly selected from the boundary of the feature space with an angular precision of 1.25°, resulting in a total of 288 Greeble variants distributed in 8 experimental sessions. Each session contained 36 object matching trials and 4 lure trials. The lure trial presented a blank screen after the fixation for 10 seconds. The MRI experiment lasted approximately 1 hour. Participants were incentivized to perform accurately through the promise of exponentially increasing monetary compensation based on their performance. After the scan, participants were debriefed, shown their performance results, and asked to discuss their subjective feeling as well as any perceived reasons for any performance issues. Importantly, participants were unaware of the existence of the Greeble space or the spatial features such as the location and direction during the task.

Behavioral performance

The behavioral performance score was calculated for each trial, with the formula T − T^′ + E/0.04, where T and T′ represent the participants’ actual and the objectively optimal trajectory length, respectively. The former was calculated by the number of movement steps in pixels from the original trajectories, while the latter was based on the number of pixels in the shortest trajectory between the Greeble variants and Greeble prototype. The optimal trajectory length was then subtracted from participant’s actual trajectory length. The rationale behind this approach was to eliminate the trajectory length difference between cardinal directions (i.e., “East”, “South”, “West”, and “North”) and intercardinal directions (e.g., “Southeast”) caused by the restricted movements along the horizontal and vertical axis of the squared Greeble space. In such a case, trajectories along cardinal directions tend to be shorter, while those along intercardinal directions become longer. The term E/0.04 quantified the error size, defined as the distance in pixels between the participant’s ending locations (adjusted Greeble) and the goal location (Greeble prototype). The constant “0.04”, representing the step size of the Greeble feature ratio, converted the error units of E from ratio to pixel space. Therefore, superior performance was represented by shorter trajectory and smaller error, characterized by smaller scores.

fMRI data acquisition

Imaging data were collected using a 3T Siemens Prisma scanner, equipped with a 64-and 20-channel receiver head coil at Tsinghua University and Peking University Imaging Center, respectively. Functional data were acquired with a Multi-band Echo Planer imaging (EPI) sequence with the following parameters: Acceleration Factor: 2, TR = 2000 ms, TE = 30 ms, matrix size = 112 × 112 x 62, flip angle = 90°, resolution = 2 × 2 × 2.3 mm³, 62 slices with a slice thickness of 2 mm and a gap of 0.3 mm, in a transversal slice orientation. Additionally, high-resolution T1-weighted three-dimensional anatomical datasets were collected for registration purposes using MPRAGE sequences, detailed as follows: TR = 2530 ms, TE = 2.98 ms, matrix size = 448 x 512 x 192, flip angle = 7°, resolution = 0.5 × 0.5 × 1 mm³, 192 slices with a slice thickness of 1 mm, in the sagittal slice orientation. At Tsinghua University, stimuli were presented through a Sinorad LCD projector (Shenzhen Sinorad Medical Electronics) onto a 33-inch rear-projection screen. At Peking University, stimuli presentation was managed using a visual/audio stimulation system (Shenzhen Sinorad SA-9939) and an MRI-compatible 40-inch LED liquid crystal display, custom-designed by Shenzhen Sinorad Medical Electronics Co., Ltd. The screen resolution for both sites was 1024 x 768. Participants viewed the stimuli through an angled mirror mounted on the head coil.

fMRI data preprocessing

The BOLD signal series of each scanning session was preprocessed independently using the FSL FEAT toolbox from the FMRIB’s Software Library (version 6.0, https://fsl.fmrib.ox.ac.uk/fsl/fslwiki)(Smith et al., 2004; Woolrich et al., 2009; Jenkinson et al., 2012). For each scanning sessions, the BOLD signals were corrected for motion artifacts, slice time acquisition differences, geometrical distortion using fieldmaps, and were applied with a high-pass filter with a cutoff of 100 seconds. Spatial smoothing was performed using a Gaussian kernel with a full width at half maximum (FWHM) of 5 mm. For group-level analysis, the preprocessed BOLD signals were registered to each participant’s high-resolution anatomical image and subsequently normalized to the standard MNI152 image using FSL FLIRT (Jenkinson and Smith, 2001). During normalization, the functional voxels were resampled to the resolution of 2 x 2 x 2 mm.

Sinusoidal analysis

Sinusoidal analysis was employed to investigate the 6-fold hexagonal activity in the EC following previously established protocols (Doeller et al., 2010; Constantinescu et al., 2016; Bao et al., 2019; Wagner et al., 2023; Raithel et al., 2023). First, the grid orientation in each voxel within the EC from each participant was calculated using half of the dataset (odd-numbered sessions: 1, 3, 5, 7). Specifically, a GLM was created for each session to model the BOLD signals of the EC with two parametric modulators: sin(6θ) and cos(6θ), convolved with the Double-Gamma hemodynamic response function. θ represents the movement directions, derived from drawing a line between the starting and the ending location. The factor “6” represents the rotationally symmetric 6-fold neural activity. The factors “3”, “4”, “5” and “7”, misaligned with the primary axes of grid cells, were used as control parameters. The estimated weights β_sine and β_cosine were then used to compute the grid orientation φ, ranged from 0° to 59°, where . The φ maps of each participant were averaged across the voxels of a hand-drawn bilateral EC mask and odd-numbered sessions to represent the grid orientation. Second, the neural representation of hexagonal activity, modulated by movement directions in Greeble space, were examined by another GLM using the other half of the dataset (even- numbered sessions: 2,4,6,8), with movement direction calibrated by grid orientation (even-numbered sessions: 2, 4, 6, 8). The GLM was specified with a cosine parametric modulator, cos(6[θ - Φ]). Third, the hexagonal activity maps were averaged across the even-numbered sessions for each participant before submitting to second-level analysis. To present the hexagonal pattern of EC activity in 1D directional space, we reconstructed the BOLD signals within the EC using parameters estimated from sinusoidal model for each participant.

Spectral analysis of MRI BOLD signals

Spectral analysis was employed to examine the spatial periodicity of BOLD signals. A GLM was created to extract the activity maps of movement directions. Participants’ movement directions were down-sampled into 10° bins (e.g., the movement directions ranging from 0° to 10° were labeled as 0°; and those from 350° to 360° were labeled as 350°). This process generated in a vector of 36 resampled movement directions (i.e., 36 directional bins), providing sufficient sample size for spectral analysis (De Valois et al., 1982). In the GLM, we modelled the BOLD signals with the 36 directional bins as regressors, yielding a vector of direction-dependent activity maps for each participant. These activity maps were sorted in ascending order (from 0° to 360°), detrended and processed with a Hanning window before being transformed from the spatial to the frequency domain using the Fast Discrete Fourier Transform (FFT) (implemented using the stats package of R platform, version 4.0, https://www.r-project.org/)(Landau and Fries, 2012). The spectral magnitude maps were generated from the absolute values of the complex FFT outputs for each of 0- to 18-fold and each participant. The two-tailed 95th percentile of the pseudo t-statistic distribution, computed as the mean spectral magnitude divided by the standard error across participants, was used as the initial threshold for the family-wise error correction in multiple comparisons.

Anatomical masks

The mask of the bilateral EC was manually delineated on the MNI T1 brain with 1 mm resolution packaged by FSL, using established protocols (Insausti et al., 1998), and the delineation software ITK-SNAP (Version 3.8, www.itksnap.org). The EC mask was then resampled to a 2 mm resolution. The anatomical mask of the bilateral HPC and the Medial Temporal Lobe (MTL) were derived from the AAL atlas, an automated anatomical parcellation of the spatially normalized single-subject high-resolution T1 volume provided by the Montreal Neurological Institute (Version: AAL 3v2; https://www.gin.cnrs.fr/en/tools/aal/)(Rolls et al., 2020).

Spectral analysis of human behavior

FFT analysis was employed to examine participants’ behavioral periodicity, with the procedure mirrored from detecting the brain activity periodicity. A performance vector was created from movement directions for each participant. These vectors were further down-sampled using 10° bins, sorted in ascending order (from 0° to 360°), detrended and processed with a Hanning window, resulting in a resampled vector of 36 sample points before sending to the FFT. The two-tailed 95th percentile of the spectral magnitude was used as the initial threshold for the family-wise error correction in multiple comparisons.

The “EC-HPC PhaseSync” model

The spatial code G of entorhinal grid cell was simulated using the cosine grating model (O’keefe and Burgess, 2005; Burgess et al., 2007; Blair et al., 2007; Bush and Schmidt-Hieber, 2018) in Python (version 3.9). At each location r = (x, y) of the Greeble space, three cosine gratings, oriented 60° apart, were generated by rotating around the anchor point r according to the rotation vector k, and then linearly summed (Equation 1). The parameter A represents the amplitude of the cosine gratings, C denotes the spatial phase at coordinates x and y, and C indicates the scale ratio (i.e., grid modules). This procedure generated a 45-by-45 grid cell population, with a constant grid orientation maintained across the population.

The trajectory code V of grid cell population excited in each movement was defined as the linear summation of grid codes across locations along the trajectory (Equation 2), where i and k denote the starting location (i.e., Greeble variant) and the ending location (Greeble prototype).

V exhibited a “planar wave” pattern (Welday et al., 2011; Krupic et al., 2012) when the movement direction φ_r aligned with the primary axes of grid cells (i.e., orientation ; aligned conditions). In contrast, when the movement direction φ_r deviated from the grid axes (i.e., orientation ; misaligned conditions), irregular multi-peak patterns were emerged from the trajectory code V. Note that the trajectory code V is exactly identical between the movement direction φ and φ + 180°, where φ was defined by arctan2 given the starting location r and ending location (Equation 3). In other words, the Greeble space is tessellated by a vector of orientations over Ψ ∈ [0, π) in the orientation domain.

The hippocampal representations of spatial orientation were defined as a vector of δ values, denoted as δ = [δ₁, δ₂,…, δ_n], where each δ_i was derived from the trajectory code V by linearly summing the EC activity across the location vector R from the trajectory code V_R,Ψ (Equation 4). The magnitude of each δ_i indexed the degree of alignment between the orientation Ψ and the grid axes, with larger values indicating closer alignment. In the orientation domain, the δ vector exhibited a 3-fold periodicity, characterized by a repeating “aligned-deviated” pattern driven by the periodic structure of the trajectory codes V(Ψ) projected from grid cell population, with the spatial phase δ_ε inherited from grid orientation.

To construct the vectorial representation C of the HPC for driving movements towards the goal, the spatial directions Φ between the location vector R and the goal location were embedded by the δ periodicity, with the spatial fold and spatial phase denoted by δ_c and δ_ε, respectively (Equation 5). C was then defined as the linearly summation between δ and a Gaussian-based distance term centered at the goal location q (Equation 5). The distance term is invariant to spatial directions but encodes spatial proximity to the goal. The parameter σ represents the spatial spread, defined by the radius of the Greeble space.

Statistical analysis

The uncorrected clusters derived from sinusoidal analysis and BOLD signal periodicity were derived using an initial threshold of p < 0.05 based on a two-tailed t-test. The reliability of these uncorrected clusters was assessed by the cluster-based permutation using FSL randomise (version 2.9, http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Randomise), which is a non-parametric statistical inference, and it does not rely on assumptions about data distribution (Nichols and Holmes, 2002). Specifically, 5,000 random sign-flips were performed on the sinusoidal beta and the FFT magnitude images. Clusters exceeding 95% of the maximal suprathreshold cluster size from the permutation distribution were considered significant. The uncorrected behavioral periodicity, derived from human participants and the EC-HPC PhaseSync model, was initially thresholded at p < 0.05 two tailed based on the shuffled magnitude distribution. The spectral magnitudes were corrected by the maximum of initial thresholds across the spatial folds from 0- to 18-fold for multiple comparisons. The EC-HPC activity coherence and the behavior-HPC activity coherence were evaluated using the amplitude-phase coupling analysis (Canolty et al., 2006) and the phase-lag index (Stam et al., 2007), respectively. For both analyses, the raw signals in the direction domain were band-pass filtered for the 3- and 6-fold using a two-way least-squares FIR filter (eegfilt.m from the EEGLAB toolbox; Delorme & Makeig, 2004). To evaluate the robustness of the coherence strength, a set of surrogates were created by offsetting the original signals by all possible time lags. The actual coherence strength was then compared with the mean of surrogates using ANOVA and t- test (Canolty et al., 2006).

Data and code availability

The analyses in the present study were carried out using custom scripts written in Python (version 3.9), Pygame (version 2.0), Matlab (version 2019b), R (version 4.0), and GNU Bash (version 3.2.57). The neuroimaging and machine learning toolkits employed include FreeSurfer (version 7.1, https://surfer.nmr.mgh.harvard.edu/), FSL (version 6.0, https://fsl.fmrib.ox.ac.uk/fsl/fslwiki), and Torch (version 1.9, https://pytorch.org/). The MRI dataset and behavioral data can be accessed via the Science Data Bank (https://www.scidb.cn/s/NBriAn). The original code for the “EC-HPC PhaseSync” model is openly available at (https://github.com/ZHANGneuro/The-E-H-PhaseSync-Model). Additional scripts are available upon request from the corresponding author (J.L.).

Supplemental figures

Fig. S1.

Fig. S2.

Fig. S3.

Fig. S4.

Fig. S5.

Fig. S6.

Fig. S7.

Acknowledgements

We sincerely thank Yue Wu, Yuannan Li, and Ao Li for their thoughtful discussions that greatly contributed to our data analysis. We are also grateful to Russell Epstein for his insightful feedback on the fMRI analysis and computational modeling, and to Dr. Michael Tarr for generously sharing the code used for generating the Greeble stimuli. This work was supported by the High-Performance Computing Platform of Peking University and the Resnick High Performance Computing Center of California Institute of Technology.

Additional information

Funding

The work was supported by the following funding sources: National Key R&D Program of China 2020AAA0105200 (J.L.), China Postdoctoral Science Foundation 2022M710470 (B.Z.), Beijing Municipal Science & Technology Commission & Administrative Commission of Zhongguancun Science Park Z221100002722012 (J.L.), Tsinghua University Guoqiang Institute 2020GQG1016 (J.L.), and Beijing Academy of Artificial Intelligence (J.L.).

Author contributions

Conceptualization: J.L.

Methodology: B.Z.

Investigation: B.Z., X.G.

Visualization: B.Z.

Resources: D.M., J.L.

Data Curation: B.Z.

Writing—Original Draft: B.Z., J.L.

Writing—Review & Editing: B.Z., J.L.

Project Administration: J.L.

Funding Acquisition: J.L.

Supervision: J.L.

Funding

National Key R&D Program of China (2020AAA0105200)

China Postdoctoral Science Foundation (2022M710470)

Beijing Municipal Science & Technology Commission & Administrative Commission of Zhongguancun Science Park (Z221100002722012)

Tsinghua University Guoqiang Institute (2020GQG1016)