Introduction

The hippocampus is known as its central role in encoding and storing experience and thus is essential for forming episodic memories (Buzsaki, McKenzie, & Davachi, 2022; Cohen & Eichenbaum, 1993). In the context of active navigation in an environment, hippocampal place cells are organized to fire sequentially within individual cycles of theta oscillations (4-12Hz) and form time-compressed representations of behavioral experiences, referred to as theta sequences (Colgin, 2020; Comrie, Frank, & Kay, 2022; Dragoi & Buzsaki, 2006; Skaggs, McNaughton, Wilson, & Barnes, 1996; Wallis, 2018). By sweeping ahead from the past to the upcoming locations, theta sequences represent information about potential future paths that could guide goal-directed behaviors (Gupta, van der Meer, Touretzky, & Redish, 2012; Joshi et al., 2023; Zheng, Hwaun, Loza, & Colgin, 2021). These predictive theta sequences exhibited experience-dependent development when rats were exploring a novel environment or learning a novel goal location (Feng, Silva, & Foster, 2015; Gobbo et al., 2022; Wikenheiser & Redish, 2015; Zheng et al., 2021). However, how the hippocampal network modulated the development of theta sequences’ temporospatial-compressed organization remains incompletely understood.

A common assumption was established by previous studies that theta sequences were manifestations of theta phase precession, in which an individual place cell fired at progressively earlier theta phase as the animal traversed the place field (Aoki et al., 2023; Foster & Wilson, 2007; O’Keefe & Recce, 1993; Qasim, Fried, & Jacobs, 2021; Skaggs et al., 1996). However, when examined the theta sequences with experience at single trial level, their sequential structures were not developed on the first experience of a novel environment even if the phase precession of individual place cells has been observed (Colgin, 2020; Feng et al., 2015). This suggests that the dynamic process of theta sequences development could be modulated by another network-level mechanism with neuronal firing organized at a finer timescale within theta cycles. Gamma (25-100Hz) could be a candidate that is thought to interact with theta to temporally organize theta sequences (Dragoi & Buzsaki, 2006; Fernandez-Ruiz, Sirota, Lopes-Dos-Santos, & Dupret, 2023; Lisman & Jensen, 2013). As two subtypes of gamma, fast and slow gamma rhythms (Colgin et al., 2009; Schomburg et al., 2014; Zhang, Lee, Rozell, & Singer, 2019) occur at distinct theta phases and coordinate theta sequences in different manners (Zheng, Bieri, Hsiao, & Colgin, 2016). This indicates that fast and slow gamma may act as different tuners but collaborate together in modulating theta sequence development.

Fast gamma (∼65-100Hz), associated with the input from the medial entorhinal cortex, is thought to rapidly encode ongoing novel information in the context (Fernandez-Ruiz et al., 2021; Kemere, Carr, Karlsson, & Frank, 2013; Zheng et al., 2016). During fast gamma rhythms, place cell spikes occurred across all positions within a place field and displayed theta phase precession (Bieri, Bobbitt, & Colgin, 2014). This type of cells was identified recently as “phase-precessing” cells, in which the phase-precession pattern emerges with strong fast gamma (Guardamagna, Stella, & Battaglia, 2023). Their discrete distribution of theta phases is likely to contribute to the sequential structure formation when the sequences have been developed (Wang, Foster, & Pfeiffer, 2020). However, a question is raised about how these fast gamma-modulated cells be coordinated for the developmental process of theta sequences.

Slow gamma (∼25-45Hz) in the hippocampal CA1 reflects information received from CA3, which is responsible for the integration processing of learned information and experience (Colgin et al., 2009; Zhu et al., 2023). In contrast to fast gamma, spiking during slow gamma exhibited dominant theta phase locking and attenuated theta phase precession (Guardamagna et al., 2023; Zhang et al., 2019). This concentrated distribution of theta phase could be interpreted as they are coordinated within slow gamma cycles, as a finer time-scale than a theta cycle. Place cells spiking exhibited slow gamma phase precession, to generate mini-sequences in an individual theta sequence, for highly time-compressed representation of ongoing and future information (Zheng et al., 2016; Zheng et al., 2021). Thus, how fast and slow gamma rhythms coordinate place cells spiking to modulate the dynamic process of theta sequence development remains unclear.

Here, we proposed a scheme that fast gamma rhythms coordinate a subgroup of place cells for rapidly encoding and processing sensory information since from the early phase of theta sequence development (Fig.1). The developing process would be then dependent on the slow gamma modulation, that the spikes are fired at slow gamma phase-precessed pattern for spatial information compression within theta cycles. To test this hypothesis, we investigated how fast and slow gamma separately modulated place cells activities during the whole process of theta sequence development. We found that the fast gamma phase-locked place cells (FG-cells) were crucial for the development of theta sequences, without which the temporospatially compressed structure of theta sequences could not be developed. Furthermore, these cells exhibited slow gamma phase precession within theta cycles associated with the sweep-ahead structure of theta sequences, which occurred earlier than those cells not phase locked to fast gamma (NFG-cells). Our findings highlight the dynamic network-modulation by neural oscillations fast and slow gamma, at finer timescales within theta, in development of theta sequences which may further facilitate memory encoding and retrieval.

Fig. 1

Results

A subgroup of hippocampal place cells phase locked to fast gamma rhythms during active behaviors

To investigate the temporal organization of place cells’ firing associated with the gamma rhythms during the development of theta sequences, we recorded neuronal spiking of place cells and local field potentials simultaneously in hippocampal CA1 from 4 feely behaving rats (Fig.2-S1). As shown in a previous study, 32% of hippocampal pyramidal cells were phase-locked to gamma rhythms during waking theta states (Senior, Huxter, Allen, O’Neill, & Csicsvari, 2008). We further tested how the hippocampal place cells were modulated by slow and fast gamma rhythms respectively. The slow and gamma episodes were detected in time periods when the power of each gamma type (Fig.2A-B, fast gamma: 65-100 Hz, slow gamma: 25-45 Hz) was above a preset threshold (see Materials and Methods). We then detected the slow or fast gamma phase-locked place cells by quantifying the gamma phase distributions of all spikes across all trials. In this study, 488 place cells were included with place field distributed on the circular track but not located at the resting box. We found that 120 (24.6%) out of these cells were significantly phase-locked to gamma rhythms, a majority of which (113 cells, 23.2%) exhibited significant phase-locking to fast gamma rhythms (Fig.2C). Their preferred fast gamma phases (202.38±2.01 deg) were concentrated around the peak of fast gamma cycle, which was consistent with the findings in Schomburg et al. 2014 (Schomburg et al., 2014). However, a very small number (3.7%) of cells were found to phase-locked to slow gamma rhythms. Thus, these place cells were subdivided according to whether they were significantly phase-locked to fast gamma (FG-cells) or not (NFG-cells) during active behaviors in the subsequent analyses. Furthermore, the fast gamma phase-locking of FG-cells was stable across trials, with a significantly higher mean vector length of fast gamma phase-locking than those of NFG-cells (Fig.2D, repeated measure ANOVA, main effect of cell type, F(1,22)=33.9, p=7.5×10^-6, partial η² =0.6; main effect of trials, F(4,88)=1.3, p=0.3, partial η² =0.06, no cell type x trials interaction, F(4,88)=0.5, p=0.8, partial η² =0.02). Meanwhile, we also found that a larger proportion of neurons in FG-cells had significant theta phase precession than NFG-cell (Fig.2-S2, chi-squared test, ratio of theta phase precession cell among the FG-cell vs. NFG-cell, n=457cell, χ²=13.7, p=2.1×10^-4, Cramer V=0.2). These results suggested that the FG-cells were stably modulated by fast gamma rhythms during active behaviors.

Fig. 2

The FG-cells were crucial for theta sequence development

As theta sequences were developed dependent on experience, we next investigated whether FG-cells and NFG-cells played different roles in the development of theta sequences. We then used a Bayesian decoding approach (see Methods) to identify trajectories represented as rats running on a circular track, and detected significant theta sequences in each single lap (Fig.3A). Another two decoders were built by excluding either FG-cells spikes or down-sampled NFG-cells spikes respectively, in order to compare the structures of theta sequences without each type of cells. It can be observed that posterior probability distribution of theta sequences was changed by excluding spikes from FG-cells (Fig.3B). The weighted correlation of each sequence (Feng et al., 2015) was then measured to quantify the temporospatial correlation of a theta sequence structure. We found that the weighted correlation of theta sequences decoded by excluding FG-cells spikes (exFG-sequences) was much lower than that of theta sequences decoded by excluding down-sampled NFG-cells spikes (exNFG-sequences) in the same theta cycles (Fig.3C, repeated measure ANOVA: main effect of cell type, F(2,4844)=128.7, p=3.5×10^-48, partial η² =0.05; post-hoc test, raw sequences vs. exFG-sequences, p=8.0×10^-56, raw sequences vs. exNFG-sequences, p=1.6×10^-26, exFG-sequences vs. exNFG-sequences, p=4.1×10^-9). This suggests that the sequential structure could be dominantly interrupted by excluding FG-cells rather than NFG-cells.

Fig. 3

Consistent with the previous studies (Feng et al., 2015; Wang et al., 2020), we also observed the development of theta sequences along with exploration experience (Fig.3D-E, repeated measure ANOVA, main effect of trials, F(4,44)=7.0, p=2.0×10^-4, partial η²=0.4). This could not be accounted for by running speed, as the rats were feely moving at stable running speed across trials (Fig.3F, repeated measure ANOVA, main effect of trials: F(4,44)=0.9, p=0.5 partial η² =0.07). We further investigated whether the FG-cells or NFG-cells were required for the development of theta sequence structure. Interestingly, the temporospatial structure of the theta sequences cannot be developed with experience when the FG-cell spikes were excluded, with significantly lower weighted correlation than theta sequences constructed by excluding NFG-cell spikes (Fig.3D,G, repeated measure ANOVA, trials x cell type interaction, F(4,44)=3.0, p=0.03, partial η²=0.2, main effect of cell type, F(1,11)=7.8, p=0.02, partial η²=0.4; main effect of trials, F(4,44)=4.0, p=0.008). Since the theta sequences have been developed after Lap2, we defined the early development phase of sequences as the first 2 laps and late development phase of sequences as the last 2 laps of exploration experience, respectively. We found that the weighted correlation of both exFG-sequences and exNFG-sequences was low during early development (Fig.3H, paired t-test, t(24)=0.2, p=0.9, Cohen’s d=0.04). During late development, however, the weighted correlation of exNFG-sequences was higher than that of exFG-sequences (Fig.3H, paired t-test, t(24)=- 3.7, p=0.001, Cohen’s d=0.8). These results suggest that fast gamma modulation of place cells may contribute to the development of theta sequence with exploration experience.

In order to estimate to what extent, the FG cells contributed to theta sequence development, we investigated the relationship between the percentage of FG-cells and the formation of sequential structure in different developments. In our dataset, the percentage of FG-cells out of all place cells varied across days, thus we detected theta sequences from trials with relatively low percentage (∼10%) and high percentage (∼35%) of FG-cells. The number of excluded NFG-cells for decoding was matched up with the number of excluded FG-cells by a pseudo-random method in each day. During early development, the temporospatial structure of theta sequences has not been developed by excluding either low or high percentage of each type of cells (Fig.4A). During late development, however, the sweeping-ahead structure of the theta sequence was significantly disrupted when excluding a high percentage (∼35%) of FG-cells, which was not observed when excluding a low percentage (∼10%) of FG-cells (Fig.4B). The sequential structure could be maintained by excluding NFG-cells at either low or high percentage. To quantify the relationship between sequence structure and the number of excluded neurons during early and late development, we then measured the weighted correlation difference between exFG-sequences and exNFG-sequences as a function of the fraction of excluded cells. The results showed that the weighted correlation difference was close to 0 regardless of the fraction of excluded cells during early development (Fig.4C, Pearson’s correlation, Corr=0.1, p=0.7). During late development, however, we found that the weighted correlation difference was positively correlated with the fraction of excluded cells (Fig.4D, Pearson’s correlation, Corr=0.6, p=0.03). These results suggest different roles of FG-cells and NFG-cells in the development of the sweeping-ahead structure of theta sequence with experience, that the more temporospatially compressed structure of theta sequences may require a larger percentage of FG-cells.

Fig. 4

In addition, we tested whether these results could be accounted for by the different firing characteristics between FG-cells and NFG-cells. We found that although these two types of cells took different proportion of all recorded place cells, the place fields of both types of cells could cover the entire track uniformly (Fig.2C, Fig.5A). We measured the distribution of center of mess (COM) of main place field of FG-cells and number-matched NFG-cells (Fig.5B), and found similar distribution of place field COM between cell types (Fig.5C, two-sample Kolmogorov-Smirnov test, Z=0.1, p=0.8). Furthermore, we found that the theta sequences in the above analyses were decoded by excluding similar numbers of spikes from FG-cells and NFG-cells (Fig.5D-E, Student’s t test, spike number: t(22)=1.8, p=0.09, Cohen’s D=0.7; mean firing rate in field: t(22)=2.0, p=0.06, Cohen’s D=0.8), indicating that the deferent development patterns between exFG-sequences and exNFG-sequences were not due to quality difference between two decoders. The two type of cells exhibited similar place fields size (Fig.5F, t(22)=0.7, p=0.5, Cohen’s D=0.3) and spatial information (Fig.5G, t(22)=1.0, p=0.3, Cohen’s D=0.4). These findings support the hypothesis that FG-cells, but not NFG-cells, may contribute to the development of theta sequences with experience.

Fig. 5

FG-cells exhibited slow gamma phase precession for information compression required in theta sequence development

In our previous study, we found that the place cell spikes occurred at progressively earlier slow gamma phases across slow gamma cycles in theta sequences, so-called slow gamma phase precession (Zheng et al., 2016). As sequences of locations were represented at different slow gamma phases, theta sequences exhibited temporospatially compressed structure representing long paths. We next investigated if the development of theta sequence structure was related to the slow gamma phase modulation of place cells. By detecting the slow gamma phases of spikes within all theta cycles, we observed that the slow gamma phases significantly shifted across successful slow gamma cycles (Fig.6A, multi-sample Mardia-Watson-Wheeler test, W(4)=369.4, p<2.2×10^-16; post-hoc test, Cycle -1 vs. Cycle 0, W(2)=97.5, p<2.2×10^-16, Cycle 0 vs. Cycle 1, W(2)=95.4, p<2.2×10^-16), which was consistent with the previous finding (Zheng et al., 2016). Then we separated out the spike activity from the early and late development, and tested if the slow gamma phase precession still existed during different developments. Interestingly, we observed slow gamma phase shifting, but not phase precession, across successive gamma cycles during early development. The distribution of slow gamma phases in Cycle0 was similar as that in Cycle-1, then was dispersed in Cycle1 (Fig.6B, multi-sample Mardia-Watson-Wheeler test, W(4)=18.5, p=9.8×10^-4; post-hoc test, Cycle -1 vs. Cycle 0, W(2)=0.8, p=0.7; Cycle 0 vs. Cycle 1, W(2)=15.5, p=4.4×10^-4). However, the slow gamma phase precession existed during late development (Fig.6C, multi-sample Mardia-Watson-Wheeler test, W(4)=21.2, p=2.8×10^-4; post-hoc test, Cycle -1 vs Cycle 0, W(2) = 8.8, p = 0.01; Cycle 0 vs. Cycle 1, W(2)=8.2, p=0.02). This finding suggests that slow gamma phase precession of place cells may occur with exploring experience.

Fig. 6

Furthermore, we investigated how the slow gamma phase precession occurred in theta sequences developed with experience. We detected theta sequences during slow gamma episodes and measured the distribution of phase shifts across gamma cycles within single unites. In the significant theta sequences during slow gamma episodes, we observed the slow gamma phase precession from single unit, that the place cell spikes occurred at earlier slow gamma phases in later slow gamma cycles (Fig.7A). To quantify the slow gamma phase precession at single unit level, we calculated the phase shift as the phase difference within spike pairs from adjacent active slow gamma cycles (See Methods). We found that slow gamma phase difference was significantly lower than 0 between 2 active slow gamma cycles (Fig.7B, one sample t-test, t(485)=-11.9, p=6.3×10^-29, Cohen’s D=0.5), as well as among 3 or ≥4 active slow gamma cycles for strict limitation of theta sequences modulated by slow gamma (3 active cycles, t(140)=-6.8, p=3.0×10^-10, Cohen’s D=0.6; ≥4 active cycles, t(27)=-3.2, p=0.003, Cohen’s D=0.6). Furthermore, we averaged the slow gamma phase differences across different cells within a theta sequence and compare their distribution along with learning trials. In the first lap, we found that the slow gamma phase differences were not significantly different from their surrogate data (Fig.7C, median of phase difference from lap1 was in 95% confidence intervals of those from shuffled data, p=0.3). However, the distribution of slow gamma phase difference was negatively biased compared with their surrogate data (Fig.7C, median of phase difference from lap2-5 was beyond the 95% confidence intervals of those from shuffled data; lap2 p=0.001, lap3 p=0.001, lap4 p=0.002, lap5 p=0.002). These results can also be replicated by averaging the slow gamma phase differences from a single unit across succussive theta sequences (Fig.7D, lap1 p=0.1, lap2 p=0.001, lap3 p=0.001, lap4 p=0.002, lap5 p=0.01).

Fig. 7

These findings raised a question of whether the development of slow gamma phase presession existed for both FG-cells and NFG-cells. To address this question, we quantified the distribution of slow gamma phase difference in FG-cells and NFG-cells respectively. We found that the slow gamma phase difference of FG-cells consistently exhibited a stable phase precession across trials, independent with learning experience (Fig.7E, lap1 p=0.007, lap2 p=0.006, lap3 p=0.007, lap4 p=0.002, lap5 p=0.02). In contrast, the slow gamma phase precession of NFG-cells was experience-dependent, with negative biased distribution of phase difference occurring from Lap2 (Fig.7E, lap1 p=0.9, lap2 p=0.01, lap3 p=0.001, lap4 p=0.009, lap5 p=0.1). These results suggest that FG-cells and NFG-cells have different dynamic modulated by slow gamma rhythms during exploration experience. The slow gamma modulation of FG-cells may occur earlier than that of NFG-cells, thereby driving the location information to be coded within theta cycles in a highly temporospatially compressed way.

FG-cells and NFG-cells exhibited different dynamics of slow gamma modulation along with theta sequence development

Finally, we wondered if there was correlation between slow gamma phase precession and the development of theta sequence structure, and how the FG-cells and NFG-cells coordinated during this process. Within each theta sequence, we measured its weighted correlation of the sequence structure, and plot it as a function of slow gamma phase difference across slow gamma cycles. Interestingly, we found most data were distributed in the Quadrant1 (Q1) compared to the shuffled data distributed at the center (Fig.8A-B). The theta sequences with negative slow gamma phase difference and positive weighted correlation were significantly dominant (Fig.8C, the probability of observed data from lap1 exceeded 95% confidence intervals of those from shuffled data, p=0.001). This finding supports our hypothesis that the theta sequence development could be associated with the slow gamma phase precession.

Fig. 8

We further investigated the contribution of FG-cells and NFG-cells to the slow gamma modulated theta sequence development during different exploring phases. Importantly, the distribution of FG-cell-dominant theta sequences was significantly biased to Q1 than that of the shuffled data, not only in late development when the theta sequences have been developed, but also in early development (Fig.8D-F, early development, p=0.001; late development, p=0.004). However, for NFG-cell-dominant theta sequences, their distribution was significantly biased to Q1 than that of the shuffled data only in late development (Fig.8G-I bottom, p=0.046), but not in early development (Fig.8G-I top, p=0.08). Taken together, these results support the hypothesis that FG-cells and NFG-cells may exhibit different dynamical characteristics in gamma modulation, thus their coordination may contribute to the development of theta sequences. FG-cells were consistently modulated by slow gamma phase precession, which would be precondition for theta sequence development.

Discussion

In current study, we employed a circular track exploration task to investigate the different roles of fast and slow gamma rhythms in the coordination of neurons during the experience-dependent development of theta sequences. We first found that a subgroup of place cells (FG-cells) that were phase-locked to fast gamma played a crucial role in the development of theta sequences. By excluding a sufficiently high proportion of FG-cells, the sweep-ahead structure of theta sequences was significantly disturbed compared to that of excluding equivalent number of NFG-cells. These FG-cells exhibited slow gamma phase precession within single theta cycles throughout the developing process of theta sequences. Finally, we showed positive correlation between the intensity of slow gamma phase precession and development of the sweep-ahead structure of theta sequences.

Fast gamma episodes were detected with dominantly high probability as initial exploration of each day and decreased probability as being familiar to the environment (Bieri et al., 2014). This could be related to the feedforward communication from upstream to downstream regions by using fast gamma frequencies (Vinck et al., 2023). Although the information propagation of visual stimulation was unlikely synchronized between visual cortex and hippocampal CA1 (Schneider, Tzanou, Uran, & Vinck, 2023), the MEC inputs could also be entrained by fast gamma to process ongoing sensory information to the hippocampus (Buzsaki & Schomburg, 2015; Colgin et al., 2009). Thus, the fast gamma phase-locked ensemble may act as a receiver of novel information that organized its intrinsic individual neurons in a time-compressed and spatially-ordered sequence to integrate the information under the mesoscopic cellular substrate modulated by fast gamma rhythms (Fernandez-Ruiz et al., 2023). These neurons exhibited similar firing characteristics as those of fast gamma non-phase-locked neurons, however, they showed higher probability presenting theta phase precession. Therefore, this fast gamma phase-locked ensemble likely has overlapping with “bimodal cells” involved in the formation of theta sequence structure including reverse and forward windows (Chu et al., 2023; Wang et al., 2020). In this case, regardless of these neurons, the sweep-ahead structure of theta sequences could not be developed with learning, with the extent of structure impairment positively correlated with the proportion of excluded fast gamma phase-locked neurons.

The important role of fast gamma in modulating theta sequence development was also supported by a specific perturbation experiment. Liu and colleagues optogenetically entrained the MEC GABAergic cells at an artificial gamma frequency (53Hz) and observed the loss of phase precession patterns in CA1 place cells, and thereby the disrupted development of theta sequences (Liu, Todorova, Tang, Oliva, & Fernandez-Ruiz, 2023). The failure of sequential firing organization during behaviors impaired the formation of a predictive map related to the highly time-compressed awake replay sequences executing information storage. Fast gamma was also found diminished between MEC and CA1 in transgenic model mice of Alzheimer’s Disease, linked with spatial memory impairment in AD as a circuit mechanism(Jun et al., 2020). On the other hand, our previous study showed that fast gamma-dominated theta sequences represented shorter trajectories with relative less spatial information within a theta cycle (Zheng et al., 2016). This is indeed not contradictory with the current results, as we found that slow gamma, related to information compression, was also required to modulate fast gamma phase-locked cells during sequence development. We replicated the results of slow gamma phase precession at ensemble level (Zheng et al., 2016), and furthermore observed it at late development, but not early development, of theta sequences. Thus, a scenario could be drawn out that the network modulator was transitioned from fast gamma to slow gamma throughout the sequence development with learning.

At single neuron level, we also observed that spikes of an individual place cell occurred at precessed slow gamma phases across consecutive slow gamma cycles under the time-scale of a theta window. At a finer scale within a theta cycle, slow gamma cycles coordinated place cells in a “mini-sequence” way (Zheng et al., 2016), in order to integrate information of long spatial trajectory related to the memory retrieval and predictive behaviors. This negatively-biased phase difference distribution appeared from lap2, when the sweep-ahead structures of theta sequences haven’t been fully developed, suggesting that slow gamma phase precession of place cells contributed to the information compression required to the theta sequence development. More interestingly, we found that slow gamma phase precession of a subset of neurons was absent at initial exploration when the theta phase precession has already displayed (Feng et al., 2015). This finding was in line with a previous study showing that silencing of CA3 output didn’t affect theta phase precession of a single cell, but impaired theta sequence structures (Middleton & McHugh, 2016). A possibility was that recruitment of PV-interneurons during the CA3-driven descending phase of the theta cycle, associated with the induction of slow gamma oscillations in CA1, was found reduced in the absence of CA3 input (Fernandez-Ruiz et al., 2023; Topolnik & Tamboli, 2022). We also provided additional evidence about the relationship between slow gamma phase precession and theta sequence formation, that the more negative slow gamma phase difference across cycles, the more sweep-ahead structure represented at late development of theta sequences.

It is noteworthy that the subset of neurons (FG-cells) was modulated by both fast gamma and slow gamma rhythms during the theta sequence development (Fig.1). These place cells, associated with both memory encoding and retrieval, seem to function like memory engram-like representations (Goode, Tanaka, Sahay, & McHugh, 2020; Josselyn & Tonegawa, 2020; Lamothe-Molina et al., 2022; Pettit, Yap, Greenberg, & Harvey, 2022). Indeed, the proportion of c-Fos positive place cells (putative engram cells) had higher proportion being phase-locked to fast gamma than the c-Fos negative place cells (Tanaka et al., 2018). In contrast, the NFG-cells exhibited experience-dependent phase coding, that they likely contributed to the path representation at the late learning process. Thus, further studies would be needed to investigate how memory engram cells contribute to theta sequence development and facilitate spatial memory and navigation.

Materals and methods

Subjects

Four male Long-Evans rats weighing 500 to 700g were used in this study. Rats were housed in custom-made acrylic cages (40cm×40cm×40cm) and exposed to a light cycle opposite to natural light (light exposure time: 8:00pm-8:00am). All waking activity recordings were performed in the dark (recording time: 8:00am-8:00pm). Rats were handled at the time of electrode implantation surgery and underwent at least 7 days of behavioral pretraining. Behavioral training was initiated after at least one week of recovery from implantation surgery. Rats were physically deprived during experimental data collection so that their body weight was maintained at ∼90% of their free-feeding body weight. All experiments were conducted according to the guidelines of the Animal Care and Use Committee of Tianjin University.

Behavioral paradigm

We designed a circular track apparatus for this study. The circular track with an inner diameter of 90 cm, height of 50 cm, and width of 10 cm, was placing an acrylic resting box (height of 31 cm, length of 28 cm, and width of 21cm) as beginning and ending. Each rat underwent training to navigate on a circular track. A total of eighteen small black silicone dots were utilized as context cues along the circular track, covering approximately two-thirds of its length. Each dot was at a distance of 12 cm (or equivalently 0.28 radians) away from its neighbors. Before surgery, all rats were pre-trained on the same circular track for 1 week and recovered from surgery for at least one week before behavioral training resumed. For every recording session, rats ran 5 trials (or laps) on the circular track. For each running lap, the rats started from the rest box, navigating clockwise to the rest box, and then returned to the rest box counter-clockwise without receiving any food reward. A tone (1kHz sine wave with a duration of 0.25s) was delivered at the beginning and end of the trial. According to the recorded position, the 2D movement speed of the rat was calculated by 2D coordinates and elapsed times of consecutive two points.

Surgery and tetrode implanting

Before surgery, tetrodes were built from 17μm polyimide-coated platinum-iridium (10/90%) wire (California Fine Wire, Grover Beach, California). The tetrode tip was platinum-plated using a NanoZ system (White Matter, Seattle, Washington), ensuring that the impedance of each recording channel was below 300 kOhms to facilitate single-neuron activity recording. Eighteen independently mobile tetrodes hyperdrive were surgically implanted into the hippocampal dorsal CA1 region. The rat brain stereotaxic coordinates of the center of the electrode bundle were anterior-posterior [AP] -3.5 mm, medial-lateral [ML] 3.0mm and on day of implantation the tetrode tips moved 1mm downward (dorsal-ventral [DV] 1mm). Screws were placed on the surface of the skull. The hyperdrive bundle were attached to the screws and skull with dental cement. Two screws placed in the skull above the cerebellum served as a panel ground for recording. After surgery, tetrodes were gradually downward to dCA1 stratum pyramidale. One tetrode was designed as recording reference. Recording reference was used to differentiate from other recording tetrode to monitor signal During other recording tetrode downwards, the recording reference was moved up and down to a quiet position in the cortex above the hippocampus and remained quiet throughout signal recording session.

Data acquisition

Behavioral and electrophysiological data were acquired using a Digital Lynx system and Cheetah recording software (Neuralynx, Bozeman, Montana). LFPs recorded continuously from a single channel in each tetrode, within the frequency range of 0.1-500 Hz from raw signal at 2,000 Hz sampling rate. Before each recording session starting, the input amplitude ranges were carefully adjusted to optimize resolution and prevent signal saturation. Typically, the input ranges for LFPs were within the range of ±2,000 μV to 3000μV. To detect unit activity, signals from each channel in a tetrode were bandpass filtered from 600 to 6,000 Hz. The spike detected threshold was 55 μV. All the recorded signals were obtained by differentiating with the reference tetrode (see “surgery and tetrode planting”), eliminating a large amount of background noise. The position of the rat during navigation behavior was recorded by the video tracking function within Neuralynx system. The video resolution of 720 × 480 pixels and a frame rate of 25 frames per second. The actual position was tracked with a red light-emitting-diode (LED) on one of the three HS-27-Mini-LED headstage (Neuralynx, Bozeman, Montana) attached to the hyperdrive. The actual 2D coordinates were obtained according to the pixel value where the red is brightest

Spike sorting

Detected spikes were manually sorted using graphical cluster-cutting software (MClust; A.D. Redish, University of Minnesota, Minneapolis). Spikes were sorted using two-dimensional projections of three different features of spike waveforms (i.e., energies, peaks, and peak-to-valley differences) from four channels. Clusters were judged based on waveform shape, isolation distance, and lack of inter-spike-interval violation. Units with mean firing rates of more than 15 Hz were identified as putative fast-spiking interneurons and were excluded from further analysis. The units with peak firing rates more than 1 Hz but mean firing rates of less than 15 Hz was putative pyramidal cell for further analysis.

Spatial firing maps and place cell identification

The animal position on the circular track was linearized to increasing radian values and divided into 90 bins (approximately 0.035 rad per bin and that is 3.5cm). Spatial firing maps were calculated only during locomotor periods (> 5cm/s). To calculate directional spatial firing rate-maps, the number of directional spikes was divided by the time of occupancy at each position bin on the circular track with the same direction, followed by smoothing with a Gaussian kernel (standard deviation of 2 bins). For identifying the place cell, we used all trial spikes to calculate rate-maps. A continuous span of at least three consecutive bins in the rate-map with a peak firing rate exceeding 1 Hz was calculated as a place field. A pyramidal cell contained at least one place field was defined as a place cell. The place field size was defined by 0.1Hz firing rate boundary. The firing rate in each place field is divided by the number of bins in the field to calculate the mean firing rate in the field, the peak firing rate. The primary place field used for analysis was defined as the field with the maximal peak firing rate among all fields. Finally we accepted 488 place cells with place field distributed on the circular track (region from dot 1 to 18) but not located at the resting box.

Spatial specificity was quantified by computing the spatial information score (Langston, Ainge et al. 2010). For each cell, the spatial information score in bits per spike was calculated from the recording in each period as:

Where λ_i is the mean firing rate of a unit the ith bin, λ is the overall mean firing rate, and p_i is the probability of the animal being in the i th bin (occupancy in the i th bin/total occupancies bin).

The single-lap(or trial) rate-maps was calculated using the same method and used to Bayesian decoding analysis.

Local field potential analysis

Theta, slow gamma, and fast gamma rhythms were extracted from the LFP by applying band-pass filters with frequency ranges of 4-12 Hz, 25-45 Hz, and 65-100 Hz, respectively. The time-varying phases of the theta, slow gamma, and fast gamma rhythms were computed using the Hilbert transform. As the previous study showed, the gamma power spectra were measured across different running speeds. The power spectra across running speeds and absolute power spectrum (both result were not shown) was estimated through the multitaper spectral analysis (Mitra and Bokil, 2008) in the Chronux toolbox (http://chronux.org/). We defined the gamma episode as the 160ms window with gamma power higher than 3 standard deviations to observe the firing activity of place cells.

Quantification of theta phase precession

Theta phase precessions were examined in the primary place field calculated by all trial spikes. The phase-position correlation of spikes emitted by a place cell during all lap of its place field was determined by circular-linear regression between the theta phase of spikes and the actual position of animal. The significant negative correlation was determined theta phase precession. Theta phase precession slope was calculated by least-square fit of spike locations of the animal and shifted spike phases.

Estimation of gamma phase-locking

Phase-lock analysis was implemented using custom scripts based on the circ_stats function from Matlab CircStat toolbox. The slow gamma, and fast gamma spike phase distributions for each place cell were then determined by the time-varying phases of slow gamma, and fast gamma (see “Local field potential analysis”), respectively, at the time point closest to each spike time. And each spike phase was evaluated with the LFP recorded in the tetrode where the neuron was located. Phase-locking was quantified by calculating the mean vector length of the resulting phase distributions. Significant phase-locking at fast gamma is identified only for place cells meeting the criteria of a Rayleigh test p<0.05 and firing at least 10 spikes on each lap. For easy presentation, these cells are named FG-cells. The remaining place cell are not phase-locking at fast gamma, called NFG-cells. According to our results, nearly a quarter of the neurons are significantly phase-locked at fast gamma. Only a very small number of neurons are slow gamma phase-locked.

Single-trial Bayesian decoding analysis

To translate ensemble spiking activity into radian positions on the circular track, a memory-less Bayesian decoding algorithm was implemented. Decoding was performed using a 20ms sliding time window that shifted 5ms at each step. For each decoding window, the number of spikes occurring is denoted as n_i, so that the vector of all spiking activity simultaneously recorded for the cell is n = {n₁, n₂, ⋯, n_i}. The probability of a rat being at position x, given the number of spikes n from each unit recorded in a time window, was estimated by Bayes’ rule:

where P(n|x) was approximated from the single directional firing turning curve of each unit. The approximation assumed that the turning curve of individual units was statistically independent and the number of spikes from each unit followed a Poisson distribution. Prior knowledge of position P(x) was set to 1 to avoid decoding bias to any particular location on the track. The normalizing constant P(n) was set to ensure P(x|n), or posterior probability, summed to 1. The decoding position was is the set of x that maximizes P(x|n). Decoding error was defined as the distance between the decoded position with maximum-likelihood and the rat’s true position. For excluded cell decoding, position reconstruction was performed by removing the FG-cells or NFG-cells from n, respectively. Bayesian decoding analyses were performed using custom routines in MATLAB 2020b.

Theta cycle detection

A single tetrode located in the stratum pyramidale and physical center all recording tetrodes was selected and the theta signal (including phase and power etc.) from this tetrode was utilized for further analysis. According to previous studies, almost all theta sequences fall completely into a theta cycle. Therefore, the global 0 phase (the beginning of the theta cycle) of theta rhythm was determined as the phase when the theta sequence quantization indicators (see “Theta sequence quantification”) were maximal. Briefly, the global 0 phase was moved from the peaks in steps of 10°, and in every step, the sequence of place cells was quantified in each theta cycle. The peak of the histogram of the distribution of the quantization indicators over the step size was used as the global 0 theta phase for further analysis.

Theta sequence quantification

Theta sequences were quantified only when the following conditions also hold: (1) the animal was in the middle of the track (two-thirds of the circular track length), (2) the current running direction is consistent with the expected direction with running faster than 5 cm/s, (3) the duration of the corresponding theta cycles ranged from 100 to 200 ms, (4) more than 3 place cells and emitted at least 5 spikes during the corresponding theta cycle.

Averaged sequence was constructed over the laps by aligning the probabilities over the theta cycle mid-time point and position. Each and averaged sequence were quantified within an area that decoded probabilities ±35 cm around the position of animal, and ±1⁄4 theta cycle around the mid-time point of the theta sequence.

Probability difference

The area used to quantify were divided equally into four quadrants. The Quadrant 2 represents the region both physically and temporally ahead of the rat, whereas as the Quadrant 3 represents the region physically and temporally behind the rat. Quadrant 1 and 4 represents the region that opposite with the animal’s current running direction. The summed decoded probabilities in Quadrant 1 and 4 were subtracted from the summed decoded probabilities in Quadrant 2 and 3 then divided by the summed decoded probabilities in all four quadrants. This value represents the difference between the decoding probability of the normalized current path and the opposite path in the theta sequence. Positive differences would imply theta sequences sweeping ahead in current moving direction of the rat, whereas differences close to zero would indicate a lack of sequential structure.

Weighted correlation

Weighted correlation of each or averaged theta sequence was calculated from the same area mentioned above. With the decoding probability (prob) as the weight of the position estimation, the correlation coefficient between the time (t) and the decoding position (p) is as follows:

where weighted covariance between time and decoded position is as follows:

and weighted means of time and decoded position are as follows:

and

The apparent sweep-ahead structure in the current running direction of animal produced large positive correlations (close to one), whereas the lack of sequence structure produced close to zero correlations. For further analysis, we use weighted correlation to quantify the theta sequence from different decoding result.

Flip correlation difference

We defined an Indicator to measure the front-sweep structure, named Flip correlation difference, which of each or averaged theta sequence was calculated also from the same area mentioned above. In the temporal dimension, the correlation coefficient between the decoding probability and the decoding probability after performing both temporal and spatial flips was calculated, from which the correlation coefficient between the decoding probability and the decoding probability after performing only temporal flips was subtracted. Positive difference represents intact sequential structure in the positive, whereas differences close to zero or negative would indicate a lack of sequential structure or defective sequential structure.

Gamma phase shift analysis of cell ensemble

The gamma cycle with maximal spiking across all simultaneously recorded cells was defined as cycle 0. Gamma cycles occurring before cycle 0 are represented by decreasing integer values, after by increasing integer values. Incomplete cycles at the beginning or end of the sequence or during immobile (speed<5cm/s) were excluded from analyses. We estimated the gamma phases and theta phases of spike times for place cells that exhibited spikes within at least two slow gamma cycles or at least three fast gamma cycles within a single theta cycle. The number of cycles analyzed within each theta sequence was limited to three (cycles -1 to 1) for slow gamma and five (cycles -2 to 2) for fast gamma (data not shown). Two-dimensional histograms of gamma phases and theta phases for spikes from each cycle number were binned into 40 bins and smoothed across nine bins.

Slow gamma phase shift analysis of individual cell

We employed Morlet’s wavelet method to calculate the time-varying power of slow gamma. The slow gamma energy in each theta cycle was estimated, considering only cycles with energy higher than 1.5 times the standard deviation as slow gamma events. The global 0 phase of slow gamma, set at the trough, determined the phase shift. For each cell, spikes in different slow gamma cycles within each theta cycle were identified. Spike pairs, representing adjacent spikes from different slow gamma cycles, were used to calculate the average phase difference as the slow gamma phase difference. In cases of multiple spikes in the same slow gamma cycle, their average phase was taken. Negative values indicated phase precession, while positive values indicated procession. To maintain firing location constancy, spikes used for phase difference calculation were assigned a random slow gamma phase within their cycle, repeated 1000 times for a distribution. In the analysis, we calculated one-sided, corrected empirical p-values [(r+1)/(n+1)] by comparing real data to its corresponding shuffle data. Where r is the number values from the shuffled distribution that are either smaller or larger (depending on the hypothesis) than the observed value and n=1000 is the number of shuffling.

Correlation between the slow gamma phase precession and weighted correlation

The theta cycle satisfying the criteria (see Slow gamma phase shift analysis) for calculating the phase difference was chosen to estimate the relationship between slow gamma phase precession and weighted correlation. On the one hand, the Pearson correlation coefficient was used to estimate the linear correlation between them. On the other hand, the correlation between slow gamma phase precession and theta sequence structure was analyzed by two-dimensional histogram. We binned slow gamma phase difference and weighted correlation into ten bins to construct a two-dimensional histogram, and smoothed it across three bins for visualize. The Two-dimensional histogram was divided into four quadrants and the sum of probability in Quadrant 1 (Q1) was calculated. The statistical significance was determined by calculating the one-sided, corrected empirical p-value. The shuffle data was same as slow gamma phase shuffle. The true data was considered to be significantly concentrated in Q1. The null hypothesis was the probability in Q1 of shuffled data was larger than the true data. If the p value was less than 0.05, the null hypothesis was rejected.

Histology

At the end of the experiment, histological sections were prepared as follows to verify the tetrodes position. Rats were intraperitoneally injected with a lethal dose of urethane solution. The heart was subsequently perfused with phosphate buffer and formalin to fix the brain tissue. Rat brains were then cut into 30-μm coronal sections and Nissl Staining was performed to determine the location of tetrode in CA1.

Quantification and statistical analysis

Data analyses were performed using SPSS (IBM) and custom Matlab (the Math Work) and R scripts. Paired t-test was performed to test between exFG-sequences and exNFG-sequences during early and late development (Fig.3H). Student’s t test was performed to test the spike number, mean firing rate in field, place fields size and spatial information between the FG-cell and NFG-cell (Fig.5D-G). One-sample t test was conducted to assess whether the slow gamma phase shift was significantly lower than 0 (Fig.7B). Repeated measures ANOVA was performed to test the weighted correlation of theta sequences (Fig.2D, Fig.3C,E,G) and running speed (Fig.3F) as a function of laps. Kolmogorov-Smirnov test was performed to test the distribution of place field COM (Fig.5C) on the track using the Matlab function ‘kstest2’. Multi-sample Mardia-Watson-Wheeler test was performed to test the slow gamma phase shift of place cells across successive slow gamma cycles (Fig.6) using the R ‘circular-package’. Pearson’s correlation was performed to test the correlation between weighted correlation difference and the proportion of excluded cells (Fig.4C-D).

Acknowledgements

We thank Shuang Meng for help with rat behavioral pre-training; and all Zheng Lab members for helpful discussions. This work was supported by the grants from National Science and Technology Innovation 2030 Major Project of China (2022ZD0205000), National Natural Science Foundation of China (T2322021, 82271218 to C.Z., 12271272 to L.W., 81925020 to D.M., 82371886 to J.Y., 82202797 to L.W.), The Space Brain Project from Lingang Laboratory (LG-TKN-202204-01 to J.Y.), and the China Post-Doctoral Science Foundation (2022M712365 to L.W.).

Author contributions

Conceptualization, C.Z. and L.W.; Methodology and Software, N.W., L.W. and C.Z.; Investigation, N.W., Y.W., M.G., L.W. N.Z., X.W. and C.Z.; Formal Analysis, N.W. and L.W.; Writing – Original Draft, N.W. and C.Z.; Funding Acquisition, C.Z., M.D., L.W., J.Y., and L.W.; Supervision, C.Z. and D.M.

Competing interests

The authors declare that they have no competing interests.