Intrinsic theta oscillation in the attractor network of grid cells

Summary Both grid-like firing fields and theta oscillation are hallmarks of grid cells in the mammalian brain. While bump attractor dynamics have generally been recognized as the substrate for grid firing fields, how theta oscillation arises and interacts with persistent activity in a cortical circuit remains obscure. Here, we report that the theta oscillation intrinsically emerges in a continuous attractor network composed of principal neurons and interneurons. Periodic bump attractors and the theta rhythm stably coexist in both cell types due to the division of labor among interneurons via structured synaptic connectivity between principal cells and interneurons. The slow dynamics of NMDAR-mediated synaptic currents support the persistency of bump attractors and restrict the oscillation frequency in the theta band. The spikes of neurons within bump attractors are phase locked to a proxy of local field potential. The current work provides a network-level mechanism that orchestrates the bump attractor dynamics and theta rhythmicity.


Theta rhythm emerges intrinsically in a continuous attractor network (CAN) of grid cells
Labor division among interneurons ensures both the persistent and rhythmic activity Dominance and slow dynamics of NMDA receptors are required at recurrent synapses Neural firings are phase locked to the local field potential (LFP) and external rhythmic input

INTRODUCTION
Grid cells in the medial entorhinal cortex (mEC) constitute a major cell type in the brain's navigation system; their hexagonally distributed firing fields provide a spatial metric. 1,2 Their temporal characteristic, typically theta oscillation (4)(5)(6)(7)(8)(9)(10)(11)(12), is another important aspect of functioning. 3 This theta rhythm is present in spike trains of individual grid cells and the local field potential (LFP). 3,4 Functionally, theta oscillation is strongly associated with the movement of animals, with its frequency possibly reflecting the speed or acceleration information. 5,6 Sensory information received by the hippocampal-entorhinal circuit can be integrated within temporal windows provided by theta oscillation. 7 The uncertainty about the animal's position can be significantly reduced by taking into account the timing of neural spiking relative to the LFP theta phase. 8 Indeed, the hexagonal lattice of firing fields and theta rhythm are closely correlated. For instance, cells with the same grid scale and orientation constitute a grid module and exhibit more coherent oscillations than intermodular ones; 9 there is a significant negative correlation between the theta frequency and grid spacing of individual cells. [9][10][11] When a rat moves through one firing field of a grid cell, either phase locking or phase precession appears between its spike train and the LFP, and the degree of precession depends heavily on the rat's position. 3,[12][13][14] Thus, exploring both the features within the same framework is essential to unraveling the dynamics and function of grid cells.
Two major kinds of computational models have been proposed to account for spatiotemporal firing features of grid cells: the continuous attractor network (CAN) model and oscillatory interference (OI) model, which are at the network and single-neuron levels, respectively. With the CAN model, the hexagonally distributed activation bumps on a neural sheet with toroidal topology can be topographically mapped to the two-dimensional (2D) space where the animal moves. [15][16][17] Moreover, given the specific boundary condition of the neural connectivity profile, the CAN model admitting a single bump attractor can also account for the grid firing pattern. 18 With the OI model, the interference between the theta oscillation with constant frequency and that with the frequency varying with the animal's velocity underlies the periodic firing fields of single cells. 5,11,19 Although many experimental results can be understood in the framework of either the CAN model or the OI model, few studies tried to combine them, 20,21 partly due to their radical differences in assumption and dynamics. Notably, both the assumptions and predictions of the CAN model and its extensions have been extensively validated in experiments. 4,12,17,22,23 Two mechanisms underlying theta oscillation in mEC have been suggested. The theta rhythm could result from the input from the medial septum (MS). Interneurons in superficial layers of the mEC receive strong

Sustained firing activity in the continuous attractor network model
The superficial layers of the mEC contain large numbers of excitatory principal cells and inhibitory interneurons. In previous 2D CAN models, cells are uniformly distributed on a sheet and recurrently connected, 16,35 while their firing activities are approximated by gating variables of specific synaptic receptors such as NMDA receptors (NMDARs). As limited mutual excitation between principal cells was found experimentally, 35,36 there is usually no direct connectivity between excitatory neurons in these models, and the network is dominated by recurrent inhibition. As neurons on opposite sides of the sheet are interconnected as neighbors, this sheet forms a torus in topology, and the synaptic strength between neurons relies on their distance on the torus. Given strong excitatory input, a hexagonal lattice of bell-shaped activation bumps can be elicited and is able to shift as a whole when the animal runs in the 2D space. Accordingly, individual cells display a hexagonal lattice of firing fields. It is the specific synaptic connectivity profiles, inhibition-dominated recurrent interaction, and strong external input that are essential to the CAN model of grid cells. As the current study aims to elucidate the roles for neurons at different angular positions in maintaining persistent and rhythmic firing, a 1D ring model incorporating the above essential features serves as a good starting point ( Figure 1A; STAR Methods). We do not anticipate that a 2D version of our model would act in a qualitatively different manner. Indeed, similar simplifications have been made to explore iScience Article the dynamics of synapses and membrane potentials of grid cells, and the conclusions with 1D models are similar to those with 2D models. 34,37 Here, each unit is described by the leaky integrate-and-fire (LIF) model, and neurons are connected in a manner analogous to that in the CAN model of gating variables 35,37 ( Figure 1B; STAR Methods). In the default simulation protocol, each cell is subjected to a constant excitatory background input I E Back (I I Back ) plus independent noise, and each interneuron also receives a constant inhibitory current from the MS (I MS,0 ). Two bump attractors can spontaneously arise and persist in the network without mutual excitation between principal cells; to initialize the bump attractors centered around 90 and 270 across trials, an additional brief current of 200 pA is applied to some principal cells for 500 ms (Figure 2A). Activated principal cells have a maximal mean firing rate of 3.6 Hz over the period of 1-8 s. The two attractors are separated by 180 and roughly identical in activity; they can be regarded as a minimal periodic structure on the ring, analogous to a hexagonal bump lattice on the torus in the 2D network. 37 The cells within each attractor discharge synchronously to various extents ( Figure 2B). Owing to the structured connectivity from principal cells to interneurons, the interneurons with similar angular locations to the principal cells firing persistently iScience Article display elevated firing ( Figure 2C), with higher firing rates and weaker synchronization than those principal cells ( Figure 2D). Together, the network is endowed with sustained activity as in 2D networks.
This persistent activity is evoked and maintained due to the strong external input and slowly decaying NMDAR-mediated currents from principal cells to interneurons ( Figure 2E). When a cluster of principal cells are first activated, they excite the interneurons at similar angular positions, and the excitation persists for a relatively long time period. At m I/E = 90 (half the angular distance between the two maxima in the profile of W I/E ; see Equation 10), these activated interneurons strongly suppress the principal cells 90 apart, which stabilizes the firings of the previously activated principal cells, giving rise to bump attractors. As a result of overall competition among principal cells, two clusters separated by 180 fire persistently, and the others fire sporadically. The number of concurrent bump attractors is mainly governed by m I/E ; more attractors can be elicited by decreasing m I/E and adjusting other parameters (Figures S1A-S1C).
On the other hand, the stability of bump attractors is affected by the characteristics of NMDAR-mediated current, e.g., its relaxation time t NMDA (Equation 6). Owing to the stochasticity of external input, the angle q PAV of the population activity vector (PAV), which is the Rayleigh vector of angular positions of all principal cells weighted by their firing rates (see STAR Methods), can drift randomly over time, but the drift is markedly repressed for sufficiently large t NMDA ( Figure 2F). The magnitude of the PAV fluctuates slightly around the mean: its coefficient of variation over 1-8 s is 0.142 G 0.022 and 0.146 G 0.012 (mean G S.D., here and thereafter) for t NMDA = 100 ms and 150 ms, respectively. The standard deviation of q PAV tends to rise over time, but the whole trace downregulates with increasing t NMDA in the range of 50-150 ms ( Figure 2G). Overall, the structured synaptic connectivity and slow NMDAR-mediated currents underlie stable attractor dynamics.

Theta rhythm in neural firing
To verify the rhythm in neural activity, we performed power spectral density (PSD) analyses of multiple timeseries signals (see STAR Methods). We first calculated the mean of t GABA -delayed GABAR-mediated synaptic currents over all principal cells (Equation 18). This quantity can be taken as a proxy of the LFP (pLFP), since it was shown that this synaptic current in the LIF network closely matches the LFP signal produced in a morphologically realistic 3D network model. 38 The pLFP obviously exhibits a rhythm ( Figure 3A), and the PSD of its temporal evolution over 1-8 s is calculated and averaged over 100 independent trials, showing a characteristic peak around 9.18 Hz ( Figure 3B). The primary peak frequency on individual trials always falls within the theta range (9.17 G 0.76 Hz) ( Figure 3C), while the primary peak power, defined as the power in a 2-Hz bin centered at the primary peak frequency divided by that in the range of 0-50 Hz, is 4.12 G 0.62 ( Figure 3D). Notably, this theta rhythm is coherent across a long time period ( Figure S2).
Owing to the limited simulation time, individual neurons may discharge a small number of spikes, resulting in a rather sparse auto-correlogram ( Figure S3A) and making it hard to determine the rhythm in their spike trains. Thus, we counted the spikes from all principal cells in time bins of 5 ms ( Figure 3E) and took this series for PSD analysis ( Figure 3F); the primary peak frequency is 9.36 G 0.81 Hz. Moreover, the sharply valued spike count indicates that activated principal cells are strongly synchronized. We also took the membrane potential of the principal cell at 90 ( Figure 3G) for PSD analysis ( Figure 3H), and the primary peak frequency is 8.69 G 0.98 Hz. Other principal cells, including those beyond the attractors, exhibit the similar rhythm in membrane potential ( Figures S3B and S3C). Similarly, the theta rhythm is manifested in the spike counts from all interneurons in bins of 5 ms and membrane potentials of individual interneurons ( Figure S4). Remarkably, the mean theta frequency is similar across different measures ( Figures 3I and S4F).
These results validate the theta rhythm in the neural activity. Of note, the firings of neurons within bump attractors are mainly localized near the troughs of the pLFP. This effect is quantified using the Rayleigh vector F of the neural firing phase relative to the pLFP (see STAR Methods). For the principal cell at 90 , its average arg(F) and |F| equal 228.9 and 0.883, respectively ( Figure 3J); the values of arg(F) are close to the experimental data, 3,13,39 while |F| is always greater in simulation. This discrepancy is possibly associated with the limited number of spikes a cell discharges in simulation and the difference between the pLFP and LFP measured experimentally. Interneurons also exhibit phase-locked firing; for the interneuron at 90 , its mean arg(F) and |F| separately equal 301.7 and 0.385 ( Figure S4G), suggesting a delayed phase and weaker modulation than the corresponding principal cell on average. The theta modulation of interneurons was also reported experimentally. 40 iScience Article Similar phenomena are observable in the case of more coexisting attractors where the primary peak frequency tends to be higher (Figures S1D-S1H). Together, the theta rhythm and phase locking can arise endogenously in the network.

Negative feedback between principal cells and interneurons underlying the theta rhythm
To unveil the mechanism by which the theta rhythm emerges, we ran simulations where initially all neurons have the same membrane potential of À60 mV and receive a constant external input without noise. In this deterministic scenario, all principal cells discharge simultaneously and their firing rate is around 5.07 Hz ( Figure 4A). Interneurons also fire synchronously at 101.4 Hz, and their inter-spike intervals vary periodically, with the same rhythm as the firing of principal cells. Thus, the theta rhythm remains despite the absence of attractor dynamics.
A comparison in dynamics of membrane potential between principal cells and interneurons reveals the formation of theta oscillation ( Figure 4B). The external input promotes the persistent depolarization of iScience Article principal cells. The synchronous firing of principal cells delivers a strong NMDAR-mediated current to interneurons, initiating their firing. Each time interneurons discharge, GABAR-mediated inhibitory currents tend to strongly repolarize principal cells, counteracting their depolarization by the external input. With the decay of NMDAR-mediated currents, the total currents are finally too small for interneurons to fire, and then principal cells depolarize to reach the firing threshold and discharge again; a new round of oscillation resumes. Therefore, the theta rhythm results from the negative feedback between principal cells and interneurons.
To further examine the underlying mechanism for theta rhythm in the presence of bump attractors, we suppressed some interneurons by artificially adding inhibitory currents. iScience Article and their firing rates increase while the bump attractors become more stable. But the theta rhythm in neural firing weakens as the involved negative feedback between excitatory and inhibitory neurons is reduced. Indeed, with increasing w, more interneurons around 0 /360 and 180 are suppressed, resulting in a monotonic decrease in theta power in the PSD of the pLFP, defined as the mean power in 4-12 Hz divided by that in 0-50 Hz ( Figure 4D). These results highlight the importance of the negative feedback in mediating theta rhythmicity.
On the other hand, if mutual inhibition between interneurons is markedly strengthened compared with the default settings, the firings of interneurons between the two attractors decrease prominently, leading to enhanced firing of principal cells within the attractors and less drift of the attractors (Figures S5A-S5C).
With increasing the inhibitory synaptic conductance, the theta rhythm first declines and then is sequentially replaced by the beta oscillation (12-30 Hz) and gamma oscillation (30-70 Hz) (Figures S5D-S5I). The gamma oscillation in this case is due to strong interactions between interneurons within bump attractors. 43,44 Collectively, the theta rhythm is mediated by the interneurons in-between bump attractors.
Experimentally, the typical values of t NMDA and t GABA are separately around 100 ms and 10 ms. 45,46 Given the structured connectivity between principal cells and interneurons in the model, the theta rhythm is maintained over a relatively wide range of t NMDA and t GABA (Figures 4E and 4F). Moreover, when the slow NMDAR-mediated currents to interneurons are increasingly replaced with the fast AMPAR-mediated currents (t AMPA = 2 ms), the theta rhythm declines fast and the gamma oscillation gradually arises ( Figure S6; see Supplemental materials and methods). Thus, the predominance of NMDARs at recurrent synapses and the long relaxation time of NMDAR-mediated current restrict the oscillation frequency in the theta range.

Orchestrated variation of the attractor stability and theta rhythmicity
We have proposed a mechanism for concurrent sustained and rhythmic activities in the network without direct connectivity between principal cells. Under the default parameter settings, principal cells within bump attractors (E in ) deliver strong excitation to interneurons at similar angular positions (I in ), which then primarily inhibit principal cells 90 apart, i.e., beyond the attractors (E out ), but E out is unable to influence E in via interneurons owing to their low activity, supporting the persistent activity ( Figure 5A). Interneurons in between the attractors (I out ) are also excited weakly by E in and then inhibit E in in return, promoting rhythmic activity. That is, different pools of interneurons serve distinct functions.
Accordingly, the self-sustained firing and theta rhythmicity can vary in an orchestrated manner under diverse conditions. The pLFP theta rhythm is enhanced with increasing t NMDA , for example, while the random drift of bump attractors weakens ( Figure 5B). The pLFP theta rhythm is also improved with increasing t GABA , whereas the random drift of bump attractors enlarges ( Figure 5C) mainly owing to the strengthened inhibition from I out to E in . More variation modes are shown in Figure S7 where different parameters are tuned. These orchestrated variations in the persistent and rhythmic activity still originate from the division of labor among interneurons.
When parameter values are adjusted in suitable regimes, the profiles of bump attractors, characterized by their height and full width at half maximum (FWHM) (see STAR Methods), do not vary much in most cases, and the primary peak frequency in the PSD of the pLFP remains in the theta range ( Figure S7). Thus, it is appropriate to use only the extent to which bump attractors drift and the theta power to separately characterize the stability of sustained activity and strength of theta rhythmicity.
As there exists limited recurrent excitation among principal cells in superficial layers of the mEC, 36 we also introduced it to our model. In this case, only NMDAR-mediated currents are considered, and the E-E connectivity strength W E/E is peaked at zero angular difference (Supplemental materials and methods). With increasing the synaptic conductance, the bump attractors become more stable and the spike times are more concentrated relative to the pLFP rhythm, while the primary peak frequency and theta power in the PSD of pLFP just slightly decrease ( Figure S8). Similar changing trends can be observed with increasing I E Back , except that the primary peak frequency rises with I E Back . Overall, the recurrent excitation strengthens the persistent firing but slightly impacts the theta rhythm in our model.
Of note, in the CAN model with structured E-E connectivity and uniform E-I and I-I connectivity, the persistent activity depends on the local excitation and global inhibition. 47  iScience Article along with the bump attractor due to the positive feedback between principal cell, 47,48 in contrast to its absence here, while fast oscillation may arise owing to the interplay between the fast AMPAR-mediated recurrent excitation and slower GABAR-mediated inhibition. However, strong rhythmic firing would diminish the persistent activity, because both the self-sustained and oscillatory activity depend on the whole population of interneurons; the self-sustained activity requires the persistent global inhibition from interneurons to stabilize the bump attractor, which is always disturbed by the rhythmic firings of interneurons. 47 Thus, distinct mechanisms are at play in these two types of networks.
Interneurons are uniformly connected in our default model setup. We further explored two cases of structured I-I connectivity: the connectivity strength W I/I between two interneurons is either a unimodal or bimodal function of their angular difference ( Figure S9 and Supplemental materials and methods). In either case, both the self-sustained activity and theta rhythm are maintained. In the unimodal case, interneurons in the same activated cluster primarily inhibit each other, which may weaken their inhibition of principal cells beyond bump attractors and enlarge the drift of bump attractors. In the bimodal case, interneurons in between the bump attractors are more suppressed by activated interneurons, which tends to reduce their inhibition of principal cells within the bump attractors, leading to more stable attractor dynamics iScience Article and weaker theta rhythmicity. It is also worth noting that g I/I GABA , maximal I-I synaptic conductance, is much less than the others, and thus changing the concrete form of W I/I does not cause large variations. Together, all these results highlight the importance of structured synaptic connectivity between excitatory and inhibitory cells in maintaining the persistent and rhythmic firing.

Modulation of neural activity by the background input
It has been reported that the MS input exerts a control over neural firing in mEC; blocking the MS input to the mEC largely disrupts both the rhythm and grid patterns of grid cells, 25,31,32 while tuning the MS input evokes distinct firing properties. 18, 39 The MS primarily delivers GABAergic projection to interneurons. 24,26 We first examined the case where the MS input to each interneuron is a negative constant I MS,0 . The stability of bump attractors weakens with decreasing |I MS,0 | ( Figure 6A), and the network fails to maintain persistent firing for |I MS,0 | < 50 pA ( Figure 6B). On the contrary, the stability is enhanced with increasing |I MS,0 | (Figure 6C), with the standard deviation of q PAV at 7 s remaining at relatively low levels for |I MS,0 | > 100 pA (Figure 6D). The theta frequency, defined as the peak frequency within the theta range, basically equals the primary peak frequency calculated above in most cases. The theta frequency drops monotonically with decreasing |I MS,0 |, whereas the theta power first rises, varies slightly for 130 > |I MS,0 | > 50 pA, and then declines. With |I MS,0 | deviating from 100 pA, the drop in theta power for |I MS,0 | < 50 pA is mainly caused by the breakdown of the attractors, while that for |I MS,0 | > 130 pA mainly results from that the primary peak frequency on some trials exceeds 12 Hz. Moreover, when the excitatory background and the inhibitory MS input are simulated with synaptic currents from Poisson spike trains, 47 the network behavior does not change markedly ( Figure S10). Together, the strength of inhibitory input from the MS markedly impacts both the persistent firing and rhythmic activity.
These results are generally consistent with the experimental observations: suppressing the input from the MS weakens the theta rhythm ( Figure 6F) and breaks the grid pattern. 31,32 Indeed, when |I MS,0 | varies from iScience Article 100 pA (default case) to 10 pA (suppressed case) in our model, both the theta power and theta frequency decrease ( Figure 6E), while the bump attractors drift remarkably and finally break down ( Figure 6D), resulting in the broken grid pattern. Overall, the deteriorating grid pattern and theta rhythm under suppressed MS input can be accounted for by a decrease in MS input strength.
Furthermore, we explored how the rhythmicity of MS input could influence the neural activity. The MS input is assumed to vary periodically-I MS (t) = I MS,0 + Acos(2pf MS t), with A and f MS denoting the amplitude and frequency of the rhythmic MS input. 18,33,34 At A = 0, there exists only one major peak in the PSD of the pLFP on individual trials and the mean primary peak frequency equals f p0 = 9.17 Hz-intrinsic oscillation frequency ( Figure 3). Figures 7A, 7B, S11, and S12 illustrate examples of neural activity for different A and f MS . At A = 40 pA, the intrinsic oscillation predominates for f MS < 4 Hz, and the mean primary peak frequency f p in the PSD of pLFP is around f p0 (Figures 7C, S11A, and S11B); by contrast, the external input rhythm dominates and f p remains around f MS for f MS > 11 Hz (Figures S11C and S11D). For 4 % f MS % 11 Hz, there emerges a resonance effect between the external input and intrinsic oscillation, and peaks in the PSD of pLFP on individual trials may appear around 0.5f MS and its integer multiples (Figures S11E-S11H); under strong resonance (6 < f MS < 9 Hz), the peak around f MS is much higher than others in the PSD (Figures 7D and S11I-S11L). Consequently, the average of primary peak frequency over 100 trials exhibits a complicated dependence on f MS ( Figure 7C). Similarly, the peak power has a local maximum around f MS = 7 Hz because of the strong resonance effect ( Figure 7E). By contrast, at A = 100 pA the network displays a different behavior, always acting as a forced oscillator; peaks always appear at multiple integers of f MS or 0.5f MS in the PSD of pLFP on individual trials, and the major peak is located around f MS ( Figure S12). iScience Article At A = 40 pA, the bump attractors have minimal random drifts and maximal sharpness when f MS is between 7 and 8 Hz, which may be called resonance frequency evoking the maximal resonance effect (Figures 7F and  7G); owing to the involved complex interactions, this frequency may be different from f p0 . For a much larger amplitude of the MS input, such as A = 100 pA, bump attractors cannot persist for most f MS except for f MS around the corresponding resonance frequency (Figures S12C, S13C, and S13D). Overall, the external rhythm disfavors the robustness of bump attractors, but the resonance effect can weaken the deterioration.
Although the phase precession was suggested to arise from the interference between an external and an endogenous theta oscillation with slightly different frequencies, 5,11,19 it is absent in our model. Instead, here the phase locking is rather robust; irrespective of f MS , spikes of both principal cells and interneurons within bump attractors are always locked to the pLFP theta rhythm although the theta rhythm is not dominant in pLFP in some cases ( Figures 7H and S13E-S13H). Moreover, when the rhythmicity of MS input is not too weak (otherwise, it hardly influences the network dynamics), spikes of both principal cells and interneurons are also locked to the input rhythm but with a varying locked phase ( Figures 7H and S13I-S13L). Thus, activated principal cells and interneurons are phase-locked to both the pLFP theta rhythm and MS input rhythm.
Notably, part of the results above agrees with the experimental observation. 39 Experimentally, the parvalbumin-positive (PV + ) interneurons in MS are stimulated at the frequency of 11 Hz or 30 Hz; the resulting peak frequency and peak power in the PSD of mEC LFP and the spiking of neurons locked to the input rhythm can be accounted for by our results for f MS R 11 Hz ( Figures 7C, 7E, and 7H). Changes in neural behavior within other f MS ranges remain to be experimentally validated.

DISCUSSION
In this study, we investigated the intrinsic theta rhythm in the CAN model with a grid-cell-like connectivity profile. Without recurrent excitation among principal cells, the connectivity between principal cells and interneurons has to be structured for generating bump attractors, which necessitates the labor division among interneurons. Interneurons within attractors receive strong excitation from principal cells within the attractors and then inhibit principal cells beyond. Given the slow dynamics of NMDAR-mediated currents, such excitation and the resulting inhibition can persist, stabilizing the bump attractors. Interneurons beyond attractors also receive sufficient excitation from principal cells within attractors and then inhibit them in return, and the resulting negative feedback allows for the theta oscillation. Together, both the structured synaptic connectivity and slow NMDAR dynamics are indispensable for the CAN model with theta rhythm.
Our work reproduces several experimental observations. First, the theta rhythms in the pLFP and spike trains of activated principal cells are nearly anti-phased, similar to the experimental reports that the spikes of grid cells are locked near troughs of the LFP. 3,13,39 Second, the declining theta rhythm and broken bump attractors due to reduced MS input agree with the experimental results when inputs from the MS are blocked. 25,31,32 Finally, given the rhythmicity of MS input at the upper end of or far outside the theta range, neural firings are locked to the MS input while bump attractors are maintained. 39 All these results support a network-level mechanism for theta rhythm in a local circuit of the mEC.
The periodic bump attractors and rhythmic activity are linked inherently in our model, different from previous studies where the theta rhythm was induced externally. 18,33,34,49 Here, the theta rhythm is endogenous and coexists stably with bump attractors, independent of specific pacemakers or rhythmic inputs. 26 The current theoretical framework could be exploited to further probe the spatiotemporal firing patterns of grid cells and their role in spatial navigation. The conclusions drawn here may be enlightening for exploring cognitive processes involving both self-sustained and rhythmic activity, such as working memory. 50,51 Two testable predictions can be made here. First, we underscore the roles for interneurons and NMDARs in generating the theta rhythm besides grid patterns. 35,52 Thus, suppression of interneurons or blockade of NMDARs in the mEC probably leads to breakdown of both the grid pattern and theta rhythm. Part of the prediction has been validated: the stability of grid pattern is reduced after the ablation of NMDARs in the retro-hippocampus region, and the LFP theta rhythm significantly weakens when NMDARs in mEC are blocked. 53 iScience Article from the MS input. Second, there exists a frequency around the intrinsic theta frequency, with which stimulating the interneurons in mEC may elicit strong theta oscillation in the LFP (e.g., the peak power in its PSD takes a local maximum) and stable bump attractors with a locally minimal drift. Experimentally, such modulation of rhythmic input might be realized by tuning the MS PV + cells via optogenetic methods. 39 The validation of these predictions would advance the understanding of grid cells.
Two assumptions are critical for the presented results: the NMDAR-mediated current and synaptic connectivity profiles. NMDARs have proved a dominating receptor mediating excitatory synaptic currents in layer II of the mEC 45 and underlying persistent activity in attractor networks. 47,55 Different from W E/I , W I/E , and W I/I adopted here, other connectivity profiles also contribute to grid patterns. 18, 33,34 The main difference lies in the angular bias in E/I and I/E connectivity profiles. This will only impact which pools of interneurons are excited to suppress principal cells beyond the bump attractors, and thus division of interneurons via structured synaptic connectivity can still be realized in these models to support neural oscillation. Overall, our assumptions have a solid experimental and theoretical basis.
Although the MS is involved in modulation of phase precession in mEC, 56 our results suggest that it may not directly arise from the interplay between the intrinsic rhythm and external input rhythm; we speculate that additional inputs related to spatial information might play a role. Possibly, real-time inputs from speed cells may bring about a shift in bump attractors over time, underpinning the phase precession. 57 The deviated ramps in the membrane potential of single neurons when a rat runs through one firing field might reflect a current tuning. 4,20,58 Theoretical works also suggested that the phase precession was related to movements of bump attractors. 34,59 Note that other kinds of oscillation exist in mEC and interact with the theta rhythm, such as the theta-nested gamma oscillation. 18 Although gamma oscillation arises in our model due to fast AMPAR-mediated currents 55 or GABAR-mediated currents, 43,44 its relationship with the theta rhythm was not explored in detail here. Different types of interneurons have distinct properties of firing and projection, affecting the spatiotemporal features of grid cells in different manners. 60,61 It is worthy to further explore the functions of diverse interneuron types. The random drift of bump attractors over time is inevitable in CAN models. It could be suppressed by sensory inputs such as those from environmental boundaries. 62 Our results suggest that it can be repressed by tuning various parameters, such as the time constant of NMDARs and applying a rhythmic input with specific frequency to interneurons. It would be interesting to systematically explore how bump attractors could be stabilized. These issues may prompt further work for more advanced models.
To conclude, in the CAN model with structured connectivity between principal cells and interneurons, the theta rhythm can naturally arise along with periodic bump attractors, and their variations are orchestrated under various conditions. A remarkable phase locking emerges between the pLFP and spike trains of both principal cells and interneurons within bump attractors. The theta rhythm gets stronger and the bump attractors become more stable when the frequency of the MS input is near the intrinsic frequency. This study sheds fresh light on the association of attractor and rhythm dynamics in grid cells.

Limitation of the study
The current model is a 1D network of grid cells. To fully interpret the spatiotemporal firing pattern of grid cells, it is still necessary to construct a 2D network model.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

AUTHOR CONTRIBUTIONS
F.L., T.W., and W.W. conceived and designed the project. Z.W. and T.W. performed the study. All authors analyzed the data. Z.W., T.W., and F.L. wrote the paper.

Article
The excitatory background inputs are constant, i.e., I E Back = 750 pA and I I Back = 325 pA. Two cases of MS GABAergic inputs to interneurons are considered: 1) I MS (t) = I MS,0 = -100 pA, which is the default setup; 2) I MS (t) = I MS,0 + Acos(2pf MS t), where A and f MS are the amplitude and frequency of the oscillatory input. The noise current is independent of each other and described by the Ornstein-Uhlenbeck process: 67 where n i (t) is the Gaussian white noise with zero mean and unit standard deviation, h Noise = 150 pA (standard deviation of the noise) and t Noise = 2 ms (time constant of the noise). Another form of input and noise currents is also modeled to exclude the dependence of main conclusions on the concrete form of external input (see Supplemental materials and methods).

Population activity vector and profile of a bump attractor
A PAV is calculated to track the central location of a bump attractor. 47 It is defined as the Rayleigh vector of angular positions of all principal cells weighted by their firing rates. There can exist two bump attractors separated by 180 with the default parameters. To track them simultaneously, the angular position of each neuron is reset by doubling (i.e., q i ´ = 2q i ), with the range converting from [0 , 360 ) to [0 , 720 ), which renders the centers of two attractors differ by 360 . Thus, the two attractors overlap when illustrated in a polar coordinate, and the PAV involving all principal cells reads where j is the imaginary unit and f i is the firing rate of neuron i over a time bin of 400 ms. The argument of P denotes the central location of the 'overlapping' bump attractors. Thus, the original angular location of either bump attractor, q n (n = 1, 2), before doubling is q n = arg ðPÞ + ðn À 1Þ$360 2 : (Equation 16) For simplicity, q 1 is always chosen to represent the PVA location, q PAV = q 1 . With the time bin sliding every 5 ms, a series of q PAV can be obtained to illustrate the drift of bump attractors over time.
Three quantities are introduced to characterize the profile of a bump attractor. Its height is defined as the maximal firing rate of neurons within the attractor, its width refers to the FWHM, and its sharpness is defined as follows: 69 where N 0 E = N E =2, q i and F i are the angular position and firing rate of neuron i, and q 0 is the center of the bump attractor. The sharpness can be regarded as an integrated feature of the height and FWHM of the attractor profile to some degree. For bump attractors of a standard Gaussian profile and a plateau profile, the sharpness is 3.0 and 1.8, respectively.
Power spectral density analysis PSD analysis is used to verify the presence of oscillation in network activity and to determine the primary frequency. The analysis is performed with Python using the signal.welch function in scipy package, involving a hamming window with the window length L, number of overlapping samples (n s ), and number of DFT points (n D ). Here, L = 256, n s = 128, and n D = 1024.
To justify the rhythm in neural activity, we took three signals for PSD analysis. First, the mean of I I/E GABA;i ðt À t GABA Þ over all principal cells is taken as a proxy of the LFP (pLFP), 38 i.e.,