Research ArticleNEUROSCIENCE

Optogenetic pacing of medial septum parvalbumin-positive cells disrupts temporal but not spatial firing in grid cells

See allHide authors and affiliations

Science Advances  05 May 2021:
Vol. 7, no. 19, eabd5684
DOI: 10.1126/sciadv.abd5684


Grid cells in the medial entorhinal cortex (MEC) exhibit remarkable spatial activity patterns with spikes coordinated by theta oscillations driven by the medial septal area (MSA). Spikes from grid cells progress relative to the theta phase in a phenomenon called phase precession, which is suggested as essential to create the spatial periodicity of grid cells. Here, we show that optogenetic activation of parvalbumin-positive (PV+) cells in the MSA enabled selective pacing of local field potential (LFP) oscillations in MEC. During optogenetic stimulation, the grid cells were locked to the imposed pacing frequency but kept their spatial patterns. Phase precession was abolished, and speed information was no longer reflected in the LFP oscillations but was still carried by rate coding of individual MEC neurons. Together, these results support that theta oscillations are not critical to the spatial pattern of grid cells and do not carry a crucial velocity signal.


Medial entorhinal cortex (MEC) expresses multiple spatially modulated cell types. Among them, grid cells represent the environment by a characteristic hexagonal activity pattern (1). Characterized by their regular spacing, orientation, and phase, grid cells are believed to aid navigation based on self-motion cues (1, 2). Grid cells form distinct modules that share similar spacing and orientation (3). This population activity is coherent across environments (4), strongly suggesting that grid representations are not produced by single neuron computations but rather emerge from a continuous attractor network (CAN) (5). Notably, individual grid cells are not directly connected but communicate through a dense network of inhibitory neurons (68), supporting the idea that grid cells emerge through self-organizing principles dictated by underlying network connectivity. How grid cells support self-localization is, however, highly debated, and a prominent hypothesis builds on the notion of oscillatory interference (OI) (9, 10).

Rhythmic neural activity is proposed to be essential to synchronize activity for efficient population coding. Oscillations in the theta frequency range (6 to 10 Hz) are particularly prominent in the medial temporal lobe, including the MEC, of animals engaging with their environments. Theta activity has a proposed role in memory processing, as sequences of grid cells and place cells in the hippocampus and entorhinal cortex are consistently modulated by theta activity, possibly facilitating plasticity and strengthening connections within neural ensembles (11). In the entorhinal cortex, this can be observed in a subset of grid cells in superficial layers (LII and III) that show phase precession, where spikes fall progressively earlier in the theta cycle as the animal traverses through a grid field (12, 13). Phase precession was first described in hippocampal place cells (14) and is proposed to be a facilitating mechanism for sequence learning and an essential component for spatial information processing.

Theta oscillations in MEC are associated with the rhythmic activity in the medial septum and the diagonal band of Broca, together forming the medial septal area (MSA) (15, 16). The MSA contains multiple cell types that drive the pace-making activity (17), which is transmitted to the MEC by projections from glutamatergic, cholinergic, and GABAergic neurons (18). Lesions or pharmacological inactivation of MSA disrupt theta oscillations in MEC and impair grid cell activity (1921), suggesting a casual role for theta in grid cell activity. Inactivation of MSA affects, however, all downstream areas, not only the projections to MEC. Consequently, direct evidence linking theta oscillations driven by MSA to grid cell firing in freely moving animals is still lacking.

Several efforts have been made to establish computational models that may explain how grid cell activity arises. In one prominent model, the relationship between theta oscillations and phase precession has been suggested to underlie the formation of grid cell spatial representations through OI (9, 10). However, such single-cell models exclude the population representation that is evident in grid cells, while population representations are included in grid cell models based on CANs (22, 23). Lately, OI models have been merged with CANs in the so-called hybrid models, encompassing both population representation and phase precession (24, 25). When combined, these models are able to explain both ramping in membrane potential and phase precession as predicted by CANs and OI models, respectively (24, 26). To map neural activity to stable spatial representations, some form of path integration is required. One of the main disagreements between grid cell models is how location is encoded, updated, and decoded (27, 28). In hybrid models, path integration is accomplished by means of phase precession, given by a phase difference between velocity-controlled oscillators and a baseline oscillation typically thought of as being driven by MSA input (24, 25). As a consequence, the model predicts a causal relationship between phase precession and the grid cell spatial representation, but this remains to be tested.

The GABAergic part of the MSA projection predominantly terminates on inhibitory interneurons in MEC and is suggested to be the main generator of rhythmic activity through coordination of the local inhibitory circuitry in MEC and the hippocampus (2931). To test the prediction that phase precession underlies path integration of grid cells, we disrupted the endogenous theta [6 to 10 Hz, containing the peak of the power spectrum density (PSD)] and phase precession in MEC using optogenetic perturbation of the MSA input. We hypothesized that by controlling the spiking activity of grid cells through disinhibition, thus locking the spike-phase relation to the stimulation frequency, phase precession would also be abolished. We restricted the perturbation to the GABAergic part of the MSA projections by taking advantage of the Cre-Lox system using a PvalbCre rat line (32). This allowed us to selectively target inhibitory cell types expressing parvalbumin (PV), the main GABAergic cell type in MSA (33, 34). To test a relation between phase precession and path integration, we chose stimulation frequencies at 11 and 30 Hz, which are at the high end and far outside the endogenous theta ranges, respectively. Optogenetic activation of PV+ neurons in MSA induced artificial local field potential (LFP) oscillations in MEC and significantly reduced endogenous theta power. Activation of MSA PV+ cells caused robust activation of grid cells through disinhibition. Thus, both LFP oscillations and single-unit phase relations followed the stimulation frequency and effectively disrupted grid cell phase precession. The perturbation did not affect the location of the grid fields. Moreover, there was a prominent reduction in information rate of grid cells during stimulation, suggesting a possible impairment in the spatial readout downstream of the grid cell circuitry. Last, spike-speed correlations remained stable during stimulation, while both frequency and power of LFP theta were disrupted. The findings show a dissociation between theta oscillations and grid cell spatial activity patterns. This supports theory in which theta oscillations do not cause grid cell firing patterns.


Specific targeting of PV+ cells in MSA

To target the GABAergic projection from MSA to MEC, channelrhodopsin (ChR2) was selectively expressed in PV+ cells in MSA by injecting a viral construct in four animals from a PvalbCre rat line (Fig. 1A) (32). Immunohistochemical labeling of ChR2 confirmed that virus expression was largely restricted to MSA (Fig. 1B and fig. S1), and costaining with a PV antibody showed that expression was selective to PV+ cells. We did not observe any ChR2-expressing cells that were not immunopositive for PV. Terminals expressing ChR2 were prominent in both MEC and the hippocampus, and we observed ChR2-expressing boutons contacting PV+ cell soma locally in MEC (Fig. 1C).

Fig. 1 PV+ cells in MSA selectively express channelrhodopsin (ChR2) after injection of virus in PvalbCre rats.

(A) Illustration of a rat brain seen slightly from above on the left side highlights the MSA in purple and MEC in blue/cyan, with long-range projections from MSA to MEC (green). These projections target all layers of MEC, from the dorsal to the ventral area. A viral construct carrying ChR2 was injected in MSA. (B) Three coronal sections from a representative animal show the extent of the virus expression (green) in MSA at three different anterior-posterior positions. Expression covered large parts of MSA. (C) Virus expression in MSA was restricted to PV+ cells (red). White arrowheads mark overlap between virus (ChR2) and PV+ cells (PV) (top row). Virus-expressing projections were found in MEC, parasubiculum (PaS), and all regions of the hippocampus [CA1, CA3, and dentate gyrus (DG)] (bottom row). The area of MEC chosen for the magnified images is indicated by a small square. ChR2-labeled septal projections target PV+ cells in MEC (small outline). (D) Illustration of experimental setup. PvalbCre rats were implanted with optic fiber in MSA and recording electrodes in MEC (two animals also had optic fibers together with recording electrodes in MEC). Blue laser light (470 nm) was used to activate PV+ cells in MSA at two different frequencies, 11 and 30 Hz.

To determine whether we could pace oscillations in MEC, we stimulated PV+ cells in MSA at two different frequencies, 11 and 30 Hz, while recording LFP and single-unit activity in MEC (Fig. 1D). All single-unit recordings were located in superficial layers of MEC (LII and III) (fig. S2); see table S1 for cell counts.

Stimulating PV+ cells in MSA paces LFP oscillations in MEC and drives grid cells through disinhibition

We stimulated either PV+ cell bodies in MSA (Fig. 2A) or PV+ terminals in MEC (Fig. 2D) while recording LFP and single units in MEC from PvalbCre rats during exploration of a 1 m–by–1 m open field. Experiments where we stimulated in MSA started with a 10-min baseline session (Baseline I), followed by 10-min stimulation at 11 Hz, 10 min without stimulation (Baseline II), and ending with 10-min stimulation at 30 Hz. These four sessions were separated from each other by a 10-min break where the animal rested in their home cage. Immediately after stimulus onset, the endogenous theta (6 to 10 Hz) was reduced and LFP oscillations shifted to the stimulation frequency (Fig. 2A and fig. S3). Thus, the activity of stimulated PV+ cells in MSA controlled the frequency oscillations in MEC. Average PSD for all baseline or stimulation sessions showed a similar pattern with a peak at 8 Hz in baseline sessions and a new peak corresponding to the stimulation frequencies of either 11 or 30 Hz (Fig. 2B). Endogenous theta power, measured by the relative height of the peak PSD, was significantly reduced during stimulation [pooled sessions, means ± SEM; Baseline I versus 11 Hz: Baseline I, 6.3 ± 0.67; 11 Hz, 0.23 ± 0.052 (n = 44, W = 2.0, P = 8.7 × 10−9); Baseline II versus 30 Hz: Baseline II, 6.3 ± 0.74; 30 Hz, 0.049 ± 0.054 (n = 32, W = 0.0, P = 8.0 × 10−7, Wilcoxon signed-rank tests)]. The relative power in the stimulated band was larger for 30-Hz stimulation [pooled sessions, means ± SEM; 11 Hz versus 30 Hz: 11 Hz, 18 ± 2.4; 30 Hz, 85 ± 14 (n = 32, W = 30, P = 1.2 × 10−5, Wilcoxon signed-rank test)] (Fig. 2C).

Fig. 2 Optogenetic pacing of PV+ cells in MSA abolish endogenous theta in MEC.

(A) Illustration of a presumed MSA-to-MEC connection showing that PV+ cells (PV) in MSA terminate on interneurons (IN) in MEC that together with grid cells (GC) form a recurrent circuit. Optogenetic activation of MEC projecting PV+ cells in MSA (blue shading) caused a frequency shift in oscillations in MEC as shown by time-frequency representations of LFP during open field exploration. (B) Average PSD for all baseline and stimulus sessions. (C) Cumulative density plots of relative theta power and relative power in the stimulated band show that the power of the endogenous theta was strongly reduced for both stimulation frequencies. (D) Illustration of optogenetic stimulations in MEC. Optogenetic activation of MSA PV+ cell terminals in MEC did not reduce the endogenous theta frequency. Time-frequency representations of LFP in two simultaneously recorded hemispheres where one is stimulated with 11 Hz. (E) Similar to (B). The endogenous theta remained when stimulating PV+ cell terminals in MEC at 11 Hz. (F) Cumulative density plots of relative theta power show that the power of endogenous theta was reduced when stimulating PV+ cell terminals in MEC. (G) Raster plot of single-unit spiking responses of one grid cell and one narrow spiking inhibited (NSi) cell. Probability density for the whole population of grid cells and NSi cells shows fast inhibition of NSi cells followed by grid cell activation. Thirty-hertz stimulation led to a second activation of both cell types after approximately 20 ms. (H) Raster plots from one grid cell and one NS cell during optogenetic stimulation in MEC. Probability densities of the two units are overlaid in the bottom graph to illustrate difference in response times.

Optogenetic stimulation within MEC also started with a baseline session (Baseline I), followed by two 10-min sessions with 11-Hz stimulation, one in each hemisphere and finalized by 10 min without stimulation (Baseline II). Stimulus onset led to an immediate response in theta frequency (Fig. 2D) but did not disrupt the endogenous theta (Fig. 2, E and F). Nevertheless, the relative power of the endogenous theta band was significantly reduced between baseline and stimulation [pooled sessions, means ± SEM; Baseline I versus 11 Hz: Baseline I, 9.7 ±1.4; 11 Hz, 4.8 ± 7.8 × 10−1 (n = 8, W = 1.0, P = 1.6 × 10−2, Wilcoxon signed-rank test)].

Since inhibitory neurons in MEC constitute the main postsynaptic targets of GABAergic projections from MSA (30), we investigated the response to stimuli in putative inhibitory neurons together with responses in grid cells (Fig. 2G). We separated narrow spiking (NS), putative inhibitory neurons, from broad spiking (BS), putative principal cells based on spike waveform parameters (fig. S4). A quarter of NS units (31 of 126) responded with a rapid decrease in firing rate within 3 to 6 ms after stimulation onset, followed by a visible increase in spiking activity [we termed these cells NS inhibited (NSi)]. The majority of NS units (87 of 126) were not inhibited by the stimulus [termed NS not inhibited (NSni)] (fig. S5). Both NSi and NSni cells showed a second phase of activation, likely due to the putative disinhibition feeding back to inhibitory cells after stimulation (fig. S5A).

To label neurons as grid cells, we used the 95th percentile of gridness and information rate for each unit from a shuffled distribution of spikes as the threshold together with a minimum gridness score of 0.2. Shortly after NSi cells were inhibited, grid cells showed consistent, robust excitatory responses with an average peak response of 8 ms after stimulation onset. These results strongly indicate that grid cells were activated through disinhibition (Fig. 2G) since all observed grid cell NSi pairs (recorded on same drive) showed NSi inhibition preceding grid cell response (means ± SEM; 0.004 ± 0.0002 s). Following the activation of grid cells, NSi and NSni cells were also activated, suggesting that the widespread recurrent connectivity in MEC leads to large-scale activation of the network in response to the initial inhibition. The 11- and 30-Hz stimulation elicited different responses in grid cells, and 11 Hz generated stronger response in peristimulus time histogram (PSTH) (Fig. 2G) when compared with 30 Hz [probability of spiking: weighted means ± SEM, effect size β and P value from linear mixed effects model (LMM); 11 Hz, 0.12 ± 0.01 (n = 58); 30 Hz, 0.047 ± 0.006 (n = 33, β = 0.07, P = 0.0045)]. Moreover, 30 Hz tended to generate a longer phase of increased activity [full width at half maximum: weighted means ± SEM, effect size β and P value from LMM; 11 Hz, 0.005 ± 0.00025 s (n = 58); 30 Hz, 0.01 ± 0.001 s (n = 33, β = 0.004 s, P = 0.1)]. After the putative disinhibition in grid cells, a brief period of significant inhibition (15 to 25 ms after stimulus onset) occurred during 11-Hz stimulation compared to a baseline period (−5 to 5 ms relative to stimulus onset), while 30 Hz elicited a second bout of increased activity, although not significantly different from baseline [PSTH weighted means ± SEM, effect size β and P value from LMM; Baseline 11 Hz, 18 ± 0.65 (n = 58); postresponse, 8.3 ± 1.0 (n = 58, β = 10, P=4.5 × 10−5); Baseline 30 Hz, 8.3 ± 0.38 (n = 33); postresponse, 9.2 ± 053 (n = 33, β = 0.9, P = 0.63)].

Because of the interconnectivity between GABAergic, glutamatergic, and cholinergic neurons in MSA (17), stimulating PV+ cells could potentially alter the activity of glutamatergic and/or cholinergic neurons locally within MSA. The increased responses of grid cells after stimulation could therefore be a result of delayed activation of glutamatergic or cholinergic neurons in MSA caused by PV+ cell stimulation. To test whether grid cell responses resulted from disinhibition through local interneurons in MEC or increased excitatory input from MSA, we compared responses in grid cells and NS units following stimulation of PV+ cell terminals locally within MEC. Stimulating terminals in MEC produced responses in grid cells with similar delay as stimulating in MSA (Fig. 2H). This similarity in response supports that grid cells were likely activated through disinhibition following stimulations of PV+ cells and not by increased excitation of glutamatergic or cholinergic projections from MSA.

Grid cells follow the stimulation frequency, abolishing phase precession

Theta activity has been linked to grid cell spatial representations in both experimental (20, 21, 35) and theoretical work (9, 24, 25). By reliably controlling oscillation frequencies in MEC when stimulating in MSA, we investigated the effect on grid cells’ response patterns.

First, we looked at the PSD of spikes to assess the theta modulation of the firing rate during stimulation (Fig. 3A). The relative power of this theta rhythmicity was significantly reduced during stimulation [weighted means ± SEM, effect size β and P value from LMM; Baseline I versus 11 Hz: Baseline I, 1.3 ± 0.1 (n = 63); 11 Hz, 0.2 ± 0.02 (n = 56, β = 1.0, P = 3.1 × 10−10); Baseline II versus 30 Hz: Baseline II, 1.3 ± 0.1 (n = 46); 30 Hz, 0.5 ± 0.1 (n = 35, β = −0.67, §§P = 5.6 × 10−10)] (Fig. 3B), suggesting that the spike train was no longer modulated by theta during the stimulation.

Fig. 3 Grid cells stop phase precessing during MSA stimulation.

(A) Theta rhythmicity estimated by the power spectral density of the neuronal firing rate obtained by kernel estimation with a Gaussian kernel of width 10 ms. (B) Cumulative density of the peak theta rhythmicity (6 to 10 Hz) divided by the average power of 1-Hz wide adjacent bands. (C) Spike LFP coherence in grid cells for baseline sessions compared to 11- and 30-Hz stimulation. Grid cells showed strong spike coherency to the endogenous theta frequency (around 8 Hz) in baseline sessions, which was reduced during stimulation. (D) Cumulative density plots show that grid cells significantly reduced spike coherency at the endogenous theta frequency during both 11- and 30-Hz stimulation. Both stimulated bands caused similar peak coherence. (E) Polar plot displaying theta phase preference and vector length of all grid cells during baseline (endogenous theta) and stimulation sessions (stimulated band). (F) Swarm plot of P values (only neurons with P value below 0.10 are shown) for neurons phase precessing (circular linear correlation r < 0) during Baseline I, 11 Hz, Baseline II, and 30 Hz sessions. We found no significant phase precessing neurons during stimulation as indicated by no points below the black line (P < 0.01). (G) Rows show color-coded rate maps from two grid cells during baseline followed by an 11-Hz stimulation session. Rat and unit numbers are indicated to the left, and peak rate and gridness are indicated above. Raster plots represent spike phase relative to LFP (band-pass filtered at 6 to 10 Hz) versus distance traveled through grid fields. The regression line indicates that the spike phase predicts the position of the animal within the field. P value of phase precession is shown to the right. Both example units ceased phase precessing during stimulation, shown by the loss of significant phase precession (P < 0.01).

Next, 11- and 30-Hz stimulation caused a reduction in grid cell spiking coherence with the endogenous LFP theta frequency [Fig. 3, C and D; weighted means ± SEM; Baseline I versus 11 Hz: Baseline I, 0.24 ± 0.026 (n = 63); 11 Hz, 0.069 ± 0.0063 (n = 56, β = 0.18, P = 8.8 × 10−9); effect size and P value, respectively, from LMM; Baseline II versus 30 Hz: Baseline II, 0.28 ± 0.026 (n = 46); 30 Hz, 0.12 ± 0.037 (n = 35, β = 0.16, P = 0.09); effect size and P value, respectively, from LMM]. When comparing stimulated frequency, similar peak coherence was found [11 Hz versus 30 Hz: means ± SEM; 11 Hz, 0.47 ± 0.03 (n = 58); 30 Hz, 0.43 ± 0.04 (n = 33, β = 0.054, P = 0.48); effect size and P value, respectively, from LMM]. Both populations of NS cells showed significant reduction in spike LFP coherence with endogenous theta and increased coherence with the stimulation frequency (fig. S6, B to E).

To further assess how grid cells spiked relative to LFP theta phase, we computed the resultant vector length of phase preference [Fig. 3E; weighted means ± SEM; Baseline I versus 11 Hz: Baseline I, 0.23 ± 0.013 (n = 63); 11 Hz, 0.038 ± 0.0033 (n = 56, β = 0.2, P = 2.5 × 10−6); effect size and P value, respectively, from LMM; Baseline II versus 30 Hz: Baseline II, 0.24 ± 0.016 (n = 46); 30 Hz, 0.018 ± 0.0031 (n = 35, β = 0.19, P = 1.1 × 10−6); effect size and P value, respectively, from LMM]. When comparing stimulated frequency, similar resultant vector length was found [11 Hz versus 30 Hz: means ± SEM; 11 Hz, 0.25 ± 0.016 (n = 58); 30 Hz, 0.26 ± 0.018 (n = 33, β = 0.022, P = 0.7); effect size and P value, respectively, from LMM]. Here, n indicates pooled unique cell count, and baselines were only selected when preceding a stimulation session. Results indicate that both 11 and 30 Hz caused grid cells to lose their phase locking to theta and a large fraction of grid cells to phase lock to the stimulation frequency (Fig. 3E, right), although not consistently to a preferred angle across neurons.

These findings suggested that grid cell phase precession was disrupted during stimulation sessions. We therefore identified all grid cells with significant phase precession or recession (spikes appearing progressively later in each theta cycle) relative to the endogenous theta frequency band (6 to 10 Hz) during baseline sessions (Baseline I and Baseline II). We then compared this to the same frequency band in the stimulation sessions (11 and 30 Hz, respectively) (Fig. 3G) using the methods described in (13, 36). To reliably identify loss of phase precession/recession, we performed two separate analyses, using either spatial field detection or firing rate threshold to detect field entry (described in Methods). From these detection methods, we assessed pooled field crossings and single run crossings in accordance with (13). Moreover, to verify the analysis, we reanalyzed publicly available data from the work of Sargolini and co-workers (37) and performed an additional control to assess phase precession relative to a “random” band in the LFP (20 to 25 Hz). With these data at hand, we could reliably detect loss of phase precession by using the P values outlined in (36) and reproduce the results from the work of Reifenstein et al. (13). Comparing phase precession between the theta band and the random band showed that in pooled field crossings, we reliably detected loss of phase precession by thresholding the P value. In contrast, with single run trials, this detection method was hampered by false positives, likely because of the low spike count (fig. S7).

For all recorded units, no grid cells displayed significant phase precession during any of the stimulation sessions (field detection, Fig. 3F; rate detection, fig. S7F), showing that phase precession ceased during optogenetic pacing [the number of significant phase precessing neurons (P < 0.01, R > 0.1): Baseline I, 9 (as shown in figure 3 f) of 63; 11 Hz, 0 of 56; Baseline II, 9 (as shown in figure 3 f) of 46; 30 Hz, 0 of 35].

Grid cells remain spatially stable during MSA stimulation but show increased out-of-field spikes during disinhibition

At first observation, we saw that the spatial firing pattern of grid fields remained remarkably stable between baseline and stimulation sessions (Fig. 4A and fig. S8). However, when inspecting rate maps, we noticed that many grid cells showed either increased or decreased firing rates during stimulation. To test whether this alteration was present, we first compared pairwise relative changes in spatial shift for each grid cell from a baseline session to a stimulation session (Fig. 4B). To assess the magnitude of change, we created control populations using the first half of a baseline session compared to the second half of the session (Baseline Ia versus Baseline Ib and Baseline IIa versus Baseline IIb). The comparison of spatial shift amounted to the following (weighted means ± SEM, effect size β and P value from LMM): Baseline Ia to Baseline Ib, 0.022 ± 0.0013 (n = 63); Baseline I to 11 Hz, 0.022 ± 0.0021 (n = 21, β = 0.001, P = 0.87); Baseline I to Baseline II, 0.022 ± 0.002 (n = 21, β = 0.023, P = 0.22); Baseline IIa to Baseline IIb, 0.019 ± 7.0 × 10−5 (n = 46); Baseline II to 30 Hz, 0.018 ± 0.0026 (n = 16, β = 1.4 × 10−4, P = 0.99). Note that relative measures require the same neurons to be present in subsequent sessions, therefore reducing the number of units (n) included in this analyses.

Fig. 4 Grid cells are spatially stable during optogenetic stimulation of MSA.

(A) Color-coded rate maps (top) and running path with spikes superimposed (bottom) from three grid cells recorded in all four sessions of 10 min each. Rate maps are adjusted with color coding corresponding to the heat color of the first recording session. Peak rates and gridness are denoted above each rate map. Rat number and unit ID are indicated to the left of the rate maps. (B) Scatter plots showing relative changes in spatial shift of grid cells followed across multiple recording sessions. There were no detectable significant changes compared to Baseline Ia or IIa versus Baseline Ia or IIb (first half versus second half of Baseline I or Baseline II, respectively). (C) Cumulative density plots of average rate, maximum rate, gridness, spatial information, and spatial information specificity. There were no detectable differences in either rate measures; however, there was a nonsignificant change in gridness between baseline and stimulation and a significant change in spatial information rate when comparing Baseline I and 11 Hz. Spatial information specificity showed no significant change. (D) Raster plot from an example grid cell where 0 marks 11-Hz stimulation onset. Spikes in the preresponse (Pre: green/purple; −5 to 5 ms) and response period (Resp: orange/pink; 5 to 11 ms) and the between pulses (Between: blue; 15 to −5 ms) are marked. The rat’s running trajectory of the recording session with spikes from Pre, Resp, and Between periods superimposed. Black circles indicate outline of identified grid fields. (E) Similar to (D) but showing 30-Hz stimulation. Violin plots [minimum to maximum and median (black line)] show increased out-of-field spiking activity during response periods. Width of violin corresponds to number of samples for each value. *P < 0.05 and ***P < 0.001; LMM test. ns, not significant.

Next, we calculated average rate, maximum rate, and gridness of all recorded grid cells (Fig. 4C). We found no significant differences from Baseline I to 11 Hz or Baseline II to 30 Hz [weighted means ± SEM, effect size β and P value from LMM; average rate: Baseline I, 9.8 ± 0.93 (n = 63); 11 Hz, 11 ± 1.0 (n = 56, β = 0.6, P = 0.69); Baseline II, 10 ± 1.0 (n = 46); 30 Hz, 8.5 ± 1.1 (n = 35, β = 0.51, P = 0.6.8); Baseline I to Baseline II, β = 0.31, P = 0.82; maximum rate: Baseline I, 36 ± 2.1 (n = 63); 11 Hz, 36 ± 2.2 (n = 56, β = 1.3, P = 0.76); Baseline II, 42 ± 2.4 (n = 46); 30 Hz, 37 ± 2.5 (n = 35, β = 1.7, P = 0.72); Baseline I to Baseline II, β = 3.3, P = 0.34; gridness: Baseline I, 0.37 ± 0.06 (n = 63); 11 Hz, 0.4 ± 0.05 (n = 56, β = 0.028, P = 0.85); Baseline II, 0.53 ± 0.04 (n = 46); 30 Hz, 0.57 ± 0.05 (n = 35, β = 0.027, P = 0.86); Baseline I to Baseline II, β = 0.084, P = 0.74]. See also the Supplementary Materials for more grid cell statistics (figs. S8 and S9 and tables S2 and S3).

Unexpectedly, neither stimulation frequencies caused any significant change in spatial position and gridness or change in mean or maximum firing rate, indicating that the artificially induced oscillations did not alter spatial position or the firing properties of grid cells. We observed a small, nonsignificant effect on gridness between baseline and stimulation. We also performed these analyses for BS nongrid cells and the NSi and NSni cells (see fig. S10 for details).

Despite spatially stable fields during stimulation, there was a reduction in spatial information between Baseline I to 11 Hz [weighted means ± SEM, effect size β and P value from LMM; Baseline I, 1.4 ± 0.086 (n = 63); 11 Hz, 0.93 ± 0.069 (n = 56, β = 0.4, P = 0.017); Baseline II, 1.3 ± 0.1 (n = 46); 30 Hz, 1.1 ± 0.11 (n = 35, β = 0.027, P = 0.95); Baseline I to Baseline II, β = 0.097, P = 0.043] (Fig. 4C). However, when controlling for rate (spatial information specificity), there was no detectable reduction [weighted means ± SEM, effect size β and P value from LMM; Baseline I, 0.25 ± 0.034 (n = 63); 11 Hz, 0.22 ± 0.045 (n = 56, β = 0.037, §P = 0.76); Baseline II, 0.23 ± 0.034 (n = 46); 30 Hz, 0.24 ± 0.038 (n = 35; β = 0.02, P = 0.33)]. This lack of effect in information specificity indicates that the spatial information in each spike was preserved (most spikes fall within fields) but that the information per time was decreased; however, we did not ascertain the difference in effect size.

Since stimulation induced a putative disinhibition in grid cells followed by a period of inhibition, we suspected that more spikes would fall outside grid fields during periods of increased activation and reduce the spatial information during stimulation sessions, similar to the findings from Zutshi et al. (38). To assess this suspicion, we separated each stimulation phase into a preresponse phase, termed Pre (−5 to 5 ms relative to stimulus onset); a response phase, termed Resp (5 to 11 ms after stimulus onset); and lastly a period between pulses, termed Between (15 ms after stimulus onset to 5 ms before the next stimulus). We then calculated the relative percentage of out-of-field spikes during each of these phases (Fig. 4D). We found increased out-of-field spikes during the response phase for both stimulation frequencies [weighted means ± SEM, effect size β and P value from LMM; 11 Hz: Pre, 0.57 ± 0.021 (n = 56); Resp, 0.63 ± 0.018 (n = 56); Between, 0.55 ± 0.021 (n = 56). Pre versus Resp: β = 0.06, P = 2.7 × 10−9; Pre versus Between: β = 0.024, P = 9.4 × 10−4; Resp versus Between: β = 0.081, P = 2.7 × 10−17; 30 Hz: Pre, 0.54 ± 0.019 (n = 35); Resp, 0.54 ± 0.019 (n = 35); Between, 0.49 ± 0.022 (n = 35); Pre versus Resp: β = 0.0036, P = 0.84; Pre versus Between: β = 0.044, §§P = 0.15; Resp versus Between: β = 0.052, P = 0.028]. This suggests that out-of-field spiking was caused by the short time window of enhanced spiking after each stimulation pulse. Together, this shows that the oscillatory septal input can modulate the temporal firing of grid cells within stable fields but is not determinant of field position.

Speed correlations of grid and NSi cells were intact during MSA stimulation while disrupted in LFP oscillations

Several computational models predict that to perform path integration, grid cells need information about the animals’ speed (6, 22, 23, 3941). It has been shown that there are at least two individual sources of speed information that can be present in single neurons: one based on firing rate and the other based on oscillatory frequency (42). The oscillatory frequency decreases if input from MSA is impaired; hence, septal input may be critical for intact speed representation of neurons in MEC.

On the basis of these assumptions, we tested whether pacing MSA PV+ cells would disrupt rate-based speed representation of individual units in MEC. We computed speed scores for individual grid and NSi cells, given by the linear relationship between firing rate and the animals’ running speed (43), and compared these across baseline and stimulation sessions (Fig. 5A). We found no significant change in speed scores during stimulation for grid cells [weighted means ± SEM, effect size β and P value from LMM; Baseline I, 0.11 ± 0.011 (n = 63); 11 Hz, 0.094 ± 0.012 (n = 56, β = − 0.013, P = 0.53); Baseline II, 0.087 ± 0.0080 (n = 46); 30 Hz, 0.074 ± 0.0099 (n = 35, β = 0.013, §§P = 0.38); Baseline I and Baseline II, β = 0.0094, P = 0.68]. Moreover, we found no significant effect on speed scores for NSi cells [weighted means ± SEM, effect size β and P value from LMM; Baseline I, 0.15 ± 0.04 (n = 12); 11 Hz, 0.18 ± 0.02 (n = 25, β = 0.018, P = 0.69); Baseline II, 0.17 ± 0.025 (n = 17); 30 Hz, 0.21 ± 0.019 (n = 32, β = 0.035, P = 0.3); Baseline I and Baseline II, β = 0.017, P = 0.54] (Fig. 5B, fig. S5, and tables S2 and S3).

Fig. 5 Neuronal speed modulation is stable despite disrupted speed modulation in theta.

(A) Example of speed modulation in one grid cell and one NS cell (NSi) during baseline and stimulation sessions. Speed score is denoted above each speed graph. Vertical axis represents normalized firing rate. (B) Violin plots show speed scores in grid cells (top) and NSi cells (bottom) between baseline and stimulation sessions. Neither 11- nor 30-Hz stimulation caused any significant change in speed scores for either cell type. (C) Violin plots showing paired comparisons of running speed of animals in all recording sessions. Running speed decreased slightly from Baseline I to 11 Hz, from Baseline I to Baseline II, and from Baseline II to 30 Hz. (D) Correlation between running speed and theta peak frequency was strong in both Baseline I and Baseline II but disrupted during 11- and 30-Hz stimulation. (E) Frequency score, as represented by the correlation between running speed and peak frequency. This was significantly reduced from Baseline I to 11 Hz and from Baseline II to 30 Hz. (F) Correlation between running speed and theta power was strong in both Baseline I and Baseline II but disrupted during 11- and 30-Hz stimulation. (G) Power score, as represented by the correlation between running speed and power. This was significantly reduced from Baseline I to 11 Hz and from Baseline II to 30 Hz. *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001, Wilcoxon signed-rank test.

Activity in glutamatergic cells in MSA is strongly linked to locomotor behavior and theta frequency (44). Therefore, we also tested whether the running speed of the animals was altered during sessions of theta pacing. We observed a small but significant decrease in running speed [means ± SEM, Wilcoxon signed-rank tests; Baseline I, 0.16 ± 0.0054; 11 Hz, 0.15 ± 0.004 (n = 22, W = 47, P = 0.017); Baseline II, 0.15 ± 0.004; 30 Hz, 0.14 ± 0.005 (n = 16, W = 18, P = 0.0097)]. There was also a small decrease from Baseline I to Baseline II (n = 15, W = 17, P = 0.015) and a prominent decrease from 11 to 30 Hz (n = 16, W = 0.00, P = 4.4 × 10−4). The subsequent decrease in running speed from one session to the next, including from Baseline I to Baseline II, may suggest that the animals grew increasingly tired from running or lost interest in the chocolate treats, rather than locomotor behavior being affected by stimulating PV+ cells in MSA (fig. 5C).

Self-motion–based modulation of theta is also suggested to be a central mechanism for creating grid cell firing patterns because of the observed lack of grid cell activity during passive transport (35). Since we observed no disruption of grid cell patterns during MSA pacing, we tested whether LFP theta frequency was modulated by the animals’ running speed during stimulation. As expected, we observed a strong correlation between running speed and theta frequency for all units in baseline sessions. However, the correlation was completely disrupted both during 11- and 30-Hz stimulation sessions [means ± SEM, Wilcoxon signed-rank test; Baseline I, 0.18 ± 0.017; 11 Hz, −0.0064 ± 0.0108 (W = 16, P = 2.27 × 10−8, n = 44); Baseline II, 0.22 ± 0.02; 30 Hz, 0.012 ± 0.014 (W = 9, P = 1.86 × 10−6, n = 32)] (Fig. 5, D and E).

Under normal conditions, speed information is reflected in both LFP theta frequency and the power of these oscillations. Therefore, it could be that the speed-power relationship was sustained although the speed-frequency relationship was broken as previously observed (45). While we observed a strong correlation between running speed and theta power in baseline sessions, this was completely disrupted during both 11- and 30-Hz stimulation sessions [means ± SEM, Wilcoxon signed-rank test; Baseline I, 0.15 ± 0.02; 11 Hz, −0.02 ± 0.03 (W = 110, P = 7.02 × 10−6, n = 44); Baseline II, 0.11 ± 0.02; 30 Hz, 0.05 ± 0.02 (W = 152, P = 0.0362, n = 32)] (Fig. 5, F and G, and fig. S11). These results suggest that speed modulation of theta is dissociated from speed modulation of individual MEC neurons and is not needed for stable grid cell firing patterns.


In this study, we investigated the hypothesis of a causal relationship between theta oscillations and spatial representation of grid cells, postulated by experimental findings (20, 21, 35) and grid cell hybrid models (24, 25). We show that the spatial and the temporal firing pattern of grid cells can be dissociated. Experimentally induced oscillations disrupted endogenous theta activity in MEC, impaired phase precession, and speed modulation of LFP, but the spatial representations of grid cells and speed modulation remained intact.

While grid cells were strongly modulated by the imposed oscillatory frequency, this did not disrupt their spatial firing patterns. We did however see a reduction in spatial information rate that was confined to the short time period where grid cells were activated by the optogenetic stimulation of PV+ cells in MSA. Grid cells were likely activated through disinhibition, thus providing additional support for an inhibitory circuitry where PV+ cells in MSA provide monosynaptic input to local interneurons in superficial layers of MEC (29), which, in turn, are tightly connected to grid cells (6, 39). This inhibitory circuitry seems determinant for theta oscillations in MEC. The data adds to a growing body of evidence suggesting that the GABAergic/PV+ neurons in MSA are responsible for pacing theta oscillations in MEC and hippocampus. Grid cells showed substantial phase locking to the stimulation frequency, resulting in impaired grid cell phase precession. This demonstrates that grid cells do not rely on phase precession to maintain stable spatial representations.

In this study, we aimed to fully disrupt the endogenous theta to assess its role for grid cell firing properties. To achieve this, we induced a strongly synchronized activity via optogenetic stimulation of PV+ cells, which completely abolished phase precession and theta phase locking in grid cells through putative disinhibition and disrupted the temporal modulation of the LFP theta. Because of the high amplitude and regularity of the stimulation pulses [as seen by the large harmonic components in Fig. 2 (B to E)], the resulting oscillation is not comparable to the endogenous theta oscillations. Given that we aimed to assess the robustness of grid cells, we took advantage of this nonphysiological oscillation paradigm, which allowed us to show that neither phase precession nor speed modulation of the theta oscillation is crucial for grid cells to retain their spatial properties. The abolished phase precession in grid cells, however, might not be a direct consequence of the drive of MSA but rather caused by the strong stimulation-induced synchrony, which would cancel any temporal information during forced spikes. In addition, the strong stimulation regime might change the circuit dynamics, possibly biasing interpretation of mechanistic underpinnings of neuronal oscillations presented in this study. To better assess how the MSA modulates phase precession and phase locking, a stimulation regime with nonsquare pulses, temporal jitter in interstimulus interval, and varying stimulation amplitude might be a better choice. This might help create an artificial septal non-oscillatory input or a septal input with nonstable frequencies that is observed in bats (46). Further studies involving a more physiological stimulation are therefore needed to additionally characterize oscillatory phenomena in MSA and MEC.

Previous experiments show that pharmacological inactivation of MSA disrupts both the spatial firing pattern of grid cells and theta oscillations in MEC and hippocampus (20, 21). This has been taken as evidence to support that theta oscillations are critical to grid cell activity. However, pharmacological inactivation of MSA also silences inputs to the hippocampus and the parahippocampal areas, and hippocampal inactivation has been shown to disrupt grid cell firing (47). Inactivation of MSA also disrupts parahippocampal activity, which, in turn, may impair grid cell activity. It is therefore unresolved whether the impaired grid cell patterns after MSA inactivation are caused by direct input to MEC or indirectly by disturbed hippocampal activity. Given the low temporal and spatial resolution of pharmacological manipulations, it is challenging to directly relate interventions to functional properties of grid cell activity in such previous experiments. Pharmacological manipulations also lack the necessary cell type specificity to assess the role of separate parts of the septal-entorhinal projection. Thus, the cell type–specific optogenetic activation in the PvalbCre rat line used in the presented investigation likely overcomes both these limitations.

Pyramidal neurons of the hippocampus provide rhythmic feedback to MSA and might thereby support the maintenance of theta rhythm (48, 49). One might thus argue that disturbed hippocampal activity caused by septal inactivation leads to loss of theta, which then disrupts grid cell activity. In contrast to inactivation, various alterations of MSA activity show little effect on the spatial firing of hippocampal place cells, despite disrupting theta oscillations in the hippocampus. For example, Zutshi and colleagues (38) have demonstrated that optogenetic pacing of theta does not disrupt place fields but decreases the memory performance of mice in an alternating eight maze. Furthermore, another study shows that disinhibition of MSA using a γ-aminobutyric acid type A (GABAA) antagonist leaves place fields in the hippocampus intact but affects navigation performance in rats (50). Similar to presented findings in MEC, this differentiated effect indicates that the rhythmic input from MSA may not determine spatial representations in the hippocampus but may still support navigation and memory function. It remains to be tested whether the MSA input to MEC serves similar functions for memory and navigation as is indicated for the hippocampus.

PV+ cells in MSA send direct projections to the parasubiculum (PaS) (51), which again project to MEC layer II. We can therefore not exclude that changes in MEC LFP during optogenetic stimulation in MSA are indirectly driven from inputs from the PaS. On the other hand, the short delay in single-cell response times, primarily by inhibitory neurons in MEC, suggests that the direct GABAergic projection from MSA to MEC is, in large part, responsible for coordinating the local inhibitory network and pacing the LFP. While the role of PaS inputs to LII cells in MEC remains unknown, the dissociation between LFP and grid cell firing in presented study is clear.

Although ChR2 expression and therefore optogenetic activation was restricted to PV+ cells in MSA, the activation responses that we observed in the MEC (Fig. 2G) may not be limited to terminals of PV+ cells in MEC. The MSA projection includes glutamatergic projections that terminate on both excitatory and inhibitory neurons and provides speed-correlated input, mainly to pyramidal cells in superficial layers of MEC, but has not been found to terminate on stellate cells (29, 52). The number of neurons in MEC responding with fast glutamatergic activation upon general stimulation of MSA fibers is significantly less than neurons responding with fast inhibition (29). It is therefore more likely that the activation pattern that we observed in grid cells were due to disinhibition than direct activation by glutamatergic inputs. Another possibility is activation of the cholinergic component of the MSA projection. The cholinergic projection is found to synapse onto all functionally distinguishable cell types in MEC, and the postsynaptic responses may vary according to the composition of cholinergic receptors on the target neurons (53). Furthermore, the GABAergic inhibition onto MEC interneurons might influence the incoming septal glutamatergic or cholinergic projections into MEC. However, when considering the delay of about 8 ms from stimulation of PV+ cells in MSA until unit response in MEC, it seems unlikely that the putative excitatory neurons were directly activated by the stimulation. We observed a consistent and robust temporal relation between inhibitory input and grid cell activation with inhibition of putative inhibitory units always preceding putative disinhibition. Moreover, we observed no change in running speed, which would be expected if the glutamatergic pathway was activated (44). Last, when solely stimulating PV+ cell terminals in MEC, we observed inhibition of NS units and disinhibition of grid cells similar to when we stimulated cell bodies in MSA. Moreover, if antidromic signals activated PV+ cell bodies in MSA, we would perhaps expect to see stronger effect in LFP theta during MEC stimulations, which we did not observe. We believe that the difference in effects on LFP theta power between MSA and MEC stimulation is likely due to a far larger number of GABAergic terminals activated when stimulating in MSA compared to local stimulation of terminals below the fiber tip in MEC. Thus, it is expected that we observe a stronger effect on the LFP during MSA stimulation. On the basis of current observations, we therefore believe that the main effect on grid cells is due to disinhibition by activation of PV+ projection neurons in MSA.

The increase in out-of-field spikes during 11-Hz stimulation opens the question of whether it would be possible to disrupt grid cells with a different stimulation paradigm. We observed no further increase in out-of-field spikes when stimulation pulses were increased from 11 to 30 Hz and with no further reduction in spatial information. However, the decrease in spatial selectivity during disinhibition indicates that the grid fields could be disrupted if the stimulation frequency was increased to such an extent that the interpulse interval was on the same order as the response duration, hence inducing a constant disinhibition. Whether the grid cell network would be able to maintain its activity pattern also during higher stimulation frequencies remains to be determined and may shine light on the remarkable stability observed in grid cells.

Although optogenetic stimulation initiated a brief increase in grid cell activity causing increased out-of-field spiking, the fields remained stable when examining time-averaged rate maps. We do however acknowledge that negative effects must be supported by sufficient power to reliably assess the lack of change between baseline and stimulation. However, estimating power in a mixed measures design fitted by LMM is not trivial. Typically, in a repeated measures design, sample size can be reduced to roughly half (54) compared with independent measures design. Therefore, we chose to estimate an upper bound of detectable effect sizes by disregarding that many neurons were repeated across baseline and stimulation sessions. To compute this upper bound, we performed a nonparametric estimate of the distributions per animal using a Gaussian kernel density estimate (KDE). We then resampled these distributions and shifted the mean difference per animal to a set effect size before calculating the P value. Repeating this 100 times per effect size, we assessed the power as the probability of P < 0.05, i.e., rejecting the null hypothesis while the alternative hypothesis was true. By shifting the effect size across a range of values, we were able to obtain the upper bound of detectable effect sizes at 80% power as average rate = 2.7 Hz, max rate = 6.1 Hz, gridness = 0.14, information rate = 0.23 bit/s, and information specificity = 0.11 bit/spike. Given these upper bound values, we believe that presented analysis would detect significant changes in effect sizes and thus support the notion that there are no significant spatial changes in grid cells’ patterns during stimulation.

The observed disruption of grid cell phase precession during optogenetic stimulation with seemingly no effect on spatial coding contradicts the core of OI models where phase precession is used for path integration. Phase precession is suggested as a mechanism to provide the animal with information about direction of movement, which is needed to constantly update its position (55). Moreover, it has been shown that the animal’s position on a linear track can be decoded using phase precession (56), which is maintained in two-dimensional (2D) environments (13). During the 10-min stimulation, grid cell activity was locked to the stimulation frequency, but we observed only minor temporal decay in spatial accuracy and in the spatial information of grid maps. This indicates that during stimulation, some form of mechanism for localization other than OI must be at play for the grid map to be constantly updated. This could perhaps be provided by local cues or other path integration strategies independent of phase precession. For example, we found that speed coding in grid cells was intact during stimulation; thus, the speed signal necessary to update a path integrator may, however, still be provided. Phase precession is also suggested to be important for sequence coding and memory; thus, some disruption in spatial learning that we cannot observe with presented experiments might be addressed in future studies.

It is well established that frequency and amplitude of the theta rhythm increase with running speed (57, 58), suggesting that theta oscillations play a role in speed representation. Passive transportation of an animal does not reduce overall theta rhythm but eliminates the linear velocity modulation of theta, meaning that theta frequency does not increase with increasing running speed (35). In baseline recordings the strong and positive correlation between running speed and theta frequency saturates at a modest running speed and was disrupted when we paced PV+ neurons in MSA with optogenetic stimulation. Firing rates of neurons in MEC were still modulated by speed during the same stimulation, which corresponds well with a recent finding stating that input from the MSA does not control the speed coding by firing rate of entorhinal speed cells (59). This suggests that grid cell spatial firing pattern is generated independently of the theta-speed relationship. Furthermore, the passive movement effects on grid cells are more likely due to the absence of a self-motion velocity signal, as suggested by Terrazas et al. (2).

Typically, CAN models are associated with path integration arising from velocity input where direction and speed are provided into the network. Although these models can sustain oscillatory input (39), they do not normally account for temporal activity in grid cells such as phase precession, with the exception of the work by Navratilova and colleagues (60). However, whether CAN models and recent normative models (41) would be able to sustain a stable spatial pattern during interventions similar to those presented here is yet to be determined.

There is evidence that grid cells rely on path integration, where anchoring to environmental cues is seen as an error correction strategy [e.g., (61)]. Moreover, several lines of evidence indicate that grid cells use environmental cues to stabilize the grid map (62, 61), in the absence of path information. The observed preservation of grid maps in presented data could therefore be a result of corrections by use of cues in the recording environment during stimulation. According to both the OI theory and the hybrid version, the experimentally observed pattern emerges through path integration. Therefore, we believe that we presented conclusions still hold—that the presented results contradict the established OI theory. Specifically, path integration cannot be driven by phase precession, and the emergence of the experimentally observed grid pattern must either rely on some other velocity signal or not be dependent on path integration at all. However, how theta oscillations affect the emergence of grid patterns when the animal is exploring a novel environment still remains to be explored.

Recent theoretical work (63, 64) links the temporal relation between MSA pace-making theta and the resonance and rebound spiking found in stellate cells, to underlie spatial representation, phase precession, and theta cycle skipping. Together, these models suggest that theta oscillations are closely related to the grid cell spatial representations. This is in contradiction to the experimental findings reported here, and it will therefore be interesting to see whether such models can accommodate the dissociation between temporal and spatial grid cell properties.



We used male Long-Evans PvalbCre knock-in rats (3 to 8 months old, 350 to 550 g at surgery) (32). An internal ribosomal entry site–Cre recombinase sequence cassette was inserted into the Pvalb locus using CRISPR-Cas9 gene editing. Whole-genome sequencing of the strain indicated that the donor was correctly inserted into the desired locus but also identified a random insertion of the donor cassette into an intergenic region on chromosome 1 (chr1: 265, 469, 284-265, 469, and 502). This random insertion is unlikely to cause nonspecific expression of Cre, given the lack of transcription regulators in the region. This is supported by immunohistochemistry verification, which demonstrates high coexpression of Pvalb with Cre (32). After surgeries, the animals were housed individually in GR1800 (Tecniplast, Buguggiate, Italy) individually ventilated cages, in a temperature- and humidity-controlled vivarium. A 12-hour light/12-hour dark schedule was maintained, and testing occurred in the light phase. The rats were food-deprived 18 to 24 hours before training and kept at 85 to 90% of free-feeding body weight, with water available ad libitum. Experiments were performed 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.

Surgical procedures

Surgeries were performed in an aseptic environment; rats were anesthetized with isoflurane mixed with air (5% induction, 1.5 to 2% for maintenance) and immobilized in a stereotaxic frame (World Precision Instruments Ltd., Hertfordshire, UK). They were given subcutaneous injections of buprenorphine (0.04 mg/kg) and local subcutaneous injections of bupivacaine/adrenaline (Marcain adrenaline, 13.2 mg/kg) in the scalp before surgery began, which was cleaned and shaved with ethanol and chlorhexidine. The hind paw withdrawal reflex was used to assess the depth of anesthesia, together with continuous monitoring of heart rate and core temperature throughout the operation by using a MouseSTAT system with feedback mechanism to a heating pad (Kent Scientific, CT, USA). At the end of the surgery, animals were given a subcutaneous injection of carprofen (5 mg/kg), and the edge of the wound was cleaned followed by local anesthetic ointment lidocaine. This was repeated for 3 days after surgery.

Injection of opsin clone–carrying virus and electrode implantation. Craniotomies were made above the MSA and bilaterally above the MEC using a handheld Perfecta 300 dental drill (W&H Nordic, Täby, Sweden). Tetrode implantation and injections were performed in one surgical session. A Hamilton syringe (Harvard Apparatus, MA, USA) was used to inject a viral vector, AAV5-Ef1a-DIO-hChR2(H134R)-EYFP (UNC) into MSA. Injections were done in two locations, with two depths at each location (6.5 and 7.5 ± 0.2 mm measured from the surface of the skull). The locations were at 0.6 mm and 1.0 mm anterior of bregma and at the midline, with a total volume of 300 nl at each location. The injections were made in a stepwise manner with a rate of 50 nl/min, and the needle was left at the injection site for another 5 min before retraction. Tetrodes/optrodes were implanted above MEC at 0.4 ± 0.1 mm anterior of the transverse sinus and 4.5 ± 0.1 mm from the midline. Electrodes were held by 16-channel microdrives (Axona, St. Albans, UK) and referenced by ground using jeweler’s screws fixed to the skull. Implants containing either optic fibers, tetrodes, or both were anchored in the skull using jeweler’s screws and dental cement. To avoid light from the laser leaking through the implant, it was sealed with black nail polish.

Extracellular recordings and processing

Recording setup. Axona microdrives were connected to Intan RHD2132 headstages with custom adapters (from the Axona connector of the microdrive to the Omnetics connector of the headstage). The headstages were connected to the Open Ephys acquisition board (65) via thin Serial Peripheral Interface (SPI) cables. The recordings were acquired using the Open Ephys GUI for saving the data and visualizing LFP, high-pass–filtered signals, and spikes in real time. In addition, we used a Point Grey Flea3 camera mounted on the ceiling to track the animal’s position. The tracking was monitored in Open Ephys using the Tracking plugin (66), which uses Bonsai for image analysis and Open Ephys for synchronization, storage, and visualization of the trajectories. Optogenetic stimulation was performed using a 473-nm Blue Diode-pumped solid-state laser (100 mW; Shanghai laser, China) with power between 20 and 50 mW. The laser was triggered using the Pulse Pal 2 stimulator controlled by the Open Ephys GUI. We used train of pulses from the Pulse Pal stimulator to the laser controller. The train consisted of pulses of 5 ms and interpulse intervals of 85 ms (for 11-Hz stimulation) or 28 ms (for 30-Hz stimulation).

Recording procedure. The animals were habituated to the recording environment during daily training sessions starting a week before surgery. The recording environment consisted of a 1 m–by–1 m black box containing a white A4 sheet on one of the walls, serving as a local cue. To motivate animals to actively explore the environment, chocolate crumbles (Weetos) were randomly thrown into the arena for the entire recording session. Electrophysiological recordings were initiated 3 weeks after surgery to ensure proper virus expression. During recording sessions where we searched for grid cells, animals ran in the arena for sessions lasting from 10 to 20 min (depending on how fast they could cover the entire environment). Between sessions, if there were no grid cells recorded, the tetrodes were adjusted downward in steps of 50 μm, and we allowed the animals to rest in their home cage for 20 to 30 min.

When grid cells could be detected, optogenetic experiments were conducted and consisted of a block of four recording sessions. This started with a baseline session with no laser stimulation (Baseline I), followed by an 11-Hz stimulation session where the laser stimulation lasted for 10 min. Optogenetic stimulation started 1 min into the recording session. The next recording session was another baseline session (Baseline II), followed by a last 30-Hz stimulation session. All sessions lasted 10 to 20 min and were separated by a 10- to 20-min break where the animals rested in the home cage. The electrodes were not moved between these sessions to ensure that the same cells could be recorded across baseline and stimulation sessions. However, after one block of optogenetic experiments were conducted, we again lowered the electrodes searching for new grid cells at a more ventral position. Two of the animals had optic fibers in MEC in addition to MSA. Optogenetic stimulation in MEC was conducted as described above, but only for a single hemisphere at a time and only using 11-Hz laser pulses.

Spike sorting. Recordings were spike-sorted automatically using Kilosort2 (, run with the SpikeInterface framework (67). Before spike sorting, the signals were filtered with a 300- to 6000-Hz band-pass filter, and common median reference was applied (separately for each microdrive) to reduce common mode noise. The spike sorting output was curated using Phy (68), for manually rejecting noisy clusters, merging similar clusters, and splitting and cleaning clusters based on their principal components features. After manual curation, the data were saved to Exdir format (69) for further processing.

Automatic tracking of units over different experimental sessions. The tracking of the same units over different experimental sessions was done with an automatic approach. For each unit of each session, we computed a normalized similarity score with units recorded on separate sessions. Defining TiNxT as the average waveform (or template) recorded on a tetrode of a unit i (with N channels and T samples), the similarity score to another template from a different session (on the same tetrode) was computed assij=nt(TintTjnt)2max(Ti,Tj)N

Then, we built a graph with each unit of each session as a node and the similarity scores as edges (note that units from the same recording session could not have an edge). Last, the graph was interrogated to extract the minimum paths. Nodes along these paths were considered to be the same unit over time. We used a threshold for the minimum acceptable similarity score to build an edge (0.035) and for how far apart in time units could be considered to be the same (15 days). Moreover, we ensured that tracked units were found at the same depth: When a tetrode depth was adjusted, paths along the graph were stopped. The output of the automatic tracking algorithm was visually inspected to assess the goodness of the tracking.

Histology and identification of recording position

At the end of the experiments, animals were deeply anesthetized by an intraperitoneal injection of pentobarbital sodium (50 mg/kg) and intracardially perfused with 0.9% NaCl, followed by 4% paraformaldehyde (PFA) in 1× phosphate-buffered saline (PBS). The brains were dissected out and postfixed for 24 hours in 4% PFA. Before sectioning, the brains were cryoprotected in 30% sucrose in 1× PBS for 3 days and then cut into three parts: one anterior part for coronal sectioning and two posterior parts for sagittal sectioning. Forty-micrometer coronal (for locating MSA) and sagittal (for locating MEC) sections were then cut with a cryostat (Leica Biosystems, Buffalo Grove, USA). Staining procedures were performed on free-floating sections under constant agitation unless mentioned otherwise. An antibody for green fluorescent protein (GFP) was used to visualize the virus expression and anti-PV for visualizing PV+ cells.

First, the sections were rinsed times in 1× PBS and blocked with 2% goat serum and 0.3% Triton X-100 in 1× PBS for 1 hour at room temperature. The sections were then incubated with chicken anti-GFP (1:1000; Invitrogen, no. A-10262) and rabbit anti-PV (1:2000; Swant, no. PV27) in blocking solution overnight at room temperature. On the following day, the sections were rinsed with 0.3% Triton X-100 in 1× PBS and incubated for 2 hours in 1× PBS with secondary antibodies goat anti-chicken Alexa Fluor 488 (1:1000; Invitrogen, no. A-11039) and donkey anti-rabbit Alexa Fluor 647 (1:1000; Invitrogen, no. ab150075). After incubation, sections were rinsed in 1× PBS and incubated with NeuroTrace 430/455, blue fluorescent Nissl stain (1:100; Thermo Fisher Scientific, no. N21479), for 30 min. Sections were lastly rinsed in 1× PBS, mounted on SuperFrost Plus (Thermo Fisher Scientific, no. J7800AMNT) and coverslipped with FluorSave reagent (Millipore, no. 345789-20ML). Images were acquired with an Andor Dragonfly spinning disc microscope using the Fusion software, with a Zyla 5.5 scientific complementary metal-oxide semiconductor camera covering 2048 × 2048 pixels. Overview images of sections were acquired using a 20× objective [numerical aperture (NA) 0.75]. High-magnification images of PV+ somas and virus expression were acquired within the MSA and MEC through a 60× objective (NA 1.4).

The sagittal sections containing tracks from the tetrodes were mounted directly on SuperFrost Ultra Plus slides (Thermo Fisher Scientific, no. J3800AMNZ) after sectioning and counterstained with Nissl for visualizing cell bodies using Cresyl violet. Sections were then dehydrated in ethanol baths (70 to 96%) and xylene before being coverslipped with Entellan. Tetrode tracks were identified, measured, and photographed through an Axioplan 2 microscope (Carl Zeiss, Oberkochen, Germany). High-resolution images were then stitched together using the MosaiX extension in the AxioVision software (Carl Zeiss, Oberkochen, Germany).

Quantification and statistical analysis

All tests regarding LFP were nonparametric because of distributions failing normality tests (scipy.stats.normaltest). For LFP analysis, we used the Wilcoxon signed-rank test (scipy.stats.wilcoxon) with data paired on the animal level. Cell data were tested using LMMs to account for (i) non-independence (animal and unit groups), (ii) the mixed-design nature of the dataset (some neurons were repeatedly measured while others were only measured in one of two compared sessions), and (iii) imbalance (some animals had higher cell count than others). The statsmodels (70) package was used to run the LMMs. In particular, we used the statsmodels.regression.mixed_linear_model.MixedLM.from_formula function with formulas ´variable ~ label´, where variable denotes the measured variable of interest, e.g., gridness, and the label denotes two sessions, e.g., Baseline I and 11 Hz. The groups were set to groups = ´entity´ where entity denotes animal ID, and a variance component specifying the unit ID as a categorical variable vc_formula = {´1´: ´0 + C(unit_id)´}, where the zero suppresses an intercept. Last, the model had random intercepts and random slopes, where the former is used by default in MixedLM and the latter is specified by re_formula = ´label´, allowing the groups to have random initial values and vary randomly, respectively, from session to session.

In some cases, we ran into convergence issues with the statsmodels optimizers and were unable to fit the above-described model. In these cases, we altered the model and indicated these alterations by adding a § or §§ to the P values. For §, we removed the variance component such that only animals were considered as groups, and this removes the repeated part of the model design and lead to overestimation of the P value. For §§, we both removed the variance component and the slopes model given by the random effects formula. In this case, we underestimate the P value, as the model does not allow for intra-animal variations.

Because of imbalanced neural data, where the number of recorded neurons from each animal was different (see table S1), presented averages from neural data were weighted. First, a dataset was defined as D={xi=1,j=1,xi=2,j=1,,xi=N1,j=1,xi=1,j=2,xi=2,j=2,,xi=N2,j=2,,xi=NN,j=N}, where neural activity xi, j (for example, gridness or spike rate) was recorded from neuron i and animal j, Nj denotes the number of neurons from animal j, and N denotes the number of animals. The weighted average was given by x¯w=xw, wherexw=j=1Ni=1Njwi,jxi,jj=1Ni=1Njwi,jand the weighted variance defined asσw2=(xx¯w)2wThe weight was found by wj=1Nj, where Nj denotes the number of neurons recorded from animal j.

All analyses were performed in Python. We used software developed at Centre for Integrative Neuroplasticity (CINPLA). In particular, data management was performed with Expipe (71) using Exdir (69) to store raw data. Furthermore, analysis code can be found for spatial maps (, for speed cell analysis (, for phase precession (, and for spike waveform ( All paper-specific analyses are also documented at with reference to Docker image with all software installed.

Spiking rate maps. Spiking rate maps were produced by dividing the arena into equally sized bins. For each bin, the number of spikes and the time spent within the bin were counted to produce a spike map and an occupancy map, respectively, with a bin size of 2 cm by 2 cm. The spike map and the occupancy map were smoothed individually by using the convolution of a 2D Gaussian kernel provided by Astropy (72). Each bin in the smoothed spike map was divided by the corresponding bin in the smoothed occupancy map to reduce edge distortions. Each smoothed rate map was produced with SD σ = 4 cm.

Autocorrelogram. Spatial autocorrelations used for visualization and subsequent analyses were created by correlating the smoothed rate map with itself using the scipy.signal.fftconvolve function from the SciPy package with the mode parameter set to mode = ´full´.

Gridness score. For each rate map, a gridness score was calculated by first masking the central peak with a disc having a radius of half the distance to the closest peak found among six surrounding peaks. Furthermore, the area of the autocorrelogram outside a disc centered on the center peak was masked with a radius 3/2 times the distance to the outermost of the six closest peaks. Next, the masked autocorrelogram was rotated by increments of 30° up to 150°. For each rotation, we calculated the Pearson product moment correlation coefficient against the originally masked autocorrelogram. Last, the gridness score was calculated by taking the lowest coefficient found with rotations 60° and 120° and subtracting the highest coefficient found with rotations 30°, 90°, and 150°.

Shuffling. Null distributions of gridness and spatial information were generated by a randomized spike train for each registration using the function spike_train_surrogates.dither_spike_train from the Elephant electrophysiology analysis toolkit ( For each session, n = 1000 randomized spike trains were generated using a shift of 30 s with edges = True.

Spatial shift. To compute the shift of fields of the rate map, we computed the cross-correlation created by correlating the smoothed rate maps using the scipy.signal.fftconvolve function from the SciPy package with the mode parameter set to mode = ´full´. The spatial shift was then given by the center of mass of the cross-correlogram.

Identification of firing fields. To estimate the firing within and outside fields, we first identified the individual fields in the rate map. Following the protocol of (73), we first identified a global field radius as 0.7 times half the distance from the center peak to the closest peak in the autocorrelogram. Next, we identified all the peaks in the rate map before excluding the lower of any two peaks within a distance shorter than the global field radius.

To define the extent of each field, we first used morphological dilation to enlarge bright regions and shrink dark regions using skimage.morphology.dilation from the package scikit-image (74). To extract the boundaries of fields in the dilated image, we further used a Laplace filter using Gaussian second derivatives. The Laplacian of Gaussian was calculated using the ndimage.gaussian_laplace with sigma = 2.5 function from the SciPy Python package (75).

To separate and label the remaining regions, we used the ndimage.label function. Any regions with an area less than nine bins were excluded. The areas were then sorted on the basis of the mean firing rate in each area.

Fields were defined to be any labeled region found that corresponded with a non-excluded peak from the protocol of (73). The fields were used in subsequent analyses to identify in- and out-field spikes.

Spatial information and specificity. Spatial information is an estimate of to which degree the animal’s position can be predicted on the basis of the firing of the cell, given in bits per second (14). Spatial information was found using the following equationspatial information=ipiλiλ log2λiλwhere pi is the probability of the animal being in bin i, given by the occupancy map divided by the total session time; λi is the firing rate in bin i; and λ is the mean firing rate. Spatial information specificity was calculated from spatial information divided by average rate.

Speed score. The speed score was computed according to (43), briefly described below. To compute the correlation between running speed and the firing rate of neurons, we first calculated the instantaneous speed and interpolated linearly to match the sampling frequency of the firing rate. The firing rate was computed by convolving a Gaussian kernel with the spike times using statistics.instantaneous_rate and kernels.GaussianKernel from the Elephant package (76). Last, the speed score was computed as the Pearson correlation coefficient between the instantaneous speed and the firing rate for speed between 0.02 and 0.5 m/s.

Since we had strong disinhibitory responses both in grid cells and NSi cells, a window of 5 to 11 ms and 3 to 11 ms, respectively, was masked out when computing speed scores (larger window for NSi cells to account for the initial inhibition).

Waveform analysis. To separate between NS and BS units, we calculated the time from trough to peak and the time to cross the half-width amplitude of the largest-amplitude trough from the mean waveform. In peak-to-trough calculations, the sampling period of each spike was increased 200-fold by cubic interpolation to get an accurate measure of peak times. The half-width crossing time was refined by a linear interpolation between crossings of the constant line of half amplitude. All interpolations were done with scipy.interpolate.interp1d. Last, we separated clusters with scipy.cluster.vq.kmeans on the baseline and stimulation sessions.

LFP spectrum analysis. In all LFP-related analyses, we removed transients most likely due to chewing artifacts by thresholding the signal and setting the signal to zero (LFP >2 mV → 0); for analysis where spikes were compared with LFP, we also removed spikes that fell within the time of these artifacts.

The time-frequency spectrum was calculated by means of a continuous wavelet transform on the z-scored LFP signal using a Morlet wavelet with nondimensional frequency ω0 = 80 (77) and the PyCWT library ( Then, a weighted histogram of speed v with bin size 0.02 m/s was calculated in the range v ∈ [0.02,1] m/s. The weights were either the mean power of the wavelet spectrogram over frequencies f ∈ [4,12] or the frequency at the maximum power within the same frequency range. The PSD was calculated on LFP signal using the Welch method given by mlab.psd from the Matplotlib package (78).

The confidence intervals presented in Fig. 2 (A to C) were calculated by bootstrapping each data point from each recording at the 95% confidence level. To compute the spectrum as obtained from the time frequency representation versus running speed v, we first calculated the instantaneous speed and interpolated linearly to match the sampling frequency of the LFP signal at 1000 Hz.

Spike phase estimation. To estimate the phase preference of spikes relative to LFP, we filtered the LFP signal in a preferred frequency band and performed with a Butterworth filter of order 3 using scipy.signal.butter_bandpass and scipy.signal.filtfilt. The phase was acquired from the Hilbert transform scipy.signal.hilbert of the filtered LFP, and the angle was obtained using numpy.angle. From each recorded neuron, a phase was thus found for each spike, and to estimate its angular preference, we performed a KDE using the von Mises model given byf(ω)=12πI0(κ)exp (κ cos (ωμ))where I0 is a modified Bessel function of order 1 using scipy.special.I0 with κ = 100, f(ω + 2π) = f(ω), and02πf(ω)dω=1

From the probability density function generated with the von Mises KDE, the mean resultant vector length R=R2+C2/n with R = ∑n cos (ai) and C = ∑n sin (ai) for angular data a of size n (computed with pycircstat.resultant_vector_length) and mean angle (pycircstat.mean) was used to compute the spike-phase score and direction, respectively.

Spike LFP coherence. To evaluate spiking coherence to LFP (Pxy), we used the sta.spike_field_coherence from Elephant (76), which measures the coherence between a binned spike train and the LFP using scipy.signal.coherence. Signal peaks were extracted using scipy.signal.find_peaks, and the energy was computed by (Σ (Pxyf), where Δf is the corresponding frequency resolution.

Phase precession. Phase precession was computed according to (13, 36) where the normalized distance traveled within a field and the corresponding phase it fired relative to the band-pass–filtered LFP (6 to 10 Hz). To measure the linear circular correlation, we translated a MATLAB script obtained from Kempter et al. (36) and translated to Python using SciPy, NumPy, and the PyCircStat (79) package. The phase precession quantification for Python is available at the CINPLA GitHub site

The LFP signal was band-pass–filtered in the theta frequency range (6 to 10 Hz). The instantaneous theta phase was found using the Hilbert transform of the filtered signal. Each spike was assigned an instantaneous theta phase by interpolation using the scipy.interpolate.interp1d. Phase precession was quantified by two measures: First, assuming that the association between phase and position can be described by a linear model of the form ϕ̂=2πax+ϕ̂0 (where a is the phase, x is the distance traveled through the firing field, ϕ0 is the phase offset, and ϕ̂ is the approximation of the true phase ϕ), we obtained the slope by maximizing the goodness of fit a=arg maxa*R(a*)R=(1nk=1ncos (ϕk2πaxk))2+(1nk=1nsin (ϕk2πaxk))2and computing the phase offset withϕ̂0=arctan*(ksin (ϕk2πaxk)kcos (ϕk2πaxk))where arctan * is the quadrant-specific inverse of the tangent (numpy.arctan2). Second, a circular-linear correlation was computed byr=k=1nsin (ϕkϕ¯)sin (θkθ¯)k=1sin (ϕkϕ¯)2k=1sin (θkθ¯)2with ϕ¯=arctan*ksin (ϕk)kcos (ϕk) and θ¯=arctan*ksin (θk)hcos (θh), where θj = 2πaxj(mod 2π).

For large sample sizes, P values can be derived (36, 80) withp=1erf(z2)where erf is the error function and z=rλ20λ02/λ22 with λij=n1k=1nsin (ϕkϕ¯)isin (θkθ¯)j. When assessing whether a neuron phase precess r < 0 or recess r > 0 two parameters were used, we required a P value less than P < 0.01 and a goodness of fit R > 0.1.

Stimulus response probability and cell-type characterization. To compute stimulus response (PSTH), we first aligned spikes to each stimulus onset and cut out spikes between −10 and 50 ms with respect to the stimulus onset. A kernel density estimation was then performed with the stats.gaussian_kde with sigma = 0.1 ms to get the probability density of spiking. Then, the 99th percentile of the density before the stimulus was used to determine whether the stimulus response was significant. Response latency was further obtained as the time-to-peak or time-to-trough response for excitation and inhibition, respectively. We considered responses as inhibitory only when the time to trough of the response preceded the time to peak.

Units that were characterized as NSi once by the significance of inhibitory response were labeled as NSi regardless of response profile in other sessions. Full width at half maximum was computed from the maximum peak within an interval between 5 and 30 ms of the PSTH.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We wish to thank L. Frank for donating the PV-Cre rat line, S. Grødem and T. Klungervik for technical assistance, and S. Leutgeb for valuable input during the project, both on the experimental setup and during analysis of results. We would also like to thank M. M. Frey, M. B. Wigestrand, A. Moen, and M. B. Røe for helping out in the outset of the project and for being involved in pilot experiments. Funding: This work was funded by the Research Council of Norway (grant nos. 217920 and 248828 to M.F. and grant no. 231248 to T.H.) and the University of Oslo’s Strategic Research Initiative CINPLA. Author contributions: M.E.L., M.F., and T.H. conceived the project, performed pilots, and determined the experimental protocol. A.C.C. performed electrophysiological experiments. A.C.C. and K.K.L. performed immunohistochemical procedures. M.E.L. and A.P.B. developed the analysis platform. M.E.L. and A.C.C. performed the analysis. J.Y. developed the rat line. All authors contributed to writing the manuscript. Competing interests: All authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The data are available at, and all analyses are available at where you also will find a link to a Docker container with all the necessary software installed.

Stay Connected to Science Advances

Navigate This Article