Subjects
The data were obtained using 19 Long Evans rats (18 males and 1 female; 300–500 g at the time of implantation). Data from five of the animals have been used for other purposes in published data5,68. The rats were group-housed with three to eight of their littermates before surgery and were thereafter housed singly under enriched circumstances in large two-storey metal cages (95 × 63 × 61 cm) or in smaller Plexiglas cages (45 × 44 × 30 cm). They were kept on a 12 h:12 h light:dark schedule in humidity- and temperature-controlled rooms. Experiments were approved by the Norwegian Food Safety Authority (FOTS ID 18011 and 29893) and done in accordance with the Norwegian Animal Welfare Act and the European Convention for the Protection of Vertebrate Animals used for Experimental and Other Scientific Purposes.
Surgery and electrode implantation
In 19 rats we implanted Neuropixels silicon probes targeting either the MEC–parasubiculum region (ten rats, of which three were implanted bilaterally), the hippocampus (two rats) or both regions (seven rats). Neuropixels prototype phase 3 A single-shank probes26 were used in eight of the rats, and prototype 2.0 multi-shank probes27 were used in the other 11 rats. Probes targeting the MEC–parasubiculum were implanted 4.2–4.7 mm lateral to the midline and 0.0–0.3 mm anterior to the transverse sinus, at an angle of 18-25° in the sagittal plane, with the tip of the probe pointing in the anterior direction. Probes were lowered to a depth of 4,100–7,200 μm. Hippocampal probes were positioned vertically mediolaterally 1.4–3.0 mm from the midline and anteroposteriorly 1.9–4.0 mm posterior to bregma. In rats with probes in both the MEC–parasubiculum and hippocampus, the two probes were implanted in different hemispheres (except for one animal). Two of the rats with probes in the hippocampus had also received injections of an adeno-associated virus (AAV8) carrying the fluorescent marker mCherry and the HM4Di DREADD receptor bilaterally in the hippocampus 3–4 weeks before probe implantation. The recordings included in the present study were obtained before any administration of the DREADD agonist Descholoroclozapine. The implants were secured with dental cement. A jeweller’s screw in the skull above the cerebellum was connected to the probe ground and external reference pads with an insulated silver wire. The detailed procedure for chronic Neuropixels surgeries has been described elsewhere27. Postoperative analgesia (meloxicam and buprenorphine) was administered during the surgical recovery period. Rats were left to recover until they resumed normal foraging behaviour, at least 3 h after surgery.
Electrophysiological recordings
Instruments and procedures were similar to those described for Neuropixels recordings used in the lab5,26,27,68. In brief, neural signals were amplified (gains of 500 for phase 3A and 80 for 2.0 probes), filtered (0.3–10 kHz for phase 3A and 0.005–10 kHz 2.0 probes) and digitized at 30 kHz by the probe’s on-board circuitry. Signals were multiplexed and transmitted to the recording system along a tether cable. SpikeGLX software (https://billkarsh.github.io/SpikeGLX/) was used to control acquisition and configure the probes. A motion-capture system, based on retroreflective markers on the implant, OptiTrack Flex 13 cameras and Motive recording software, was used to track head position and orientation in 3D. The 3D tracking coordinates were subsequently projected onto the horizontal plane for estimation of 2D position and head-direction azimuth. Back markers were attached in one animal. Another camera (Basler acA2040-90umNIR) was used to capture overhead infrared video in a subset of the recordings. Overhead video frames were aligned to OptiTrack tracking data with an affine transformation between corresponding points in the video and tracking data. Timestamps from each data stream were synchronized as previously described5,68, by generating randomized sequences of digital pulses with an Arduino microcontroller and sending them to the Neuropixels acquisition system as direct TTL input and to the OptiTrack system and video camera by means of infrared LEDs placed on the edge of the arena.
Behavioural procedures
Recordings were obtained while the rats foraged in an open field, while they navigated for rewards on a linear track or on a wagon-wheel circular track, or during sleep. Owing to the previously reported gradual decay in signal quality over the first 7–14 days after probe implantation70, most recordings were made in the first week after surgery (full range, 0–151 days postoperatively). All behavioural tasks for a particular animal were performed in the same recording room (except for one recording in a novel room; Extended Data Fig. 9k,l), often consecutively on the same day. Recording sessions were sometimes interrupted to remove twists from the Neuropixels tether cable. During presurgical training, some of the rats were food restricted, maintaining their weight at a minimum of 90% of their free-feeding body weight. Food restriction was not used in any of the animals at the time of recording.
The open-field foraging task
Eighteen of the rats foraged for randomly scattered food crumbs (corn puffs or vanilla cream cookies) in a square open-field box with a floor size of 150 × 150 cm and a height of 50 cm. The floor was made of black rubber, and the walls were made of black expanded PVC plastic. The arena was placed on the floor centrally in a large room (16 or 21 m2) with full visual access to background cues. A large white cue card was affixed to one of the walls (at the same height as the wall, a width of 41 cm, and a horizontal placement at the middle of the wall). In all illuminated trials, at the time of the surgery, each rat was already highly familiar with the environment and the task (having experienced 10–20 training sessions before surgery, each lasting at least 20 min). Recording sessions lasted 23–141 min.
In one exceptional case, a rat foraged during recording in a dark room encountered for the first time (Extended Data Fig. 9k,l). In this experiment, a circular arena 150 cm in diameter was used, as described in a previous study68. The arena was encircled by thick, dark-blue curtains. All light sources in the recording room were turned off or occluded before the recording started.
In another exceptional case, a rat foraged in a 1 m x 1 m square environment with a transparent plexiglass floor. As well as the overhead camera, a Basler camera was placed under the floor to record the rat’s footsteps. DeepLabCut69 was used to track the positions of each paw, snout and tail base. Body-part coordinates were interpolated in time to match the 10-ms bins of marker-based tracking and neural data, and were spatially aligned with the marker-based tracking data using an affine transform between snout coordinates from the video and head coordinates from the marker-based tracking.
The linear track task
Five rats with MEC–parasubiculum implants, of which two also had hippocampal implants, shuttled back and forth on a 200 cm linear track with liquid rewards delivered at each end (chocolate-flavoured oat milk dispensed by tube). When the rat consumed a reward at one end of the track, the reward port at the opposite end was refilled. Before surgery, rats were trained on the track task until they consistently completed around 40 laps in one training session. Recording sessions lasted 45–66 min, with 18–48 min of running between reward sites.
The wagon-wheel task
Two rats with MEC–parasubiculum implants were tested on a wagon-wheel track, an elevated circular track 10 cm wide with two perpendicular cross-linking arms spanning the diameter of the circle5. The track was fitted with eight reward wells, placed halfway between each of the five junctions. Each of the wells could be filled with chocolate oat milk by attached tubing. At any given time, a pseudorandom subset of 1–4 of the wells was filled to encourage steady exploration of the entire maze. Before surgery, the rats were trained to asymptotic performance levels, for which they obtained at least 30 rewards per 30-min session.
The m-maze task
One rat with probes implanted in the hippocampus and the MEC–parasubiculum was tested on an m-shaped maze, consisting of a central arm connected to two side arms through a T junction (arm lengths of 1.2 m, a T length 1.1 m and a track width of 10 cm)10,71. Reward wells were placed at the ends of each arm and were remotely filled with a liquid reward. The rat was trained to visit each of the side arms alternately, returning to the central arm between each visit to the side arms (centre–right–centre–left–centre–right and so on), receiving rewards at the end of each successful trajectory. Before being exposed to the m-maze, the rat was trained to shuttle back and forth on a linear track daily for around a week with multiple sessions each day. After 2–3 days with repeated exposure to the m-maze, the rat displayed almost perfect task performance and was implanted. The data shown in Extended Data Fig. 8 were obtained during the night after surgery.
Natural sleep
Sleep recordings, including both REM and SWS epochs, were obtained from nine rats with MEC–parasubiculum implants (two of them had combined MEC–hippocampal implants). Sleep was promoted by putting the rat in a black acrylic box (40 cm × 40 cm floor, 80 cm high), lined with towel on the floor, before recording. The box walls were transparent to infrared, allowing the rat’s position and orientation to be tracked through the walls. Water was available ad libitum. During recording, the room lights were on and pink noise was played through the computer speakers to mask any background sounds. Sleep sessions typically lasted 2–3 h (see Extended Data Fig. 11a for an example recording).
Spike sorting and single-unit selection
Spike sorting was done using KiloSort 2.5 (ref. 27), with customizations as previously described5. To exclude low-firing units, fast-firing interneurons and contaminated clusters, units were excluded if they had a mean spike rate of less than 0.1 Hz or greater than 10 Hz (less than 0.025 Hz or more than 5 Hz for hippocampal units), or if their waveforms had a large spatial footprint (a similar waveform amplitude across a wide range of channels), because this seemed to be a reliable indicator of poor cluster quality27. The waveform footprint was expressed as the anatomical spread of recording channels in which the unit was detected (with detection defined as at least 10% of the maximal amplitude of the unit across all channels), weighted by its waveform amplitude on each channel27:
$$\mathop{\sum }\limits_{i=1}^{N}{{\bf{w}}}_{i}({({{\bf{x}}}_{i}-\bar{{\bf{x}}})}^{2}+{({{\bf{y}}}_{i}-\bar{{\bf{y}}})}^{2}),$$
where xi, yi refers to the x, y location of the ith recording channel, \(\bar{{\bf{x}}}\) and \(\bar{{\bf{y}}}\) are the centres of mass of the recording channel positions where the unit was detected, \({{\bf{w}}}_{i}\) is the amplitude of the waveform at the ith recording channel, and N is the number of channels where the unit was detected. Units were excluded if the spatial footprint of their waveforms exceeded 35 μm (or 50 μm for hippocampal units). Units that were recorded on sites located outside the regions of interest (the MEC–parasubiculum or hippocampus) were excluded from further analysis.
Preprocessing and temporal binning
During awake sessions, only time epochs in which the rat was moving at a speed greater than 5 cm s−1 were used for spatial analyses. Spike times were binned in 10-ms time bins for all population analyses (unless otherwise specified), and tracking data were resampled at the same time intervals to align them with the spike-count data. For computational reasons, awake sessions were truncated in length to the nearest multiple of 100 s (given by the chunk size for the LMT model) by trimming the last part of the behaviour session.
Rate maps and angular tuning curves
To generate 2D rate maps for the open-field arena and wagon-wheel track, position estimates were binned into a square grid of 2.5 × 2.5 cm bins. For each bin, we calculated the firing rate of each cell (the number of spikes in the bin divided by time spent in the bin). Rate maps were smoothed with a cross-validated smoothing procedure. In brief, the recording was split into ten folds of equal duration, and the firing rate y during each fold was compared with the expected firing rate \(\widehat{y}\) based on the rate map calculated over the remaining nine folds and smoothed with a gaussian kernel of width σ. The value of σ (1 cm < σ < 50 cm) that minimized the mean squared error of the firing-rate prediction (using the MATLAB function fminbnd) was chosen to smooth the rate map. The same procedure was used to compute spatial rate maps with respect to latent position signals from the LMT model (see the ‘The LMT model’ section). For PV decoding analyses, a fixed-width gaussian kernel was used to smooth the rate maps (σ = 7.5 cm). Spatial autocorrelations and grid scores were calculated as described previously37, based on the individual cells’ rate maps.
Angular tuning curves with respect to head direction, theta phase or internal direction (see the ‘Decoding of internal direction based on PV correlations’ section) were calculated by binning the angular variable into 60 evenly spaced angular bins. For each 6° bin, the spike rate was calculated as the number of spikes divided by time spent in the bin. Angular tuning curves were smoothed with the same cross-validated smoothing procedure as the spatial rate maps (0.01 rad < σ < 1 rad), except in PV-decoding analyses, for which a fixed-width gaussian kernel of σ = 12° was used.
Identification of grid cells and grid modules
Grid cells were detected as groups of cells corresponding to grid modules by finding clusters of co-recorded cells that expressed similar spatially periodic activity in the open field, based on a similar procedure described previously5. In brief, 2D autocorrelograms were calculated from the coarse-grained spatial rate maps of each cell (10 cm × 10 cm bins, no smoothing across bins). Autocorrelogram bins within a central radius of two bins or beyond an outer radius corresponding to the rate-map size were masked, before vectorizing and concatenating the autocorrelograms in a matrix. Considering the spatial autocorrelograms of all cells as a point cloud, for which each point (autocorrelogram) has an N-dimensional position representing the value of each of its N spatial bins, the Manhattan distances between all points were calculated, and each point’s 30 nearest neighbours were identified. The resulting neighbourhood graph was given as input to the Leiden clustering algorithm, which was used to partition the spatial autocorrelograms into clusters, using a resolution parameter of 1.0 (1.5 for sessions with more than 1,000 units). Clusters that contained cells with clear and consistent grid patterns were classified as candidate modules of grid cells. For each cluster, grid periodicity was measured by the grid score of the median autocorrelogram across all cells in the cluster. Grid-pattern consistency was measured by computing the Pearson correlation between the average autocorrelogram of the cluster and the autocorrelogram of each individual cell. The median consistency across all cells in the cluster was defined as the grid consistency of the cluster. For a cluster to be classified as a grid module, three criteria had to be fulfilled: a cluster grid score greater than 0.3; a grid pattern consistency greater than 0.5; and the cluster needed to contain a minimum of ten cells. In some recording sessions, single grid modules seemed to be split into two clusters with similar spacing and orientation. For this reason, we added a step that merged grid clusters if the correlation between their average autocorrelograms was greater than 0.7.
Subclasses of grid cells defined by differential bursting
To measure the tendency of cells to fire in bursts, we devised a burst score (BS) based on the firing-rate autocorrelograms of each cell (time range ±50 ms, bin width 1 ms). The autocorrelogram value of the centre bin was set to zero, and the autocorrelogram counts were normalized by the mean. Next, we compared the autocorrelogram values at short time lags (2–10 ms) with those found at longer time lags (13–50 ms):
$${\rm{B}}{\rm{S}}=\frac{1}{b-a}\mathop{\sum }\limits_{i=a}^{b}{{\bf{y}}}_{i}-\frac{1}{d-c}\mathop{\sum }\limits_{i=c}^{d}{{\bf{y}}}_{i},$$
where yi is the value of the ith bin of the mean-normalized autocorrelogram, and a, b, c and d are the indices of autocorrelogram bins corresponding to time lags of 2, 10, 13 and 50 ms. Thresholds for classification of bursty versus non-bursty cells were determined by inspection of the distribution of burst scores across the sample (bimodal with a trough around −0.2; Extended Data Fig. 5a). Cells with a burst score greater than 0 were classified as bursty, cells with burst scores less than −0.4 were classified as non-bursty, and cells with burst scores between −0.4 and 0 were left unclassified.
Grid cells could be classified into bursty and non-bursty subclasses on the basis of their burst scores, but a third subclass of grid cells was revealed by closer inspection of autocorrelogram shape or theta phase modulation (Extended Data Fig. 5a). To reliably identify these three subclasses in an unsupervised manner, we applied an unsupervised clustering algorithm to the temporal autocorrelograms of all grid cells5. Temporal autocorrelograms were computed for each cell by calculating a histogram of the temporal lags between every spike and all surrounding spikes within a ±100 ms window, using 2-ms bins. The histogram was then divided by the mean value and concatenated in a matrix, discarding any cells with less than 100 counts in the autocorrelogram. Principal components analysis (PCA) was applied to the matrix of autocorrelograms, treating individual cells as observations and keeping the first ten components for further analysis. Next, a neighbourhood graph was constructed by computing the Manhattan distance between all pairs of points in the 10D point cloud and finding each point’s 150 nearest neighbours. The graph was used as input to the Leiden clustering algorithm (resolution parameter 0.2). The clustering algorithm detected three distinct clusters with unique autocorrelogram shapes, in agreement with the three subclasses of grid cells identified on inspection of the theta phase modulation and burst scores across all grid cells (Extended Data Fig. 5a). Most cells in the first two clusters had positive burst scores (referred to as bursty type I and bursty type II), but cells in the third cluster had negative burst scores (referred to as non-bursty).
Classification of direction-tuned cells
Cells were classified as tuned to HD if their HD tuning curves differed significantly from a uniform distribution (P < 0.001, Rayleigh test for non-uniformity) and was stable across the first and second half of the recording session (P < 0.01, Pearson correlation between tuning curves from the first and second half). No further criteria on tuning strength were applied to classify HD tuning, because many internal direction-tuned cells display weak HD tuning and there was no clear cut-off between HD-tuned cells and non-HD-tuned cells (Extended Data Fig. 6a). Cells were classified as internal direction cells if their tuning curves, with respect to the decoded internal direction (see the ‘Decoding of internal direction based on PV correlations’ section), passed the above criteria for non-uniformity and stability, and had a mean vector length greater than 0.3. Tuning width was defined as two standard deviations of the tuning curves.
Theta phase estimation
Theta phase was extracted from the population spiking activity of all units (including fast-firing putative interneurons) in the MEC–parasubiculum region. Spike times were binned into 10-ms bins and the resulting spike counts were bandpass-filtered with a second-order Butterworth filter in the theta range (5–10 Hz). PCA was applied to the matrix of bandpass-filtered spike counts (with units as variables and time bins as observations). The first two principal components typically contained an oscillating circular representation corresponding to the theta rhythm. Theta phase at time t was defined as the direction of the projection of the PV at time t onto the plane defined by PC1 and PC2. The phase with minimal firing activity was defined as zero. In the six rats with dual entorhinal–hippocampal implants, hippocampal sweeps were referenced to the theta phase extracted from MEC–parasubiculum activity. In the two rats with probes only in the hippocampus, the theta phase was estimated from the activity of the hippocampal population.
Theta-cycle skipping
Theta skipping (the tendency for cells to fire only on alternate theta cycles) was quantified by computing a theta-skipping index (TSI), as in previous descriptions of the phenomenon10,41,42. In brief, firing-rate autocorrelograms were generated for each cell (time range ±500 ms, bin width 5 ms, gaussian smoothing σ = 10 ms), and the relative height of the second theta peak, p2, compared with the first theta peak, p1, was determined as
$${\rm{TSI}}=\frac{{p}_{2}-{p}_{1}}{\max ({p}_{1},{p}_{2})},$$
where p1 is defined as the maximum autocorrelogram value between lags 90 and 170 ms, and p2 is the maximum value between lags 180 and 300 ms. Cells were classified as theta skipping if the TSI was positive, that is, if the second theta peak was higher than the first theta peak. The TSI has a clear interpretation only for cells that are theta-rhythmic and is therefore reported for cells that were classified as modulated by theta phase. Cells were classified as theta-phase modulated if the theta-phase tuning curves differed significantly from a uniform distribution (P < 0.001, Rayleigh test for non-uniformity) and were stable across the first and second halves of the recording session (P < 0.01, Pearson correlation between tuning curves from the first and second halves). Cells were classified as non-rhythmic if they did not meet the above criteria and had a theta-phase mean vector length of less than 0.2.
Decoding sweeps according to PV correlations
This section describes how we visualized and quantified sweeps. We first decoded the position from the activity of the entire set of MEC–parasubiculum neurons. For each temporal bin, we correlated the instantaneous population activity (the PV) with the session-averaged population activity for each location in the environment (a reference PV, rPV). The inputs to the PV correlation decoder were as follows: first an N × T matrix of temporally smoothed firing rates (gaussian kernel σ = 10 ms), with T columns corresponding to time bins and N rows corresponding to neurons; and second, an N × M matrix of spatial tuning curves, with M columns corresponding to position bins and N rows corresponding to neurons. The spatial tuning curves were normalized by dividing the tuning curve of each neuron by its mean value. To decode position across successive time bins, we computed the Pearson correlation coefficient between each PV (columns of the spike- count matrix) and each rPV (columns of the tuning-curve matrix). This yielded a vector of correlation values for each time step, with elements corresponding to each [x, y] location in the environment. The decoded position was taken as the position bin with the highest correlation value, and the resulting decoded trajectory was smoothed with a σ = 8 ms gaussian kernel. To filter out unreliable estimates, the peak correlation value for each time bin was compared with a shuffled distribution of PV correlation values (computed by shuffling rows of the tuning-curve matrix). Decoded estimates were discarded if their PV correlation did not exceed the 99th percentile of the shuffled distribution. Decoding estimates were also discarded for time points for which fewer than five cells were active.
In each theta cycle, the decoded position swept outwards from a location slightly behind the tracked head position of the animal (Fig. 1). Thus, the starting locations for a series of sweeps formed a slowly evolving trajectory that roughly followed the animal’s running trajectory. This trajectory, referred to as the lowpass-filtered decoded trajectory and estimated by decoding position from spikes emitted in the beginning of each theta cycle, was used as a reference signal for measuring sweeps. Spike counts from the first half of each theta cycle were smoothed with a wide gaussian kernel (σ = 1.7 theta cycles) before the PV-correlation decoding method was used to decode position across all bins. The resulting decoded trajectory was smoothed with a σ = 10 ms gaussian kernel. Sweeps seemed to be more reliably anchored to the lowpass-filtered decoded trajectory than to the tracked head position of the animal. This difference was particularly evident when rats navigated in darkness, when sweeps and lowpass-filtered signals could deviate substantially from the rat’s actual trajectory (Extended Data Fig. 9k,l).
Extracting sweep trajectories
Individual sweeps, defined as smooth spatial trajectories in each theta cycle, were extracted from the decoded position trajectory by a simple sequence-detection algorithm. In each theta cycle, candidate sweeps were identified as the longest sequence (the highest number) of consecutive valid time bins for which the decoded position jumped less than 20 cm and changed direction less than 90° between consecutive 10 ms bins. Candidate sweeps were truncated to maximize the net Euclidian distance from beginning to end, which effectively removed folds at either end of the sweep. A sweep vector, s, was defined as the vector between the lowpass-filtered decoded position at the beginning of the theta cycle and the most distal point of the candidate sweep. To identify sweep trajectories that were fairly straight and well described by the sweep vector, we measured, for each candidate sweep, the goodness of fit, r2, between the sweep-vector axis and the collection of [x, y] points in the sweep, represented by the vectors x and y:
$${r}^{2}=1-\frac{{\rm{var}}({\bf{e}})}{{\rm{var}}({\bf{x}})+{\rm{var}}({\bf{y}})},$$
where e is a vector of residuals of all points with respect to the sweep-vector axis, such that ei is the residual of point [xi, yi]. Candidate sweeps were kept for further analyses if they included at least four samples and had r2 > 0.5 (meaning that most of the variance is explained by the sweep axis). Sweep direction and sweep length were defined as the direction and magnitude of the sweep vector. Sweep prevalence was computed as the fraction of theta cycles (where the animal moved faster than 15 cm s−1) containing sweeps that passed the selection criteria. When comparing sweep prevalence in different populations of cells (grid cells versus non-grid cells), the number of cells in the largest population was down-sampled (repeated 100 times) to match the number of cells in the other population.
For several analyses and visualizations (for example, Fig. 1b), sweeps (in allocentric coordinates) were transformed to head-centred coordinates. This was done by first subtracting the tracked position (for example, Fig. 1b) or lowpass-filtered decoded trajectory (for example, Extended Data Fig. 3e) and then negatively rotating each [x, y] coordinate by the animal’s head direction. Session-averaged sweeps (for example, Fig. 1b) were computed by first interpolating each head-centred sweep trajectory at 50 time points linearly spaced from the beginning to the end of the trajectory. Interpolated sweeps were grouped into those that followed a right sweep or left sweep, before an average sweep was computed for each group by taking the median position at corresponding time points within the sweep (Figs. 1b, 4e and 5d and Extended Data Fig. 3e,f).
Temporal delay between hippocampal and entorhinal sweeps
To quantify the temporal delay between hippocampal and entorhinal sweeps, we first decoded the position independently from activity in the hippocampus and the MEC–parasubiculum. Next, we computed, for each region, the offset between decoded position and the rat’s position, projected onto the rat’s head axis. Finally, we computed the temporal cross-correlation between head-centred decoded position from each region. The peak of the cross-correlogram nearest zero was taken as the temporal delay between hippocampal and entorhinal sweeps. MEC theta troughs were used as a reference for plots in Fig. 1i and Extended Data Fig. 3g, but only for visualization purposes. Similar delays were observed when the number of cells was balanced across regions by down-sampling cells in the region with more cells to match the number in the other region (data not shown).
Decoding of internal direction on the basis of PV correlations
Internal direction was decoded from the activity of all MEC–parasubiculum neurons using the same PV-correlation procedure described for position decoding except using angular tuning curves for head direction instead of spatial-position rate maps. Because we observed that internal direction cells were activated in discrete pulses during each theta cycle (Fig. 2c), the time bin corresponding to the theta phase with maximal activity was used to express internal direction in each theta cycle. The decoded internal direction α at time t was taken as the circular mean of all possible decoding angles weighted by the correlation values for each directional bin at time t:
$$\alpha =\text{angle}(\mathop{\sum }\limits_{j=1}^{M}{{\bf{r}}}_{j}\exp (i\times {{\boldsymbol{\theta }}}_{j})),$$
where θj denotes the angular value of the jth bin, i is the imaginary unit, rj is the correlation value for the jth angular bin and M is the total number of angular bins. The presence of negative correlation values did not affect the location of the mean. Similar results were obtained with Bayesian decoding for which posterior probabilities are always positive (Extended Data Fig. 4c,d).
For some analyses and visualizations (for example, Fig. 2d), we rotated the decoded internal direction (in allocentric coordinates) to a head-centred reference frame by subtracting the animal’s head direction.
Left–right alternation of sweeps and internal direction
The extent of directional alternation across successive sweeps and internal direction signals was characterized in a head-centred reference frame (see above) during periods when the animals moved faster than 15 cm s−1. The data were generally thresholded at 15 cm s−1 because alternation was more reliable during running (see Fig. 6f,g, which also includes data for speed thresholds down to 5 cm s−1). The prevalence of directional alternation was computed by counting triplets of theta cycles for which sweep direction or internal direction alternated in a left–right–left or right–left–right pattern (detected as sign inversions in the angles between successive directions) divided by the total numbers of theta-cycle triplets for which sweeps or internal direction were detected. The fraction of theta-cycle triplets with directional alternation was compared with a shuffled distribution of scores for which head-centred directions were randomly shuffled (1,000 iterations). In Extended Data Fig. 3i, we repeated the same analysis except with windows of alternation extending from 3 to 20 theta cycles. Directional alternation was visualized in temporal autocorrelograms of angles between successive head-centred sweep direction or internal direction. Autocorrelograms were computed as the circular correlation between the original trace of head-centred directions and a series of lagged versions of the signal (lags from −7 to 7 theta cycles). To find and visualize directional modes of sweep and internal direction angles, we computed histograms of head-centred directions that were conditioned on the decoded direction in the previous theta cycle. First, decoded directions were classified as left- or right-directed according to the sign of their angular offset from the previous cycle. Next, two histograms were computed, one for decoded directions that followed a left-directed angle and one for decoded angles that followed a right-directed angle. Because sweeps and internal direction angles alternate consistently, this procedure resulted in two unimodal distributions, one on either side of the reference head direction. The directional modes of sweep or internal direction were taken as the peak of each of the two conditional distributions.
In the artificial-agent simulation, a time-resolved measure of directional alternation was used to quantify instantaneous alternation of the agent’s chosen sweep direction. A three-sweep sliding window was used, such that at the ith timestep, a triplet of sweep directions αi−1:i+1 was selected. In the three-sweep sliding window, the two angles between consecutive sweep pairs were calculated: a = αi − αi−1 and b = αi+1 − αi. The alternation score s was therefore computed as:
$$s=\frac{| a-b| }{2\,\max (| a| ,| b| )},$$
where || denotes the absolute value. The alternation score ranged from 0 to 1, where 1 indicates perfect alternation.
Single-module decoding
To decode position from individual grid modules, the PV-correlation decoding analysis (first applied on all cells) was next applied to subsets of cells that belonged to individual grid modules. The decoded position from each grid module was mapped onto the hexagonal unit tile of the grid module using a cross-correlation procedure. First, we computed a template grid pattern for each grid module by averaging the spatial rate-map autocorrelograms across all cells in the module. The template grid pattern was cross-correlated with the 2D distribution of PV-correlation values at each time step. Peaks in the resulting spatial cross-correlogram were detected and, for each cross-correlogram, the peak nearest to the origin was taken as the decoded position. To ensure that the decoded position was within the bounds of the central unit tile, the decoded trajectory was wrapped around the three grid axes of the template grid pattern. Single-module sweeps were detected as described for whole-population decoding (by finding consecutive time bins within each theta cycle for which the decoded trajectory formed a smooth trajectory), except that spatial and directional offsets between consecutive decoded positions were computed with periodic boundary conditions derived from the template grid pattern. For some visualizations (for example, Figs. 1e and 4b), the single-module sweeps were aligned to the behavioural trajectory of the rat by subtracting the offset between the animal’s tracked position and the lowpass-filtered decoded trajectory at the beginning of each theta cycle.
Theta phase precession
To show and quantify phase precession in grid cells during 2D open-field foraging, we adapted approaches used in previous studies72,73. In brief, we first detected vertices of the grid pattern as local maxima in the spatial rate maps. Next, the vector from the rat’s position to the nearest grid vertex was computed for each time step. The resulting vectors were projected onto the decoded internal direction (which serves as a reliable measure of the sweep axis), yielding a scalar ‘in-field position’ at each time step. In-field position is negative when the animal is approaching a grid field and positive when it is heading out of the grid field. Phase precession was measured as the circular–linear correlation between theta phase and in-field position at the time of each spike, following ref. 74. Cells with a negative correlation coefficient and P < 0.01 were defined as phase precessing.
Phase precession and theta skipping in simulated neurons
To investigate whether alternating sweeps and internal direction signals were sufficient to explain phase precession and cycle skipping in individual cells, we generated spikes from firing-rate models of grid cells and internal direction cells.
Firing-rate predictions for grid cells were generated on the basis of: first, the spatial rate map of an idealized grid cell (implemented as a hexagonal tessellation of 2D gaussian bumps with a spacing of 75 cm, σ = 0.5 cm and an amplitude of 30 Hz); and second, a time-varying position signal (a recorded trajectory of a rat in the open field, with or without added sweep-like spatial offsets with length up to 35 cm). The firing-rate prediction at each time step was taken as the value of the rate map at the location of the time-varying signal divided by the sample rate (100 Hz). Poisson spike trains were generated from the predicted firing rate. Phase-precession plots were made as described in the previous section. Only spikes from the second half of the theta cycle (corresponding to the outwards phase of each sweep and the active phase of bursty grid cells) were used to compute firing-rate autocorrelograms, because any sweep-induced theta skipping would be visible only during this part of the theta cycle.
Firing-rate predictions for internal direction cells were generated on the basis of tuning curves for internal direction and theta phase, and a time-varying internal direction signal and a time-varying theta-phase signal. In Extended Data Fig. 6h, the tuning curves from a recorded cell and decoded internal direction were used to generate firing-rate predictions. In Extended Data Fig. 6i, simulated tuning curves from a recorded cell and simulated internal direction were used to generate firing-rate predictions. Tuning curves, internal direction and theta phase were simulated as Von Mises distributions with means of 60° and 180°, and concentrations of 6 and 1.5, respectively. Internal direction was simulated by adding ±30° offsets to the animal’s head direction in an alternating or non-alternating pattern. The firing-rate prediction, λ, was estimated by combining the contributions of theta phase and internal direction tuning following a linear–nonlinear Poisson model:
$${\boldsymbol{\lambda }}=\exp ({{\boldsymbol{\beta }}}_{0}+{{\boldsymbol{\beta }}}_{id}{{\bf{X}}}_{id}+{{\boldsymbol{\beta }}}_{\theta }{{\bf{X}}}_{\theta })/{\rm{d}}t,$$
where β0 is the log mean firing rate, βid and βθ are the log tuning curves for internal direction, and theta phase and data matrices Xid and Xθ represent the time-varying internal direction and theta phase signals. Poisson spike trains were drawn from the firing-rate prediction.
Identification of putative excitatory connections between cells
Putative excitatory monosynaptic connections between pairs of co-recorded neurons were identified by detecting short-latency, short-duration peaks in the firing-rate cross-correlograms (CCGs), following previous procedures44,45,75,76,77. CCGs were computed with 1 ms bins over a ±50 ms window. A baseline CCG was computed by convolving raw CCGs with a hollowed gaussian kernel (σ = 5 ms, hollow fraction 60%), which has been shown to approximate a jittered CCG75. Neuron pairs recorded on the same probe, for which each of the neurons emitted at least 2,000 spikes and the raw CCG contained at least 1,000 counts, were tested for functional connections. The baseline CCG was subtracted from the raw CCG and subsequent peak detection was performed on the baseline-corrected CCG. Neuron pairs were classified as putatively connected if the highest positive CCG peak satisfied the following five conditions: first, the peak occurred within an asymmetric time range consistent with monosynaptic excitation (0.7–4.7 ms), in agreement with experiments in which the discharge latencies of synaptically coupled cell pairs were measured more directly76,78,79; second, the peak height exceeded five standard deviations of the baseline-corrected CCG; third, the P-value of the peak was less than 0.001 (estimated from a Poisson distribution with continuity correction, as in ref. 75); fourth, the width of the peak (defined as the set of bins contiguous with the peak whose values exceeded half of the peak height or two standard deviations of baseline and had P < 0.01) was less than 3 ms (consistent with the precise spike timing expected from a monosynaptic connection); and fifth, the width of the peak did not overlap with the zero-lag bin (suggestive of common input). Candidate connections were also discarded if any of the non-peak bins exceeded 2.5 standard deviations of the baseline-corrected CCG or if any of the bins in the anticausal direction had P < 0.01.
Overall connection probability was computed by dividing the total number of putatively connected neurons by the total number of pairs that were checked for connections. Similarly, target-specific connectivity rates were computed by dividing the number of connections from one functional cell class to another by the number total number of pairs.
Decoding internal direction with PCA or UMAP
The internal direction signal was also decoded in an unsupervised manner, with PCA or uniform manifold approximation and projection80 (UMAP). In this approach, the high-dimensional neural activity was projected down to a 2D subspace to characterize the trajectory of population activity on a low-dimensional manifold. Only rhythmic direction-tuned cells (with head-direction mean vector length greater than 0.3 and theta-phase mean vector length greater than 0.3) were included in the analysis, to avoid interference between the ring-like manifold of interest and other ensemble representations, such as spatial signals from pure grid cells. Spike counts from n neurons were binned into t time bins corresponding to individual theta cycles. Theta-cycle time bins were used instead of 10-ms time bins to prevent global within-cycle firing rate fluctuations from driving the results. PCA was done on the resulting t-by-n matrix. Next, internal direction was read out either directly from the PCA output or by applying a second dimensionality-reduction step on the data using UMAP (Extended Data Fig. 4c,d).
In the PCA decoder, internal direction was read out by projecting the neural data onto the first two eigenvectors and taking the arctangent of each [x, y] coordinate in the resultant 2D projection. Because these angles were arbitrarily rotated with respect to the environment, they were aligned to the environment before further analysis. We assumed that internal direction and head direction had equal mean directions, and subtracted the average difference between the two signals from the decoded signal to align it with respect to the environment.
In the UMAP decoder, scores of the 20 principal components with the largest explained variance were used as input to the nonlinear dimensionality-reduction algorithm UMAP with these hyperparameters: n_components=3, metric=correlation, n_neighbors=199, min_dist=0.3, init=spectral. This yielded a 3D embedding of the high-dimensional population activity. The 3D UMAP point cloud typically showed a clear circular shape. The 3D point cloud was then collapsed to 2D by projecting the points onto the best-fit 2D plane. Hence, the 2D points were converted into angles by taking the arctangent of each x and y coordinate. The resulting angles were aligned to the environment by subtracting the average offset from head direction, as described above.
Bayesian decoding of position
Sweeps could also be decoded using Bayesian reconstruction (Extended Data Fig. 4a). Position was decoded at 10-ms time steps from the matrix of firing rates and tuning curves from N neurons, with an assumption of Poisson firing and a flat position prior81:
$$P({\bf{x}}|{\bf{y}})\propto \exp (\mathop{\sum }\limits_{i=1}^{N}{{\bf{y}}}_{i}\log (\,{f}_{i}({\bf{x}}))-{\rm{d}}t\mathop{\sum }\limits_{i=1}^{N}{f}_{i}({\bf{x}})),$$
where P(x|y) is the conditional probability for the rat’s 2D location x, given the observed spike count y and position tuning curves f(x). The decoded position was taken as the position bin that maximized P(x|y).
The LMT model
The PV-correlation method has two limitations. First, it can only be used to decode positions and directions that the animal has physically sampled, because it relies on tuning curves with respect to the animal’s tracked position and head direction. This is a problem for hippocampal data in particular, because the neural correlation structure varies between environments and brain states82,83,84, making it impossible to decode position in one environment based on reference tuning curves from a different context. A second challenge is that the sharp spatial and directional tuning of grid cells and internal direction cells is obscured in standard time-averaged reference tuning curves, because sweeps and internal direction deviate from the tracked position and head direction (Extended Data Figs. 5f and 6d). To simultaneously extract sweeps through unvisited space and characterize spatial tuning directly from neural activity, we therefore adapted the LMT model in ref. 46. In this framework, sweeps may be considered as hidden or ‘latent’ trajectories on a neural manifold that is not directly observable. The goal of the LMT model is to infer both the latent trajectory and each cell’s tuning to locations on the manifold, based on neural population activity. The model assumes that the latent variable evolves smoothly with time, that individual neurons are smoothly tuned to locations on the manifold, and that neurons fire according to a Poisson process (see Extended Data Fig. 10 for a schematic and ref. 46 for details). For a multidimensional (vector-valued) latent variable x(t), the temporal evolution of component j is modelled as a Gaussian process:
$${{\bf{x}}}_{j}(t) \sim {\mathcal{G}}{\mathcal{P}}(0,{k}_{t}),$$
where kt is a temporal covariance function, \(k(t,{t}^{{\prime} })\,\triangleq {\rm{cov}}\,({{\bf{x}}}_{j}(t),\,{{\bf{x}}}_{j}({t}^{{\prime} }))\). In this case, the exponential kernel k(t, t′) = r exp(−|t−t′|/l) is used, with variance r and length-scale l, respectively controlling the amplitude and smoothness of the latent variable. The log tuning curves f(x) are also modelled as Gaussian processes, with the log tuning of the ith neuron expressed as:
$${f}_{i}({\bf{x}}) \sim {\mathcal{G}}{\mathcal{P}}(0,{k}_{{\bf{x}}}),$$
where kx is a spatial covariance function, in this case a Gaussian kernel \({k}_{x}({\bf{x}},{{\bf{x}}}^{{\prime} })=\rho \,\exp (-{\Vert {\bf{x}}-{{\bf{x}}}^{{\prime} }\Vert }_{2}^{2}/2{\delta }^{2})\) with variance ρ and length scale δ. Furthermore, an L-1 penalty parameter, λ, was included to enforce the sparsity of the fitted tuning. The value of the latent variable x(t), in conjunction with each cell’s tuning curve, fi(x), predicts the cell’s log-firing rate, which is then transformed with an exponential nonlinearity into a Poisson-distributed spike count yi,t as:
$${y}_{i,t}|\,{f}_{i}\,,{{\bf{x}}}_{t} \sim {\rm{Poiss}}(\exp (\,{f}_{i}({{\bf{x}}}_{t})).$$
The predicted and observed spike trains were compared, yielding a log-likelihood value. The model was fitted using an expectation-maximization algorithm that maximizes the log-likelihood of the spike trains by separately optimizing the latent variable and tuning curves in alternation. Over multiple iterations of this two-step optimization procedure, both evolve to capture the latent dynamics in the neural population activity, thus improving the prediction of the observed spikes.
The original LMT framework uses a single latent variable to predict the neural activity; however, as a form of Poisson regression, the LMT model can be trivially combined with other Poisson regression models, as a means of extracting other factors of interest, or to regress out noise. In the present study, we formulated the activity of each neuron as a sum of log firing-rate contributions from five input variables:
$${\bf{Y}} \sim {\rm{Poiss}}(\exp ({{\bf{M}}}_{{\rm{id}}}+{{\bf{M}}}_{{\rm{pos}}}+{{\bf{M}}}_{{\rm{theta}}}+{{\bf{M}}}_{{\rm{hd}}}+{{\bf{M}}}_{{\rm{pop}}})),$$
where Y is the time-by-neurons matrix of predicted spike counts and Mvarname is a time-by-neurons matrix of log firing-rate contributions from model varname. The first two contributions correspond to the two latent variables of interest: internal direction (Mid) and position (Mpos). Internal direction was modelled as a 1D latent circular variable (hyperparameters: ρ = 0.1, r = 100, δ = 0.5, l = 0.1 and λ = 1) that was initialized with the animal’s tracked head direction. Position was modelled as a 2D latent variable (hyperparameters: ρ = 0.1, r = 10, δ = 6, l = 0.015 and λ = 1) that was initialized with the animal’s tracked position. Importantly, grid-like periodic mapping of space is not assumed by the model, which has no periodic constraints or boundary conditions.
The final three covariates (theta phase Mtheta, head direction Mhd and population firing rate Mpop) are known to modulate MEC–parasubiculum activity4,19,37,85, and are here included as ‘noise’ covariates to regress out their substantial contributions to the neural activity, hence reducing the likelihood of them influencing the extracted latent variables. Theta phase and head direction were modelled as circular 1D variables in the LMT framework; however, the variables were respectively fixed at the values of measured theta phase and head direction, and only tuning curves were optimized (hyperparameters: ρ = 1,000, δ = 2 and λ = 0 for theta phase, and ρ = 0.1, δ = 0.5 and λ = 1 for head direction). The population firing rate model was implemented as a generalized linear model Mpop = Xpopβpop, where Xpop is a time-by-neurons matrix of log population firing rates, and βpop is a row vector of learnt coefficients for all neurons. Population firing rate was computed as the average instantaneous firing rate across neurons, smoothed with a σ = 20 ms gaussian kernel.
Each step of fitting in the composite model consisted of serially updating the parameters for each of the five submodels by maximizing the log likelihood of the model. Because the latent trajectories may be arbitrarily rotated and distorted with respect to the physical environment, an alignment procedure was performed after model fitting46. The latent internal-direction trajectory was aligned to the tracked head direction of the animal by subtracting the average angle between the two signals. The latent position trajectory was aligned to the tracked position of the animal with an affine transformation. The fitted latent variables were used for most LMT-based analyses of sweeps and internal direction. For some visualizations (for example, that in Fig. 4d) and analyses during sleep, position and direction were decoded from neural activity and fitted LMT tuning curves using the Bayesian framework described in the previous section.
The GLM-based single-cell tuning model
Theta sweeps introduce offsets between the rat’s current location and the location represented by place and grid cells, resulting in smeared receptive fields when a cell’s spikes are plotted as a function of animal location (Extended Data Fig. 5e,f). To characterize the spatial tuning of individual cells independently, in a manner that accounts for sweeps, we used a GLM to model the spike train of a cell, y, as a function of a set of explanatory variables, X, parametrized by the learnt parameters β:
$${\bf{y}} \sim {\rm{P}}{\rm{o}}{\rm{i}}{\rm{s}}{\rm{s}}(\exp ({{\boldsymbol{\beta }}}_{{\rm{c}}{\rm{o}}{\rm{n}}{\rm{s}}{\rm{t}}}+{\boldsymbol{\beta }}{\bf{X}})).$$
Construction of the data matrix X
The model included four explanatory variables: position, head direction, internal direction (from LMT) and theta phase. The variables were expressed by using a basis-expansion procedure, by which a set of smooth basis functions was used to decompose each single variable into multiple variables. The weighting of each basis function was given by a corresponding β parameter, giving the GLM the flexibility to fit any smooth function of the input variable in question. The 2D position variable was expressed by a set of 2D Gaussian basis functions (σ = 2 cm) arranged in a 10-cm-spaced triangular grid that tiled the open-field arena and a surrounding buffer zone. The angular variables head direction, internal direction and theta phase were expressed by a set of 50 Von Mises functions (κ = 10) with equally spaced mean values from 0 to 2π. Basis expansion was done by evaluating each of the basis functions for a given value of the input variable in question.
The basis-expanded representations are high dimensional and multicollinear (that is, the basis function values are correlated). In a regression model, these attributes tend to cause overfitting, so we used PCA to produce a low-dimensional, orthogonal representation of the basis-expanded data matrix. PCA acted as a form of regularization. For position, the principal components explaining 99% of the variance were retained, reducing the dimensionality from 527 to 92. For head direction and theta phase, the principal components explaining 80% of the variance were retained, reducing the dimensionality from 50 to 6. After performing basis expansion and dimensionality reduction for the input variables, the resultant matrices for all input variables were concatenated into the GLM design matrix X:
$${\bf{X}}=[{{\bf{X}}}_{{\rm{p}}{\rm{o}}{\rm{s}}},{{\bf{X}}}_{{\rm{i}}{\rm{d}}},{{\bf{X}}}_{{\rm{h}}{\rm{d}}},{{\bf{X}}}_{{\rm{t}}{\rm{h}}{\rm{e}}{\rm{t}}{\rm{a}}}].$$
Theta-phase-dependent shifting of the position covariate
To model the effect of sweeps on position-modulated firing, we added a preprocessing step that applied a theta phase-dependent shift to the animal’s tracked position coordinates. Specifically, the 2D position coordinates (x) were parametrically shifted by a distance δ along the internal direction axis (α):
$${{\bf{x}}}_{t}^{{\prime} }={{\bf{x}}}_{t}+{{\delta }}_{t}[\cos ({{\alpha }}_{t}),\sin ({{\alpha }}_{t})],$$
where x′ denotes the shifted position coordinates. The shift quantity, δt, at each time point was modelled as a function of the current theta phase. Specifically, the shift quantity δ was modelled as a function shift parameter γ, fitted by the model, and the basis-expanded theta phase Xtheta:
$${\boldsymbol{\delta }}={{\boldsymbol{\gamma }}}_{{\rm{c}}{\rm{o}}{\rm{n}}{\rm{s}}{\rm{t}}}+{{\boldsymbol{\gamma }}{\bf{X}}}_{{\rm{t}}{\rm{h}}{\rm{e}}{\rm{t}}{\rm{a}}}.$$
The shifting parameters γconst and γ were fitted together with the GLM β parameters using a gradient-based solver with the finite-difference method (MATLAB function fminunc).
Across cells and recordings, the GLM shift model yielded sharper receptive fields than standard rate maps with respect to tracked position (Extended Data Fig. 9h). The model’s estimates of position tuning were more robust than those from the LMT model, because the latter depended on large numbers of co-recorded spatially modulated cells. Therefore, rate maps based on the GLM-shifted position were deemed most appropriate to use to identify grid cells and grid modules (see the ‘Identification of grid cells and grid modules’ section).
Sweeps through unvisited space
For analyses of sweeps and spatial tuning to never-visited locations outside the bounds of the wagon-wheel track, we first defined the area of space that the animal had visited. This was achieved by binning the 2D environment in 2.5-cm bins and finding all the bins that the rat had visited (resulting in a binary 2D map with values of 1 for visited bins and 0 otherwise). The bounds of the animal’s coverage was found by applying a binary dilation operation of the occupancy map (Matlab function imdilate with a disk-shaped structuring element, radius = 1), followed by a morphological closing operation (Matlab function imclose with a disk-shaped structuring element, radius = 1). The zero-valued bins in the resulting occupancy map were defined as never visited.
Head oscillations
For analyses of coupling between left–right alternation in the internal direction signal and body movements (Extended Data Fig. 2), we extracted vertical and lateral head oscillations from OptiTrack tracking data during linear track running. Because head oscillations are most prominent during fast, straight running86, analyses were restricted to epochs with running speed greater than 40 cm s−1. Vertical head speed was taken as the first derivative of the highpass-filtered vertical head position. Lateral head speed was computed by highpass-filtering 2D head position in the x–y plane, rotating the highpass-filtered signal to align with the animal’s movement direction and taking the first derivative in the lateral direction. Spectral analysis was done on 10-ms-binned head-speed signals and multi-unit activity using the multi-tapered Fourier transform, implemented by the Chronux toolbox (http://chronux.org/, function cohgramc). Non-overlapping 2 s windows were used, with a frequency band of 1 Hz and 3 tapers. All signals were bandpass-filtered between 1 and 20 Hz. The relationship between lateral head oscillations and internal direction signals was also investigated using phase-coupling analysis (Extended Data Fig. 2h). To this end, we extracted the phase of the lateral head-swing cycle from the Hilbert transform of the bandpass-filtered lateral head speed (pass band 2–10 Hz). Phases of 0° and 180° were set to match the times of right and left extremes in lateral head position, respectively. We next computed histograms of head-swing phase at times for which internal direction pointed to the left or right side of the animal’s head axis. The mean vector of these histograms was used to measure the consistency and strength of phase locking across animals.
Forelimb footsteps
For analyses of coupling between left–right alternation in the internal direction signal and forelimb footsteps (Extended Data Fig. 2i–k), we extracted forelimb footsteps from under-floor video. The position of paws, snout and tail base were extracted using DeepLabCut. The gait cycle, in which one cycle corresponds to the period from one right forelimb plant (phase 0) to the next, was extracted by projecting the vector between left and right forelimbs onto the head axis and computing the phase of the resulting signal using the Hilbert transform. Phase coupling between alternations in the internal direction signal and the gait cycle was examined using the same method as described for head oscillations.
Sleep-stage classification
Sleep stages were identified as described in previous studies5,19. First, we identified periods of sustained immobility (longer than 120 s, locomotion speed below 1 cm s−1, head angular speed below 6° s−1). These periods were subclassified into SWS and REM on the basis of delta- and theta-rhythmic population activity in the recorded cells. Population firing rate was computed by summing the binarized 10-ms spike counts from each cell. The rhythmicity of this aggregated firing rate with respect to delta (1–4 Hz) and theta (5–10 Hz) frequency bands was quantified by applying a zero-phase, fourth-order Butterworth band-pass filter and then calculating the amplitude from the absolute value of the Hilbert transform of the filtered signal, followed by smoothing (Gaussian kernel with σ = 5 s) and standardization (z-scoring). Periods for which the ratio of the amplitudes of theta and delta activity (the theta/delta ratio) remained above 5.0 for at least 20 s were classified as REM. Periods during which the theta/delta ratio remained below 2.0 for at least 20 s were classified as SWS (Extended Data Fig. 11b).
Detection of sweep and internal direction signals during sleep
To decode sweeps and internal direction from neural activity during sleep, we used tuning curves (LMT) from open-field sessions from the same recording day as the sleep session. Position was decoded separately for individual grid modules (see the ‘Single-module decoding’ section) because the correlation structure of grid cells across brain states may be preserved within but not across modules5,19. Because the theta rhythm is absent during SWS, we used local maxima in population activity as reference points for analysis of sweeps and direction signals in all brain states. To detect local maxima in the population activity, regardless of brain state, the spike counts of all internal direction cells were summed and smoothed with a Gaussian kernel (σ = 20 ms) before applying the Matlab function findpeaks with default parameters to detect peaks in the summed activity. Local maxima occurred at theta-rhythmic intervals while awake and during REM, and irregularly during SWS (Extended Data Fig. 11). Pairs of local maxima 2–250 ms apart were used to quantify directional alternation in all brain states, and all detected maxima were used to measure alignment between sweep and direction signals. Internal direction was taken as the decoded direction at the time of local maxima. To extract sweeps, we identified smooth sequences of decoded positions (sequences for which decoded position jumped less than 15% of grid spacing and changed direction by less than 2 radians between successive time bins) that occurred in windows centred around each of the local maxima in population activity. The windows extended 50 ms to either side of local maxima or to the edge of the neighbouring window. Because spatial representations were decoupled from physical movement during sleep (Extended Data Fig. 11g), sweep trajectories were referenced to the low-pass-filtered decoded trajectory (smoothed with a 100 ms gaussian kernel) and aligned to a ‘virtual head direction’ (low-pass-filtered decoded direction, σ = 1 theta cycle gaussian smoothing).
Simulation of an ideal sweep-generating agent
To test the hypothesis that alternating sweeps are controlled by an algorithm that maximizes the sampling of surrounding space, we simulated a sweep-generating agent that maximized environmental sampling by choosing sweep directions that minimized the overlap with previous sweeps.
First, we modelled the spatial coverage of a single sweep. Grid modules express sweeps at multiple spatial scales, so we reasoned that the total spatial coverage of a sweep may be considered as a sum of sweeps across individual grid modules (Extended Data Fig. 12a). A model sweep footprint was formulated, based on previous empirical observations of geometric relationships between single-module grid patterns16 and the geometric properties of sweeps in the present results (Fig. 1e–g). In brief, we summed the grid patterns of five idealized grid modules, with an inter-module scale ratio of √2 and gaussian-shaped grid fields with σ = 1/6 of each module’s spacing, at offsets from the origin corresponding to typical single-module sweep lengths of 1/3 of module spacing. The sum of sweeps across modules resembled a torch beam radiating outwards: as distance from the origin increased, the footprint broadened and decayed in intensity (Extended Data Fig. 12b). We approximated this shape by multiplying two simple spatial functions: an inverse distance function and an angular weighting function taken from a Von Mises distribution.
If we let d and θ denote the distance and direction, respectively, from the agent’s [x, y] position, xagent, to a location x in the environment (d = ||x − xagent|| and θ = arctan(x − xagent)), the intensity of the sweep footprint at location x for a chosen sweep direction α becomes:
$$f({\bf{x}})=1/{d}^{2}\exp ({\kappa }\;\text{cos}(\theta -\alpha )),$$
where κ is the angular concentration parameter of a Von Mises distribution, which determines the angular width of the sweep footprint. A value of κ = 5 was used initially (Fig. 6b) to reflect the empirically derived sweep shape (Extended Data Fig. 12b), but a parameter search revealed that stable alternation emerged across a range of values of κ (Extended Data Fig. 12d).
To run the simulation of sweeps on a linear track, we created an artificial scale-free 2D environment, binned into a 401 × 401 square grid. The agent was moved along a linear path at constant speed and instructed to generate a sweep every time step, by placing a sweep footprint in a specified direction. The cumulative trace of sweeps h at time t was computed by summing the footprints of previous sweeps:
$${h}_{t}({\bf{x}})=\mathop{\sum }\limits_{i=1}^{t-1}{f}_{i}\,({\bf{x}}).$$
The optimal sweep direction αoptimal at time t was chosen by finding the angle α that minimized the spatial overlap between the current sweep f and the cumulative trace of previous sweeps h:
$${\alpha }_{\text{optimal}}=\mathop{{\rm{argmin}}}\limits_{\alpha }\mathop{\sum }\limits_{i=1}^{N}f({{\bf{x}}}_{i},\alpha )h({{\bf{x}}}_{i}),$$
where the spatial overlap was calculated by multiplying the current sweep footprint and the cumulative sweep trace and summing across all N spatial bins.
Next, the simulation was run using the recorded behavioural trajectory of rats running in the open field. Real-world positions were mapped onto the agent’s simulated environment by setting the agent’s bin size to 1 cm and placing the open field at the centre of the bin grid. The agent’s time steps were yoked to the times of theta cycles in the experimental data, and for each theta cycle, the agent deployed the above algorithm to select the optimal sweep direction. To prevent the agent from being influenced by sweeps that occurred at similar locations far in the past, we introduced a temporal decay factor τ (range 0–1) that exponentially discounted the intensity of the cumulative coverage trace at each time step:
$$h({\bf{x}})=h({\bf{x}}){\tau }^{dt}.$$
The agent’s sweep-direction choices are determined solely by its previous decisions, so we called this the self-driving version of the model. We also formulated an empirically driven version, for which the agent was tasked with predicting the optimal sweep at each time step, given the directions of previous sweeps decoded from neural data. This was implemented by using the LMT, instead of the agent’s past sweep directions, to compute the cumulative sweep trace. Three rats were excluded from these analyses because internal direction could not be reliably estimated in these animals.
To investigate whether alternation also emerged when sweep directions were chosen independently for each grid module, we formulated a multimodule version of the simulation in which sweep footprints were simulated as three Gaussian functions (corresponding to three grid modules) whose width and sweep length increased with a geometric ratio of 1.5. At each time step, sweep directions for each module were chosen according to one of three updating rules: common, parallel or serial. In the common version, a single sweep direction was chosen for all modules, corresponding to the direction that minimized the overlap of the summed multimodule footprint with the summed multimodule trace. In the parallel and serial versions, sweep directions were chosen individually for each module, by minimizing the overlap between the single-module sweep footprint and the summed multimodule trace. In the parallel version, modules were updated simultaneously, meaning that modules were agnostic of each other’s sweep directions within same time step. In the serial version, modules were updated sequentially, by updating the summed multimodule sweep trace with the chosen sweep of one module before choosing the sweep direction for the next module.
Histology and recording locations
The rats received a lethal dose of pentobarbital, after which they were perfused intracardially with saline followed by 4% formaldehyde. The brains were extracted and stored in 4% formaldehyde, and were later cut in 30-µm sagittal or coronal sections using a cryostat. The sections were Nissl-stained with cresyl violet, and probe shank traces were identified in photomicrographs. In 14 animals, recording sites on the probes targeting the MEC–parasubiculum were aligned to the histological sections, as done previously5, by using as reference points the tip of the probe shank and the intersection of the shank with the brain surface. The aligned shank map was then used to calculate the anatomical locations of individual recording sites (Extended Data Fig. 1). Estimates of anatomical locations are subject to some degree of measurement error owing to the limited accuracy of the alignment process and the fact that units may be detected some distance away from the recording site.
Data analysis and statistics
Data analyses were performed with custom-written scripts in Matlab and Python. Clustering analyses of grid-cell modules and bursting subtypes of grid cells were conducted using the python package Scanpy87 and its dependencies (including numpy, pandas, scipy, scikit-learn and matplotlib). The LMT model and the functional-connectivity analyses were implemented by adapting publicly available code from ref. 46 and ref. 77, respectively. Statistical analysis was done using Matlab. Circular statistics were computed using the Circular Statistics Toolbox88. Results are reported as mean ± s.e.m. unless otherwise indicated. Statistical tests were non-parametric and two-tailed, unless otherwise indicated. The Mann–Whitney U-test was used for unpaired comparisons, and the Wilcoxon signed-rank test was used for paired comparisons. Pearson correlations were used unless otherwise indicated. Power analysis was not used to determine sample sizes. For each animal, the recording session with the best unit yield and behavioural performance was included in the study. The study did not involve any experimental subject groups, so random allocation and experimenter blinding did not apply and were not performed.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
