Grid cell firing patterns may arise from feedback interaction between intrinsic rebound spiking and transverse traveling waves with multiple heading angles

This article presents a model using cellular resonance and rebound properties to model grid cells in medial entorhinal cortex. The model simulates the intrinsic resonance properties of single layer II stellate cells with different frequencies due to the hyperpolarization activated cation current (h current). The stellate cells generate rebound spikes after a delay interval that differs for neurons with different resonance frequency. Stellate cells drive inhibitory interneurons to cause rebound from inhibition in an alternate set of stellate cells that drive interneurons to activate the first set of cells. This allows maintenance of activity with cycle skipping of the spiking of cells that matches recent physiological data on theta cycle skipping. The rebound spiking interacts with subthreshold oscillatory input to stellate cells or interneurons regulated by medial septal input and defined relative to the spatial location coded by neurons. The timing of rebound determines whether the network maintains the activity for the same location or shifts to phases of activity representing a different location. Simulations show that spatial firing patterns similar to grid cells can be generated with a range of different resonance frequencies, indicating how grid cells could be generated with low frequencies present in bats and in mice with knockout of the HCN1 subunit of the h current.


INTRODUCTION
An understanding of cortical function requires that network dynamics be linked to intrinsic cellular properties. The model presented here links cellular properties of resonance in neurons in the medial entorhinal cortex (MEC) to the properties of grid cells. Grid cells exhibit spiking activity as a rat forages in a regular array of locations that fall on the vertices of tightly packed equilateral triangles (Fyhn et al., 2004;Hafting et al., 2005;Moser and Moser, 2008). The model presented here links the difference in size and spacing of grid cell firing fields along the dorsal to ventral axis of MEC (Hafting et al., 2005;Stensola et al., 2012) to the intrinsic resonance properties of stellate cells in layer II of MEC along this axis Giocomo and Hasselmo, 2009;Boehlen et al., 2010;Heys et al., 2010;Canto and Witter, 2012;Pastoll et al., 2012a;Shay et al., 2012). These resonance properties depend upon the hyperpolarization activated cation current (h current) (Dickson et al., 2000;Haas and White, 2002;Erchova et al., 2004;Fransén et al., 2004;Rotstein, 2014) which also causes rebound spikes after hyperpolarization (Alonso and Llinas, 1989;Alonso and Klink, 1993). The model also addresses differences in resonance frequencies in layer II neurons from slices of MEC in the bat and the rat (Heys et al., 2013) and from mice with knockout of the HCN1 subunit of the h current ).
In one recent attractor model (Navratilova et al., 2012), intrinsic rebound spiking allowed reactivation of the attractor state on each theta cycle to generate spatial periodicity, and replicate the phase precession of grid cell spiking relative to theta rhythm oscillation. Earlier models used this mechanism to simulate shortterm memory (Lisman and Idiart, 1995) and theta phase precession (Jensen and Lisman, 1996). These models resemble the model presented here, but they used a depolarizing rebound based on spikes in an integrate-and-fire model (Navratilova et al., 2012). This previous dependence on spiking prevented explicit simulation of subthreshold properties caused by h current such as resonance and rebound from inhibition studied here, and the membrane potential oscillations seen in intracellular recordings of grid cells (Domnisoru et al., 2013;Schmidt-Hieber and Hausser, 2013). That previous paper also did not use an interaction with phases of network oscillations.
The model presented here builds on the oscillatory interference models of grid cells developed by Neil Burgess Burgess, 2008;Bush and Burgess, 2014) and variants of this model (Blair et al., 2007(Blair et al., , 2008Hasselmo, 2008;Zilli and Hasselmo, 2010;Welday et al., 2011). Oscillatory interference models (Burgess, 2008) are effective in simulating the theta phase precession of grid cells shown in one dimensional tracks (Hafting et al., 2008;Eggink et al., 2014) and two-dimensional environments (Burgess, 2008;Climer et al., 2013), building on the models of theta phase precession in hippocampal place cells (O'Keefe and Recce, 1993;O'Keefe and Burgess, 2005). In theta phase precession, grid cells start out firing at late phases of theta cycles when a rat enters the firing field of the cell and gradually shift to earlier phases of firing as the rat exits the firing field.
Oscillatory interference models could address the correlation between the size and spacing of grid cell firing fields at different dorsal to ventral anatomical positions (Hafting et al., 2005;Sargolini et al., 2006;Stensola et al., 2012) with the intrinsic spiking frequency measured with autocorrelograms of extracellular spiking in awake, behaving rats (Jeewajee et al., 2008;Stensola et al., 2012) and the intrinsic resonance frequencies measured with intracellular recording of membrane potential in slices Hasselmo et al., 2007;Giocomo and Hasselmo, 2008a). Changes in intrinsic properties due to knockout of the HCN1 subunit of the h current channel alters the spacing and size of grid cell firing fields . However, previous models do not directly address how intrinsic resonance can influence grid cell firing. Early oscillatory interference models suggested that membrane potential oscillations could interact within single neurons Hasselmo et al., 2007), but later analysis showed that subthreshold oscillations of different phase cannot coexist within a single neuron (Remme et al., 2009(Remme et al., , 2010 and do not change linearly with depolarization (Yoshida et al., 2011). Subsequent models showed oscillatory interference dependent on network interactions (Blair et al., 2008;Burgess, 2008;Hasselmo, 2008;Zilli and Hasselmo, 2010;Welday et al., 2011;Bush and Burgess, 2014). A related model used waves spreading across the cortical surface for coding location by firing phase and resulted in decoding activity with hexagonal interference patterns (Nadasdy, 2009(Nadasdy, , 2010. However, these models have not explicitly simulated how intrinsic rebound spiking could interact with network mechanisms to influence the size and spacing of grid cell firing fields. The model presented here focuses on the role of intrinsic resonance and rebound in causing phase shifts in spiking relative to network oscillations. This work extends beyond a previous model in which rebound spiking was used to simulate firing fields and theta phase precession along a one-dimensional trajectory (Hasselmo, 2013). This paper compares rebound to excitatory and inhibitory synaptic input, extends simulations to two dimensional trajectories, and addresses interactions of traveling waves with multiple heading angles.

LINK TO RESONANT PROPERTIES
A major focus was the simulation of resonance properties of single medial entorhinal neurons in response to current injection sweeping through increasing frequencies (Haas and White, 2002;Erchova et al., 2004;Fransén et al., 2004;Giocomo et al., 2007;Canto and Witter, 2012;Shay et al., 2012), as shown in Figure 1A. The current injection sweeping through frequencies was generated by the MATLAB chirp function and used to measure the impedance amplitude profile (ZAP). The resonance properties arise from intrinsic currents including the hyperpolarization activated cation current (h current) as described in layer II stellate cells of the MEC (Dickson et al., 2000). Previous models addressed these resonance properties (Dickson et al., 2000;Erchova et al., 2004;Fransén et al., 2004;Schreiber et al., 2004;Izhikevich, 2007;Engel et al., 2008;Rotstein, 2014;Rotstein and Nadim, 2014). The h current has a slow time course such that it responds to the slow hyperpolarizing phase during lower frequencies at the start of the chirp function with a more rapid depolarization that occurs at earlier phases relative to input and decreases membrane potential changes. Later when the frequencies in the chirp function match the time constant of the h current, the h current adds to the oscillation amplitude, causing a peak in the amplitude of membrane potential response. At higher frequencies of the chirp function, the h current response to hyperpolarization lags behind the input causing a decreased change in membrane potential.
Here the resonance properties of entorhinal neurons are modeled with linear coupled differential equations with oscillatory dynamics (Hasselmo, 2013). This differs from many previous oscillatory interference models that used sinusoids to represent oscillations Blair et al., 2007Blair et al., , 2008Hasselmo et al., 2007;Burgess, 2008;Hasselmo, 2008;Hasselmo and Brandon, 2008). The sinusoids in those models A C B FIGURE 1 | (A) Experimental data shows the response of a layer II stellate cell in MEC neuron to a chirp function (Shay et al., 2012). In response to oscillating current injection that increases in frequency, the neuron shows a gradual increase in amplitude of membrane potential oscillation that reaches a peak at resonance frequency and then decreases. (B) The chirp function input sweeps from 0 to 20 Hz over a period of 20 s. (C) The resonance frequency of the modeled stellate cell Equation (1) can be seen in the response to the chirp function. In this simulation, g p = 0.75, g h = 0.15, m = 1, and o h = 0.35, give resonance frequency of 10.2 Hz. could represent frequency and phase of oscillations but kept amplitude constant. Coupled differential equations allow the simulation of resonance frequency and strength in single neurons, as well as the change in response amplitude with circuit interactions.

RESONANT NEURONS
The equations of this simple model of resonance represent the change in membrane potential of an individual neuron v relative to resting potential (zero in these equations), and the change in activation h of the hyperpolarization activated cation current as follows: The membrane potential v has passive decay modeled by the parameter g p , for decay back toward resting potential. The h current is a cation current that depolarizes the cell, so when h current is positive, it depolarizes the cell proportional to the parameter g h . The h current h is turned off by depolarization, so when v goes to positive values, it decreases the magnitude of h in proportion to mv. The h current is activated by hyperpolarization, so when v goes to negative values, it increases the magnitude of h in proportion to mv. The h current decays according to parameter o h which was set to either 0.35 or 0.1. The mathematical properties of these equations are well described (pp. 89-97 of Hirsch and Smale, 1974;Rotstein, 2014;Rotstein andNadim, 2014, pp. 101-106 of Izhikevich, 2007). Here, parameters were chosen to give properties of resonance frequency that resemble the experimental data using the ZAP protocol. The dynamics of the network described below depend upon the resonance frequency of simulated stellate cells relative to the frequency of medial septal input f described below. The equations above can be algebraically reduced to the characteristic equation for a damped oscillator with forcing current: (3) This second order version of the equation could also be used to create a wave equation for transverse waves defined relative to the spatial location coded in the environment. The undamped resonance frequency of this equation is f = (g h m + g p o h )/2π . For example, with a time step of 0.01 s, the parameters g p = −0.49, g h = 0.24, m = −1, o h = −0.35 give f = 10.2 Hz. These parameters work well in Figures 3-7. However, the network dynamics also depend upon the strength of synaptic interactions, so the quantitative network dynamics cannot be determined only by Equations (1) and (2). Equations were solved in MATLAB using simple forward Euler methods, and qualitatively similar results were obtained using the ode45 solver (Runge-Kutta) in MATLAB.
The parameters were chosen to replicate resonance properties of stellate cells in layer II of MEC as shown in Figure 1A (Shay et al., 2012) in response to current injection consisting of the chirp function in Figure 1B, in which the frequency of the input current changes linearly from zero Hertz to 20 Hertz over 20 s. These functions are sometimes referred to as ZAP currents, where ZAP refers to the impedance amplitude profile computed in response to the chirp. In Figure 1C, a simulated neuron using the above equations shows a gradual increase in amplitude of oscillatory response to current injection until it reaches a peak response at the resonant frequency after which the amplitude of the oscillatory response decreases. This resembles the resonance response in the recording from a layer II stellate cell. The plot shown in Figure 1C used g p = −0.75, g h = 0.15, m = −1, o h = −0.35 give f = 10.2 Hz. However, it was more difficult to balance the network dynamics with g p = −0.75, so some network simulations used a lower value of g p . Figure 2 Column 1 shows how the model can replicate different resonance properties in experimental data with different values of parameters g h generating high (Figures 2A-C) and low resonance frequencies (Figures 2D,E). Examples 2C and 2E have the lowest resonance strength. The network model with excitatory connections below works best with the parameters shown in Figure 2A, but still works effectively with parameters shown in Figures 2B-D. The model with inhibitory connections works better across the full range of parameters. Figure 2 Column 2 shows that with depolarized initial conditions, neurons show rebound depolarization that could activate the next oscillation cycle of an attractor. These traces in column 2 resemble the sequence of afterhyperpolarization and afterdepolarization that occurs after spikes in the experimental data from intracellular recording from stellate cells Giocomo and Hasselmo, 2008b;Navratilova et al., 2012 (Navratilova et al., 2012). This rebound depolarization is weak in (C2,E2). (F) Input current for (A1-E1) consists of chirp stimulus sweeping through different frequencies.
resulting in faster afterdepolarization in stellate cells from more dorsal slices (Navratilova et al., 2012). Network attractor dynamics work well when rebound depolarization is strong, as shown in Figure 2A2 and is still effective with one cycle of rebound as shown in Figures 2B2,D2. Weaker resonance strength corresponds to weaker rebound depolarizations in these neurons. Figures 8-13 later show that network dynamics are easier to maintain with rebound due to hyperpolarizing inhibitory input rather than the depolarizing input shown here.
Another important property associated with resonance is the phase shift of the membrane potential in response to a sinusoidal current injection. When resonance frequency is higher than the input frequency, the simulated stellate cells show a phase difference in which the oscillations lead each cycle of the sinusoidal input. This difference in phase of the rebound depolarization can shift the phase of stellate cells relative to input from the medial septum to the network as described below, allowing a progressive shift in oscillatory activity within the network.

Network interactions of stellate cells
The network simulations presented here generate patterns of activity based on network interactions summarized in Figure 3. The stellate cells in the model are proposed to be grid cells, so they are labeled with G/S in the figures and their membrane potential is designated by v g . The network dynamics involve an interaction of the resonant properties of stellate cells with the activity of other neurons receiving oscillatory input at different phases that are directly or indirectly regulated by rhythmic input from the medial septum. The medial septum has been shown to regulate theta rhythm oscillations in the entorhinal cortex (Mitchell et al., 1982;Jeffery et al., 1995;Brandon et al., 2011) and contains neurons that spike coherently with theta rhythm oscillations at a wide range of phases (Bland, 1986;King et al., 1998). As noted above, the network interactions of stellate/grid cells are analyzed in two different versions of the network model: (1) The excitatory model, and (2) The inhibitory model.

Excitatory model
The excitatory model focuses on excitatory connections of the stellate cells with excitatory pyramidal cells in layer II or in deeper layers of entorhinal cortex such as layer III and layer V, as shown in Figure 3A. In this model, the pyramidal cells are proposed to be conjunctive cells, so they are labeled with C/P in Figures 3-7, and their membrane potential is designated by the v c . The dynamics involved subsets of grid/stellate (G/S) cells interacting with subsets of conjunctive (C/P) cells via excitatory connections between these two classes of neurons. The activity was maintained in one location when the timing of rebound activity of grid/stellate cells matched the peak phase of activity of conjunctive/pyramidal cells driven by rhythmic input from the medial septum. This can be considered a type of attractor, in which a wave of excitatory activity caused by recurrent excitation alternates with a period of hyperpolarization that activates the h current in stellate cells and causes rebound depolarization that reactivates the same pattern of stellate activity from the previous cycle. The activity can also shift when the rebound activity in stellate cells is shifted in phase relative to the currently active pyramidal cells.

Inhibitory model
The inhibitory model focuses on inhibitory input from interneurons with local focused connectivity to the stellate cells, as shown in Figure 3B. The simulations shown in Figures 8-13 focus on this inhibitory input to stellate cells. The network simulations use a single population of inhibitory cells that interact with the stellate cells, but these inhibitory neurons could also receive head direction input from head direction cells (Taube et al., 1990), or could be activated by conjunctive/pyramidal cells in deep layers that respond to head direction input. Therefore, the inhibitory interneurons could also be considered to have conjunctive cell properties consistent with experimental data. The inhibitory model in Figures 8, 9 was simulated with continuous firing rate equations whereas the inhibitory model in Figures 10-13 used a spiking model of layer II stellate cells in MEC developed by Izhikevich (2007) and described below.

Equations of the model
Simulations of the excitatory model in Figures 3A, 4-7 and the initial inhibitory model in Figures 3B, 8, 9 use a continuous firing rate representation of individual neurons based on the equations of resonance above. The equations representing individual neurons were replicated and combined into local circuit interactions (see Figure 3) as follows: Where v s is a vector representing the different membrane potentials of the individual grid cells/stellate cells (with index s) within a population of grid/stellate (G/S) cells with resonance properties due to their intrinsic variables h s [which represents the h current as described in Equations (1) and (2)]. v i,p is a vector representing the different membrane potentials of a different population of cells that in the excitatory model are referred to as conjunctive/pyramidal (C/P) cells with each cell indicated by index i, and different sets of C/P cells with different preferred head directions θ p designated by p. In the inhibitory model the activity of the conjunctive cells provide output not to the stellate cells but to interneurons with index i that provide one to one connectivity back to the stellate cells. The function H [.] is a Heaviside function, that is zero below threshold η and one for values above η. The matrix W si is an identity matrix, so that each conjunctive cell or interneuron i is connected back to a single grid cell at the same position with the same index s (but the input is summed over the active conjunctive cells using the index p). The matrix W is provides connectivity from each grid cell to all the cells within a single "line" of the conjunctive cells (that send output to interneurons) with an index of the spatial location coded by the neurons that overlaps with the position of a given grid cell. Figures 12, 13 showing the two-dimensional simulation in the Results section illustrate the functional effect of this connectivity pattern in which each single stellate cell contacts a "line" of conjunctive cells.

Medial septal input
As can be seen in Equation (6), the membrane potential of the second population v i lacks effects of the h current and therefore lacks resonance properties. However, this population receives a driving rhythmic (sinusoidal) input from interneurons that represents neurons regulated by oscillatory input from the medial septum (MS). Experimental data suggests that the medial septum primarily contacts inhibitory interneurons in the hippocampus (Freund and Antal, 1988) and in the entorhinal cortex (Witter, personal communication). Therefore, in the excitatory model, the MS input represents input to C/P cells from local interneurons that receive oscillatory modulation from the medial septum. In the inhibitory model, the MS input in the equation represents how local interneurons gate the activity of conjunctive cells that drive the inhibitory interneurons in the model (see below). The sinusoidal input to v i regulated by medial septum in both models is represented by I ms with magnitude μ = 2.8. The frequency and phase of the sinusoidal input is determined by the following equation with x evenly distributed between 0 and 2π so each value x i = 2π i/n, where i is the integer index of the neuron in that line and n is the total number of neurons with different evenly spaced phases along the dimension of the heading direction selectivity of the conjunctive cells.
This corresponds to a traveling wave where f is the medial septal frequency and k is the wave number (number of spatial cycles across the simulated population) . The frequency f of medial septum was set to 8.1 Hz in Figure 4 and wave number was set to 1. This equation could be seen as a solution (Elmore and Heald, 1969) to a medial septal wave equation: δt t These wave functions could occur for a wide range of heading angles as shown in equation 10 below.
The network dynamics depended on the resonant frequency determined by the parameters of the above equations relative to the frequency f of the medial septal input. When the resonant frequency was significantly different from the medial septal input, this caused a shift in activity between different neurons in the population as shown in the figures.
The speed of shift in activity could be regulated by altering the frequency of medial septal input (as one possible mechanism). As shown in the Results section, this provides a smooth influence of movement speed on the shift of activity in the network. Experimental data has shown that running speed is associated with changes in both frequency of network theta rhythm (Maurer et al., 2005) and firing rate of entorhinal neurons Wills et al., 2012). The even distribution of phases of medial septal input is motivated by recording of unit activity in the medial septum showing a broad distribution of different phases of firing relative to hippocampal theta rhythm (Bland, 1986;King et al., 1998).

Head direction input
The conjunctive cells receive depolarizing input I hd based on the current head direction θt of the animal relative to the preferred head direction θ p of each population of C/P cell. This simulates the fact that conjunctive grid-by-head-direction cells respond to both the current head direction of the animal and the position of the animal relative to grid cell firing fields (Sargolini et al., 2006;Boccara et al., 2010;Brandon et al., 2011;Stensola et al., 2012). This means that changes in the relative frequency of resonance and MS input will selectively influence the relative firing of conjunctive cells for the current head direction, providing a direction component to interact with the speed component.

Pattern of connectivity
Most of the simulations presented here utilize arrays of neurons representing interactions along a one-dimensional trajectory through space with only two directions of movement along this trajectory, as in other recent models (Navratilova et al., 2012). However, simulations in Figures 12, 13 also address firing in twodimensional environments with a full range of head direction angles as in other models (Burgess, 2008;Hasselmo and Brandon, 2012), so the description here will describe two dimensional arrays.
In simulations of the excitatory model there was a single population of G/S cells. In simulations of the inhibitory model, there were two populations of G/S cells (labeled with subscript indices "T1" and "T2"). In contrast to the G/S cells, the conjunctive cell populations are organized into head direction selective planes that consist of two-dimensional arrays of cells with dimensions x and y corresponding to position within the array (these positions differ by head direction angle relative to other arrays). In the simulations shown in Figures 12, 13, there were four populations of conjunctive cells for each of the populations T1 and T2 of G/S cells. Each of these 4 populations of conjunctive cells drove a single population of inhibitory interneurons providing input back to the stellate cell populations. The conjunctive cells received phasic input from medial septum assigned relative to their x and y position and their head direction selectivity θ p (this MS input was treated as separate populations, but would most likely contain overlap for neurons of the same phase). The total size of dimensions x and y were the same for the single population of G/S cells, the multiple populations of C/P cells and the input from medial septum (MS).
Within each plane, if the current head direction θ is close to the preferred head direction θ p of that plane, then all the conjunctive cells receive the same input for preferred head direction I hd = 1.0. The cells can be considered to be organized into lines of cells along the dimension of head direction preference, consisting of x different lines in each plane and y cells in each line. Within each line, each of the y C/P cells receive input of strength W is = 0.5 to 0.53 from each of the full array of G/S cells in that line. This input is generated whenever the cell v s is over threshold θ s = 0.2, where the brackets [] represent the Heaviside step function which has value 0 for v s < θ g and value 1 for v s > θ s . Each G/S cell v s receives excitatory feedback input from the single corresponding C/P cell in population v i with the same Heaviside step function computed for values above θ i = 4. The feedback input to grid cells from conjunctive cells arrives via an identity matrix (one to one connectivity) with uniform connection strength W si that varied dependent upon the resonant strength. When resonance strength was high the connections were set at 0.85 (e.g., Figure 5), whereas lower resonance strength required stronger connections up to 2.8.

Global inhibition
The network also included inhibition with global connectivity that helped to increase stability in the simulations. This global inhibition used a variable z i that has a slow buildup similar to the time course of GABA B receptor activation of potassium conductances. Each inhibitory unit was modeled as receiving input from the full population of G/S cells. Each G/S cell received uniform input of strength H from the inhibitory unit z within that head direction selectivity, so that this corresponded to uniform inhibition for a given line of cells with a slow time course.

Spiking neuron simulations
The simulations shown in Figures 10-13 used the model of simplified spiking neurons developed by Izhikevich (2007). In these models, excitatory entorhinal stellate cells were simulated using the parameters (Izhikevich, 2007, p. 315) that were tuned to replicate experimental data on stellate cells from John White's lab (Burton et al., 2008). These neurons were combined with more abstract representations of conjunctive cells and inhibitory interneurons. The neurons used the following equations (Izhikevich, 2007, p. 273): For stellate cells, the parameters were: Consistent with experimental data (Couey et al., 2013), there was no direct synaptic coupling between stellate cells. Based on data on theta cycle skipping (Jeffery et al., 1995;Deshmukh et al., 2010;Brandon et al., 2013), the stellate cells in the simulations in Figures 10-13 were split into two populations S T1 and S T2 that fired on alternating cycles of the theta rhythm. Synaptic coupling between the stellate cells and conjunctive cells used step functions that were active for one time step after a spike, in order to preserve the detailed timing of spikes relative to the phase of medial septal traveling wave input to conjunctive cells. The medial septal input was structured as square wave inputs that did not overlap in the depolarized phase, so only a limited number of conjunctive cells at a given time or phase were sufficiently depolarized to be able to respond to the synaptic input from the spiking stellate cells. This was done separately for two groups of four different populations of conjunctive cells representing different directions (e.g., C T1,W represents a conjunctive population associated with stellate cell population T1 with a west-moving traveling wave). These groups of four directions of conjunctive cells drove single interneuron populations that provided inhibitory feedback to the stellate cell populations. The spiking model had the interesting characteristic that running speed could be represented by the strength of feedback inhibition, because the speed of rebound spiking depended upon the magnitude of hyperpolarizing pulses. This is discussed in more detail in the Results section.

MODEL WITH WIDE RANGE OF WAVE DIRECTIONS
Understanding the role of resonance in entorhinal cortex neurons requires understanding the response of entorhinal neurons in terms of oscillations in time. Understanding how resonance could influence the representation of different spatial phases by grid cells requires the use of oscillations of different phases coding different spatial locations. The need for both temporal phase and spatial phase indicates that the network might be modeled in terms of waves of activity. Note that the spatial dimension of these waves is not the anatomical position of cells within the cortex, as in many previous models (Ermentrout and McLeod, 1993;Coombes, 2005;Meijer and Coombes, 2013), but in terms of the location within the environment coded by individual cells. The simulations in Figures 14, 15 focused on showing how the feedback interaction with populations of traveling waves with many different directions could result in a dominant influence of waves with a limited range of directions.
As an initial example of how interacting waves could generate the grid cell firing pattern within a population, Figures 14,  15 show a model in which waves with multiple different directions of propagation have feedback interaction with a population of entorhinal cells to generate a grid cell firing pattern within the population and within individual cells. The wave in this model corresponds to the models by Blair and Zhang showing that input from ring attractors representing oscillatory input from theta cells in the medial septum or thalamus can generate grid cell firing patterns (Blair et al., 2008;Welday et al., 2011). In the Blair papers, the oscillating rings are referred to as ring attractors and use one dimensional rings of neurons. In the equations below, the ring attractors are described as generating two-dimensional transverse traveling waves that propagate across the plane of grid cells representing an environment. The representation of animal velocity in these models depends upon shifts in the pattern of interference between oscillations, so the model falls in the category of oscillatory interference models developed by Burgess et al. (2007), Burgess (2008).
Modelers that use attractor dynamics with circularly symmetric synaptic feedback to generate grid cells have criticized oscillatory interference models for assuming that the velocity controlled oscillators must be constrained to have preferred heading angles at intervals of 60 • . However, the simulations presented in this section show that an oscillatory interference model can generate hexagonal grid cell patterns with velocity controlled oscillators driven by inputs with a wide range of arbitrary input angles. This was shown using a variant of the oscillatory interference model presented in an earlier publication in which multiple directions of two-dimensional transverse traveling waves interact to generate grid cells (Hasselmo and Brandon, 2012).
Several equations were used to describe the model in that previous paper (Hasselmo and Brandon, 2012), but the function of that previously published model can be greatly clarified here in a different mathematical description that allows summary with a single equation: Where g xy (t) represents the activity at time t of an array of grid cells that code a set of locations x,y in the two-dimensional spatial environment being explored by the rat. In Equation (10), the grid cell activity is updated by input from a large array of transverse traveling waves (referred to as heading angle planes in the previous paper) that provided input to the two-dimensional array of grid cells. As in the previous model, each location in Equation (10) is represented by a full set of traveling waves coding a range of heading directions φ (each with a different index i) in the rat environment. The angle φ of propagation of each traveling wave is represented by the direction of a unit vector κ, related to the direction angle φ and wave number k as follows κ This is a standard representation for the direction of two-dimensional transverse, traveling waves (Elmore and Heald, 1969). Thus, each transverse, traveling wave propagates across the plane of spatial locations in a different direction φ. The traveling waves have temporal frequency ω. The equation shows how the input to each grid cell g xy (t) is influenced by the sum of traveling waves arising from the full array of locations x y with the full set of angles φ at each location. The spatial phase of each of these waves is computed by the dot product of the vector κ(φ i ) with the vector computed from the difference between the full array of location vectors r xy = [x y] and the vector representing the central location of that traveling wave r x y = [x y ]. This results in the spatial phases κ (φ i ) ( r xy − r x y ).
Each of the arrays of traveling waves with spatial phases corresponding to the location x y with the set of angles φ had a magnitude in proportion to the level of previous grid cell activity g x y (t − 1) coding each of the locations x y at time t-1. The grid cell activity was put through a threshold linear function with threshold . Equation (10) shows that in the model in the previous paper (Hasselmo and Brandon, 2012), the set of traveling waves corresponding to each location x y are first summed over all of the angles φ i , giving concentric circles around the location x y . Subsequently, the grid cell activity g xy (t) at the current time step t was determined by summing over all of the spatial phases corresponding to summation over all locations x y (plus the current activity scaled by τ = 0.3).
In simulations of that previous model, a set of twodimensional traveling waves was generated with a range of 24 different heading angles evenly distributed over 360 • and a wave number of 4 across the full environment being coded by a rat (Hasselmo and Brandon, 2012). The full set of 24 traveling waves was generated for each coded location, and summation generated a circularly symmetric standing wave pattern. The magnitude of the standing wave was then regulated by the activity of each individual grid cell on the previous time step. Thus, the previous model performed summation over spatial locations subsequent to summation over heading angles. The interaction of the grid cells and the traveling waves resulted in population activity settling into a pattern of grid cell activity as shown in multiple figures in that paper (Hasselmo and Brandon, 2012). One important problem with this model is that the circularly symmetric standing waves generated by summation across angles would have different spatial locations of peak firing on different half cycles. This was avoided by the use of two populations with π phase difference, and the use of only half the cycle for each. The mechanism for this constraint was not simulated but could be achieved by cross-inhibition between the two populations.
The computation of this model can be made more efficient by changing the order of summation to speed the computation process and to describe the input as a smaller array of traveling waves. The new mathematical description of the previous model allows simplification by changing the order of the summation and by adding a different mechanism for movement. This allows the previous paper to be extended by showing simulations for much larger range of interacting angles and full two-dimensional trajectories that were not previously tested.
Thus, the summation can be performed as follows: In Equation (11), for the traveling waves coding a particular heading direction, the waves are summed over all the phases corresponding to the grid cells coding all locations x y . In addition, the computation can be further speeded if instead of explicitly computing and summing the full set of traveling waves, the circular mean phase and mean amplitude are computed across the full array of waves representing the angle φ i at every location x y . This different order of summation first generates the amplitude and phase of a single traveling wave for each heading direction angle φ i , and then sums these waves across the full set of direction angles. The waves generated before summation across angles can represent the input from a single cyclical traveling wave for each angle, which could be generated by a mechanism such as the previously proposed ring attractors (Welday et al., 2011), instead of requiring a full array of ring attractors for coding of each location. As noted above, after computing the phase of each traveling wave across the population, the grid cell activity g(t) at the current time step is determined by the summation over all of the angles φ i of the different traveling wave input. This allows the input from transverse traveling waves with a large range of preferred directions (heading angles) to generate an interference pattern in the grid cell population that settles into a hexagonal grid cell firing pattern, and ensures that all the waves have the correct relative phase when location is updated, allowing a shift in the driving force of velocity when the animal turns. Note that this model still requires the use of two sets of waves offset by phase π to avoid the shift in coded location. A gap of π /8 between the waves was used in simulations shown here, τ was set to 0.5, and the normalization factor was squared to reduce overall amplitude. Note that the direct feedback excitation used here and in Hasselmo and Brandon (2012) differs from the modification of synaptic input used in models of the self-organization of directional input (Mhatre et al., 2012) or place cell input (Kropff and Treves, 2008;Si et al., 2012).
Equation (11) includes an additional element H(φ i ) x(t) that represents the effect of the current velocity x(t) of the animal on the phase of individual waves. The effect of velocity on each wave is scaled by one row of the matrix H(φ i ) which computes the cosine and sine for each heading angle φ i with H(φ i ) = [cos φ i sin φ i ]. Because the network updates the phase of all waves, the location update also works with narrower functions more similar to the tuning of real head direction cells, such as rectified cubesH(φ i ) = [cos 3 + φ i sin 3 + φ i ]. Note that equations 10 and 11 are using cosine functions for waves rather than simulating them in a neural circuit with the rebound properties of neurons, but these traveling waves can be considered a potential solution to wave equations for this network.

RESULTS
The system summarized by Equations (4)-(8) provides dynamics that allow progressive shifts of activity in grid/stellate (G/S) cells and conjunctive/pyramidal (C/P) cells. Figure 4A shows the conditions without feedback between the C/P and G/S cells. The input from the medial septum (MS) shown in Figure 4A drives oscillations that remain subthreshold in a single C/P neuron v i (Figure 4A row MS). Starting from a depolarized potential, the resonance properties of the G/S cell causes a single rebound depolarization that brings it above threshold. In the absence of excitatory connections, the C/P cells do not come over threshold. Figure 4B shows the effect of bidirectional excitatory connections (arrows in circuit diagram). The G/S cell starts with depolarized initial conditions, but the addition of excitatory connections allows this activity to drive the C/P cell over threshold and the neurons drive each other. In contrast to fixed-point attractors, the activity does not persist in a static manner. The activity decreases because the MS input to the C/P cell decreases in amplitude and because the depolarization shuts off the h current in the G/S cell, allowing the G/S cell to drop below threshold and subsequently later show a rebound depolarization that brings it above threshold again. The rebound depolarization over threshold causes further output that drives the C/P cells over their threshold (gray dashed line). The C/P cell excitatory feedback drives the G/S cell to higher activity. As the MS input to the C/P cell decreases and the h current shuts off, both the C/P cell and the G/S cell fall below threshold again, but the G/S cell shows rebound depolarization on the next cycle, driving the C/P cell over threshold and pushing the G/S cell to the same depolarized level. In this manner, the two cells continue to excite each other and maintain a repeating oscillatory interaction throughout this example. This does not require continuous activity as in traditional fixed-point attractors. Instead, oscillatory activity is maintained with low firing rates across theta cycles by the use of intrinsic rebound depolarization. This resembles the rebound spiking used in Navratilova et al. (2012), but uses subthreshold dynamics that allow simulation of subthreshold resonance properties.

MOVEMENT OF ACTIVITY BETWEEN MULTIPLE STELLATE-PYRAMIDAL PAIRS
A shift in the oscillatory activity can occur when multiple C/P neurons with different phases of input from medial septum interact with multiple G/S neurons, as shown in Figures 4C, 5A. This allows a shift in the phase of maximal activity that is proportional to the difference in resonance frequency in the G/S cells relative to the MS input frequency to the C/P cells. In this example, the rebound depolarization occurs early relative to the MS input oscillations to C/P. This corresponds to the phase shift caused by the currents underlying resonance. Figure 4C shows how this difference between resonance frequency and MS input causes a shift in phase for a single G/S neuron and relative to a single C/P neuron when both are part of a network shown on the left. The black lines in the circuit diagram show the neurons illustrated in this example, and the gray lines in the circuit diagram show other neurons present in the simulation whose traces are not shown. In this network, the C/P neurons each receive the same frequency of MS input, but with a different phase. The G/S neurons each receive input from their corresponding single C/P neuron (black upward arrow), but they send output to the full line of C/P neurons at all phases (downward line that forks into many branches).
The G/S neuron in Figure 4C starts out above threshold, and interacts with the C/P neuron shown here, but on each cycle its rebound depolarization occurs a bit earlier, causing it to send stronger output to another C/P cell that receives MS input with an earlier phase. Because the C/P cell with earlier phase of MS input is closer to threshold, this causes it to generate greater output, which does not go to the G/S cell shown here, but goes to the G/S cell associated with the earlier phase C/P cell. Thus, the activity transitions to a C/P cell with an earlier phase of MS input, and the corresponding G/S cell. As the first G/S cell shown here goes out of phase with the MS input to its corresponding C/P cell, the amplitude of oscillation in both of these neurons decreases over the first 5-6 cycles.
The activity transitions to different pairs with different phases until it cycles through the full network and starts to activate the C/P and G/C cells illustrated here to start an increase in the envelope of oscillations (that appears as a beat pattern). This beat pattern looks like an oscillatory interference model, but the model presented here uses intrinsic oscillations (modeled by differential equations) rather than fixed sinusoidal functions for the interaction. In these differential equations, the beat frequency depends upon the difference in resonance frequency of the G/S cell and the frequency of MS input to the C/P cell. This links the properties of oscillatory interference models directly to the resonance properties of layer II stellate cells in entorhinal cortex.
The shift in phase of the G/S cells relative to the MS input resembles theta phase precession of MEC grid cells as shown in the oscillatory interference model (Burgess, 2008) and in experimental data from the Moser laboratory on a linear track (Hafting et al., 2008). Theta phase precession also appears in twodimensional data based on the distance since a rat enters a grid cell firing field . In the figure, precession also occurs in the C/P cell consistent with precession appearing in some conjunctive cells in deeper layers  including layer V pyramidal cells (Misuzeki et al., 2009). However, many layer III conjunctive cells do not show precession (Hafting et al., 2008). Figure 5 shows the transition of activity in a full population of six G/S cells and six C/P cells receiving input from six MS cells with different phases. A difference in resonance frequency relative to the MS input frequency is present in all the G/S cells in Figure 5B. The parameters are g p = −0.49, g h = 0.24, m = −1, and o h = −0.35, resulting in resonance frequency f = 10.2 Hz relative to MS input frequency of 8.1 Hz. The full array of G/S cells all have the same resonance frequency. Initial conditions start activity in cells G/S1 to G/S3. These cells show rebound depolarization that is early relative to their coupled C/P cells, so that their output activates the C/P cell with an earlier phase (C/P4) and this activates the associated G/S cell (G/S4). As this process continues, it progressively activates C/P and G/S pairs with earlier phase (e.g., C/P5 and G/S5), causing a shift in the amplitude of activity to neurons with earlier phase.  (6)]. The C/P cells receive oscillatory MS input with different phases [I ms in Equation (6)]. Each C/P cell sends output to only one G/S cell in the line, but receives input from all G/S cells in the line.

(B)
The resonance frequency of G/S cells (10.2 Hz) is higher than the driving frequency from MS cells (8.1 Hz), the rebound oscillations of G/S cells activate C/P cells receiving MS input of an earlier phase, causing a shift in activity from C/P1 and G/S1 to CP/2 and GS/2 and progressively to cells of earlier and earlier phase. When the activity reaches C/P6 and G/S6, the rebound activates the cells of next earliest phase, which are C/P1 and G/S1, so the cycle repeats. (C) When the resonance frequency of the G/S cells is reduced to a lower value (9.7 Hz) and closer to the frequency of the MS input (8.1 Hz), then the rebound oscillations match the phase of MS input in the same population and the activity remains stationary in the same set of cells 1-2 and 5-6.
When the activity reaches the earliest phase pair C/P6 with G/S6 at the top of each line, it starts to activate the first pair in the line which is the next earliest phase. In this manner, the cyclical nature of phases allows activity to be propagated repeatedly through the network, providing a rationale for the repeating firing fields of grid cells that requires a twisted torus connectivity in fixed point attractor models. Note that this process can work in the complete absence of inhibition, because the shift in phase of activity allows C/P neurons to drop below firing threshold, and terminates the end of the activity distribution across the population. However, addition of inhibition in the model makes the mechanism more robust across a wider range of parameters.
With lower resonance frequency that is closer to the MS input frequency (or with a different MS input frequency), the network can maintain a stable activity state in a single location, as shown in Figure 5C. The parameters are all the same except that parameter g h is set to 0.206 in Figure 5C so that resonance of G/S cells is 9.7 Hz. Note that even though there is still a difference in frequency, the network dynamics end up being stationary. In this example, cells 1, 2, 5, and 6 maintain activity for an extended time period (beyond the period shown), and the other pairs remain inactive. Note that C/P3 and C/P4 still get input from the active cells, but input shows destructive interference with MS input, and therefore G/S3 and G/S4 are primarily inactive (flat lines) except for occasional bumps of activity.
The shift in activity across the population of neurons shown in Figure 5B is considered to represent the movement of the animal in one direction. If the difference in frequency in Figure 5 represents the speed of the rat as it moves in a straight line in one direction, the activity would transition between different elements of this population of neurons to represent the shift in location. This resembles how oscillatory interference models generate bands of constructive interference due to interactions of velocity controlled oscillators Burgess, 2008). In this simulation, the relative phase inputs from medial septum provide a mechanism for competition between neurons.
The grid cells could interact with different conjunctive cells oriented in different directions, as in the model using conjunctive cells with six directions of tuning described by Burgess (2008).
As shown in Figure 2, the magnitude of resonance strength corresponds to the magnitude of rebound depolarization. Because of this, the magnitude of resonance strength corresponds to the ease of simulating oscillatory interference. However, as shown in Figure 6, the oscillatory interference beat patterns can be generated with a range of different parameters even corresponding to the magnitude of resonance strength that resembles the experimental data.
Note that the capacity for generating beat patterns does not depend upon the intrinsic resonance frequency of the neuron. The low resonance frequency shown in Figure 6D can generate a beat pattern similar to the beat pattern generated in Figures 6A,B, as long as the network frequency is reduced to interact with the resonance frequency. The spatial periodicity of the beat pattern does not depend upon the baseline resonance frequency or the frequency of input from the medial septum, but instead depends upon the magnitude of the difference between the stellate cell resonance frequency and the medial septum input frequency. Thus, the difference in resonance frequency found in bats and rats when recording from layer II neurons in MEC does not mean that resonance frequency does not contribute in some manner to the generation of grid cell spatial periodicity. However, the framework shown in Figure 6D would require lower frequency input from medial septum in the bat. The possible presence of such lower frequencies has not yet been tested in published results from bats.
The model can account for different spacing of grid cell firing fields dependent on the mapping of running speed to resonance frequency. In simulations, the change in resonance frequency can be achieved by altering the parameter g h in the first set of equations describing the continuous firing rate model. However, one difficulty is that the change in speed differs for the speed relative to different directions of movement. This is less of an issue if the conjunctive cells represent all directions of movement at high resolution, but it is a factor if there are a smaller number of directions represented. As an alternative, changes in speed can be coded by changes in the frequency of medial septal (MS) input, as shown in Figures 9, 11B, or possibly in the magnitude of the synaptic feedback inhibition within the circuit as shown in Figure 11A.
As shown in Figure 7A, if there is a large shift in resonance frequency with running speed, then the large difference between resonance frequency and MS input causes a rapid beat frequency. This would cause narrow spacing between smaller grid cell firing fields, consistent with the recording of higher resonance frequency found in dorsal MEC Giocomo and Hasselmo, 2008a,b;Boehlen et al., 2010;Pastoll et al., 2012a) where grid cell firing fields are smaller and closer together (Hafting et al., 2005;Sargolini et al., 2006). In contrast, if there is a small shift in resonance frequency with running speed, the smaller difference from MS input frequency causes a slower beat frequency. This could cause wider spacing between larger grid cell firing fields, consistent with the lower resonance frequencies found in ventral MEC Boehlen et al., 2010;Pastoll et al., 2012a) where grid cell firing fields are larger and spaced further apart (Hafting et al., 2005;Sargolini et al., 2006). The simulations in Figure 7 used global feedback inhibition to enhance stability. This involved a buildup of inhibition proportional to k = 0.05 for each time step that a stellate was above threshold η, with a decay of λ = 0.01. This influenced the network with inhibitory magnitude μ = −0.0928. Changing resonance frequency with running speed causes a shift in the beat frequency. However, shifts in resonance frequency with changes in speed do not give very stable transitions. For example, in simulations not shown, the running speed alters the parameter g h systematically, and MS input frequency is 8.1 Hz, using G/S activity with g p = −0.4585 or g p = −0.75 during a smooth increase in g h directly scaled to time (simulation not shown). The transitions are moderately consistent using g h = 0.2064 (f = 9.64 Hz) and going to g h = 0.2364 (10.03 Hz). However, the transitions become rather random starting with g h = 0.15 (10.22 Hz) and increasing to g h = 0.155 (10.28 Hz). Even with a smaller range of frequencies, the network breaks down with very small changes in the parameter g h . Thus, the speed modulation is unstable in the excitatory model with parameters that match the experimental data in Figure 1.
Because the beat patterns depend upon the difference in resonance frequency and MS input frequency, the network models in Figure 7 can generate a full range of beat frequencies that correspond to the spatial periodicity of a grid cell model. These one-dimensional spatial periodicities could potentially be combined together to generate a two-dimensional grid cell pattern in the same manner as described for sinusoid functions in Burgess (2008). As in the oscillatory interference model, the two-dimensional spatial periodicity in the model depends upon the scaling of resonance frequency to running velocity. Thus, this model could use elements of oscillatory dynamics to link the intrinsic property of resonance frequency in stellate cells to the features of grid cell responses, allowing a model to simulate the experimental data on larger spacing between firing fields in more ventral entorhinal cortex (Hafting et al., 2005) based on the lower resonance frequency of ventral stellate cells .

SIMULATION OF MOVEMENT WITH INHIBITORY FEEDBACK
The excitatory connectivity used in Figures 1-7 requires a biphasic rebound response of grid/stellate cells. The G/S cells are driven to strong firing by excitatory feedback from C/P cells. This shuts off the h current and causes the cells to hyperpolarize below threshold. This activates the h current, causing an afterdepolarization that brings the cells over threshold again. As shown in Figure 2, this rebound is more robust for a certain range of modeling parameters.
In contrast, inhibitory input can cause a monophasic rebound, as the cell can directly respond to inhibitory hyperpolarization by a rebound excitation. This could allow more effective rebound and network function with resonance parameters that do not work for excitatory feedback. Therefore, the model was tested with inhibitory feedback mediated by connections from inhibitory interneurons that connect selectively to the associated grid/stellate cells as shown in Figure 3. The connection from stellates to interneurons could occur via conjunctive/pyramidal cells to provide the head direction sensitivity. Figure 8 shows the function of a model using inhibitory feedback in a matched pair of cells as shown in Figure 3C. Because immediate feedback inhibition would counteract the spiking activity of a cell, the inhibitory feedback cannot be immediate, but instead uses the delayed rebound effect in a different cell pair to cause delayed feedback inhibition. As shown in Figure 8A, this mechanism requires coupling each G/S and interneuron pair (designated with subscript T1 in the figure, indicating one theta cycle) with a second G/S and interneuron pair (designated subscript T2 in the figure, indicating a second theta cycle). The pairs are defined by one to one coupling of the interneuron (e.g.,  Int T1 1) to the G/S cell in the same pair (e.g., G/S T1 1), whereas crossed connections link the G/S cell (e.g., G/S T1 1) to the opposite interneuron (e.g., Int T2 1). This circuitry results in a pattern of spiking on alternate theta cycles that has support from recent unit recording data in our laboratory (Brandon et al., 2013). This model also predicts that coupled pairs of neurons firing on opposite cycles should share the location of their firing fields. The dynamics of this circuit are shown in Figure 8B. In this figure, Grid/Stellate cell S T1 1 starts out hyperpolarized causing it to rebound above firing threshold (dotted line) in step a of Figure 8. The excitatory synaptic output from G/S cell S T1 1 possibly via a C/P cell depolarizes interneuron I T2 1 to bring it above threshold (step b). The output from inhibitory interneuron I T2 1 then causes hyperpolarization in the G/S cell S T2 1 (step c) causing a rebound depolarization (step d). G/S cell S T2 1 sends excitatory output (possibly via a C/P cell) to interneuron I T1 1. This brings interneuron I T1 1 over threshold (step e) causing hyperpolarization in G/S cell S T1 1 (step f). The rebound of G/S cell S T1 1 from hyperpolarization starts the entire cycle again, so that the oscillatory interaction is sustained within the circuit. This provides a form of maintenance of oscillatory activity that maintains activity using the resonance properties mediated by h current in stellate cells to provide the rebound from hyperpolarization. Note that both the G/S cells and the interneurons in the circuit show firing on opposite cycles of an underlying network rhythmicity, consistent with data on theta cycle skipping in the MEC (Jeffery et al., 1995;Deshmukh et al., 2010;Brandon et al., 2013) and input from cells in the medial septum that show theta cycle skipping (King et al., 1998;Varga et al., 2008). The maintenance of activity due to feedback inhibition resembles previous mathematical models of activity maintained by rebound spiking using different currents (Wang and Rinzel, 1993;McCarthy and Kopell, 2012).
These inhibitory dynamics provide the same capabilities for transition between attractor states as the excitatory connections in the earlier figures, as shown in Figure 9. Figure 9 shows the neural activity across a set of 6 different pairs of grid/stellate cells (both the S T1 1 and S T2 1 grid/stellate cells are shown) and across 6 different associated interneurons (again both the I T1 1 and I T2 1 interneurons). The shift in activity is caused by increasing the parameter g h in the model so that the resonance property has frequency 10.2 Hz and therefore the rebound from inhibition is faster relative to the frequency of MS input (9.75 Hz) to the pyramidal cells, causing a progressive shift in activity. Figure 9D shows that a progressive shift in the MS input frequency changes the speed of transition between different populations.

SIMULATION USING REBOUND SPIKING WITH INHIBITORY FEEDBACK
The rebound from inhibition that works in the abstract model also works in a network model using a more biologically realistic single neuron model of spiking and resonance properties in MEC developed by Izhikevich (2007). The single neuron model simulates the subthreshold resonance and oscillatory dynamics of medial entorhinal stellate cells, and shows resonance similar to the abstract model. Figure 10 shows that this model can function similar to the abstract model. Activity is initially induced by a hyperpolarizing input to stellate cell S T1 1 (label a) that causes a rebound spike. The stellate cell spike activates interneuron I T2 1 (label b) which causes an inhibitory hyperpolarizing potential in stellate cell S T2 1 (label c). This hyperpolarizing input causes a rebound spike in stellate cell S T2 1 after a delay. This rebound spike activates interneuron I T1 1, which sends inhibitory input to stellate cell S T1 1 (label d). This causes a rebound spike in the first stellate cell S T1 1 to restart the cycle. In this manner, A hyperpolarizing pulse (labeled with a) given only to G/S S T1 1 causes a rebound spike (solid black). This brings Interneuron I T2 1 above threshold (dashed black line) by label b. The interneuron sends inhibitory input to G/S S T2 1 (label c). This causes a rebound spike in G/S S T2 1 that activates interneuron I T1 1. This causes an inhibitory potential in G/S S T1 1 (label d) starting the cycle again. This process allows maintenance of spiking activity within this circuit. The spikes fall on alternating cycles of a theta rhythm oscillation shown at top (note that this theta rhythm oscillation is at double the frequency of the cycle skipping input to the interneurons).
the circuit maintains spiking activity with different neurons firing on alternate cycles of the network theta rhythm oscillation. This is consistent with data on theta cycle skipping shown in unit recordings from awake, behaving animals (Jeffery et al., 1995;Deshmukh et al., 2010;Brandon et al., 2013). Figure 11A illustrates that the spiking model with inhibitory feedback can hold a stable state for a period of time and also transition between neurons at a wide range of different frequencies. In the start of this simulation, stellate cell S T1 1 maintains its activity due to the mechanism shown in Figure 10. Because the timing of the rebound spikes matches the timing of the MS input to the interneurons. As the simulation progresses in Figure 11A, the magnitude of feedback inhibition from the interneurons gradually increases proportional to each time step. As the strength of feedback inhibition increases, the time of the rebound spike occurs at earlier phases relative to the medial septum input to interneurons. This causes a progressive shift in the spiking activity of different interneurons, which activate spiking in different stellate cells. Thus, the activity spreads from stellate cells S T1 1 and S T2 1 to the interneurons I T1 2 and I T2 2 to the stellate cells S T1 2 and S T2 2, and then on to S T1 3 and S T2 3 and eventually to S T1 5 and S T2 5. When the rebound spiking occurs earlier than S T1 5 and S T2 5, the next preceding phase is interneurons I T1 1 and I T2 1, so the cycle repeats itself through the network. This network model has multiple different interneurons receiving square wave rather than a sinusoid input from the medial septum, with no overlap between the phasic input to different interneurons. Square wave input helped maintain stability by preventing the excess spread of activity in the network, but stability should be feasible using sinusoidal input in a larger scale network. Excess activity was prevented in a different one-dimensional network model with a single interneuron and sinusoidal input to stellate cells (Hasselmo, 2013). That model has not yet been extended to two dimensions, but could be used to obtain grid cell activity using a small number of interneurons that fire at all spatial When the magnitude of feedback inhibition is progressively increased, this changes the speed of transition. Initially, the spiking activity stays stationary in stellate cells S T1 1 and S T2 1, but gradually as the magnitude of feedback inhibition increases the rebound spikes fall at earlier times relative to MS input. This results in a progressive shift to S T1 2 and S T2 2 and as the frequency difference gets larger the transition occurs more rapidly. Note that when S T1 5 and S T2 5 are active, the next earlier interneuron phase activates stellate cells S T1 1 and S T2 1, so that the sequence repeats through the line. (B) A continuous change in speed of transition can also be obtained by increasing the frequency of medial septal input. At the start, the frequency difference is large and the transition between neurons happens rapidly, but as the MS input frequency increases this eventually causes the rebound spiking to match the MS input resulting in stationary activity, and when MS input frequency is higher than the rebound spiking speed this causes a reversal in the direction of transitions between spiking neurons.
locations with different phases relative to sinusoidal input to stellate cells.
As the magnitude of feedback inhibition increases continuously in Figure 11A it causes a more rapid transition between different spiking neurons, corresponding to a higher speed of movement through the environment. This simulation demonstrates that a spiking network can show a shift in the location of activity due to rebound from inhibition as shown by the more abstract model in Figure 9. The speed of the transition through the network can be altered by changing the magnitude of feedback inhibition.
Alternately, the network activity can be also shifted at different speeds by changing the relative frequency of medial septal input, as shown in Figure 11B. If the frequency of medial septal input wave starts out low (left), this gives a period of the medial septal input wave that is slower than the interval of rebound spiking and the network activity shifts rapidly. As the medial septal input frequency is increased, this results in a smaller difference between the oscillation period and the rebound spiking interval, causing a slowing of the shift. The shift eventually becomes static (Figure 11B center) when the period matches the rebound spiking interval. Subsequently, as the period of the medial septal input waves becomes shorter than the rebound spiking interval, this causes movement in the reverse direction. Thus, the speed of movement in the network can also be regulated by the temporal frequency of medial septal traveling wave input. As can be seen in Figure 11B, this gives a continuous transition between different frequencies, as it changes the duration of the phases when each neuron can be activated. This predicts that some septal neurons should show speed modulation of firing frequency, and specifically might show decreases with increased speed. Preliminary data indicates speed modulation of septal neurons. The frequency only needs to shift by 3.5 Hz relative to baseline to cover the wide range of speeds shown in Figure 11B.

SIMULATIONS USING TWO-DIMENSIONAL TRAJECTORIES
The mechanism of shift between neurons based on rebound spiking has been expanded to a representation of two dimensional movement. The example shown in Figure 12, shows how the framework from Figures 10, 11 can work to cause shifts in cell firing in two different movement directions on different theta cycles. Stellate cell rebound spiking shown in Figure 12A2 interacts with traveling wave input from medial septum spreading through the conjunctive populations to activate interneuron populations in Figure 12A1. This results in movement of activity within the stellate cell population in different directions during different cycles.
The same activity is shown as part of two-dimensional plots of population activity in Figure 12B with synaptic input from the grid cell that activates a single column of conjunctive cell spiking to bring one conjunctive cell over threshold. Section A2 shows spiking in two rows of stellate cells. Note that the simulations used for Figures 12, 13 used four different conjunctive cell populations getting distinct directional traveling waves for each stellate cell population. In the figures only the population getting the West-moving traveling wave is shown as C T1,W and only the population getting the South-moving traveling wave is shown as C T2,S . These conjunctive cell populations then contact single interneuron populations that provide inhibitory input to each stellate cell population. In this manner, the simulation represents the effect of four different excitatory conjunctive cell populations getting the traveling wave input that thereby gate the activity arriving in two interneuron populations (so that there would be a much larger population of excitatory conjunctive cells than interneurons). Figure 12, inhibitory input to one stellate cell visible in both A2 and B S T1 causes a rebound spike at time 283 visible as a spike in the membrane potential in A2 and a black square in B S T1 . This activates a column of conjunctive cells in C T2,S that converges with the row of activity due to the Southmoving traveling wave at time 284. The convergence causes a conjunctive cell spike that activates inhibitory interneurons to cause inhibition (white square) in the corresponding stellate cell in S T2 at time 286. At time step 316, the stellate rebound spike activates a row of conjunctive cells in C T1,W that converges at an earlier time with the column of cells activated by the West-moving traveling wave, activating a conjunctive cell one step to the East. This activates inhibitory neurons to cause inhibition in S T1 at a position one step to the East. Similarly, the rebound spike at time 347 converges with the South-moving traveling wave in conjunctive cells (in C T2,S ) one step to the North at time 348. This causes The example in Figure 12 shows the mechanism for the twodimensional model to correctly shift the locus of spiking activity in proportion to movement in the two-dimensional environment.

At time 251 in
As shown in Figure 13A, this mechanism is reasonably accurate for a trajectory that scans across two-dimensional space with long segments in the north-south direction and only short movements in the eastward direction. Note that the internal model shifts in proportion to movement in the east-west direction and in response to the other movement relative to a reference vector at 60 • in the southwest to northeast direction corresponding to multiplying each location (x,y) by the operator: H = cos θ 1 sin θ 1 cos θ 2 sin θ 2 where θ 1 = 0 and θ 2 = π/3. This means that a 90 • angle difference in the actual trajectory is transformed to a 60 • angle in the internal representation, resulting in periodicity with the topology of a twisted torus, that causes a hexagonal periodicity. The accurate shift in spiking activity within the population can be seen by the plot of a single neuron from the center of the population that spikes in a regular hexagonal array of locations ( Figure 13A).
For certain parameters, the network can track more complex two-dimensional trajectories with frequent movement and with changes in speed. The accuracy of multiple directions of movement requires that rebound spiking be able to interact with the wave of medial septal input at large intervals. For example, a rebound spike in a stellate cell with index (1,1) might occur at the time of a medial septal input wave that causes an conjunctive cell to spike at location (1.5) and the next rebound spike may activate a conjunctive cell at location (4,5). The larger shift in neuron index in each direction is necessary to maintain accurate angles of movement. With this property, the network can track a trajectory with continuous curves, as shown in Figure 13B. In this figure, the actual trajectory is compared with an internal representation derived from the index of the stellate cell firing a spike on each time step. Shifting of the parameters allowed a reasonable match between the actual trajectory and the internal population activity in some cases (as shown in Figure 13B). As shown in Figure 13C, the accurate movements at a range of angles required larger steps in the index of conjunctive cells activated by the rebound spike of stellate cells combined with the medial septal input to conjunctive cells. Note that the examples in Figure 13C show very narrow waves of medial septal input, but it should be feasible to obtain similar function with smoother input waves from medial septum. This would allow future simulations to extend beyond single neuron spiking within the two-dimensional population, to obtain a delimited twodimensional region of spiking activity that shifts with running velocity.

FEEDBACK INTERACTIONS SELECT SMALL NUMBER OF TRAVELING WAVES
As noted in the Methods section, the interaction of waves of activity was also explored in a more abstract model with feedback interactions between grid cells and traveling wave inputs. The goal of this model was to understand how a large population of grid cells could interact with a large number of different head direction inputs. This more abstract model does not yet include rebound spiking, but could ultimately be combined with the model of rebound spiking. This more abstract model is similar to the model published in Hasselmo and Brandon (2012)  mathematical representation that is clearer to describe and more efficient to implement, allowing much more extensive analysis of different numbers of heading angles. Figure 14 shows the grid cell firing pattern resulting when a large number of different traveling waves are combined. The feedback interaction causes traveling waves that are not at 60 • intervals to be decreased in magnitude by the dynamics of the circuit, as shown in the figure. This allows effective generation of grid cells in the model. Different random initial conditions can result in different orientations of the hexagonal pattern across the population corresponding to different maximal activation of the heading angle waves providing input. Figure 15 shows that shifting the frequency of individual traveling waves can allow the activity to be moved around the network. Note that in this simulation, traveling waves can interact to generate a stationary grid pattern, and then be shifted if the relative frequency of traveling waves is shifted (Blair et al., 2013). Figure 15 shows that a simple trajectory through the environment that can generate a two-dimensional grid cell firing pattern. Figure 15 also shows that the hexagonal pattern of neuronal The firing activity of single neurons during movement is most reliable for 55 or more heading angles (right side). There is a lack of consistent firing for lower numbers of heading angles (e.g., 36 or 42) but the pattern is more consistent for larger numbers of heading angles. activity is only consistently obtained for large numbers of heading angles above 48 in Figure 15A and above 55 in Figure 15B.

DISCUSSION
The models presented in Figures 1-13 uses simulated neurons that replicate the resonance properties of stellate cells in MEC shown in extensive experimental studies of these neurons (Haas and White, 2002;Erchova et al., 2004;Giocomo et al., 2007;Engel et al., 2008;Fernandez and White, 2008;Giocomo and Hasselmo, 2009;Boehlen et al., 2010;Heys et al., 2010). The resonance allows rebound depolarization to mediate rhythmic interactions with conjunctive/pyramidal cells or interneurons receiving theta rhythmic input from the medial septum. The rhythmic interactions allow maintenance of oscillatory activity in subsets of neurons within the population. When the intrinsic resonance frequency of the stellate cells is increased relative to the medial septum input frequency, or the medial septum input frequency is decreased, this causes rebound depolarization to cause a shift in the activity state between different sets of neurons within the model. The movement along a two-dimensional trajectory can be obtained by skipping between shifts in different directions during alternating cycles of theta rhythm oscillations. The patterns of increases and decreases of activity in different subsets of neurons resembles the beat patterns in the oscillatory interference model of grid cells Blair et al., 2008;Burgess, 2008;Hasselmo, 2008) and related models (Nadasdy, 2009(Nadasdy, , 2010. But the use of coupled differential equations in Figures 1-9 allows a systematic change in amplitude of oscillations, rather than just changing the phase and frequency as in previous implementations of the oscillatory interference model. This provides a framework for describing how differences in resonance frequency at different dorsal to ventral positions Giocomo and Hasselmo, 2009;Boehlen et al., 2010;Pastoll et al., 2012a) can cause the difference in spacing of grid cell firing fields along the dorsal to ventral axis (Hafting et al., 2005;Sargolini et al., 2006;Stensola et al., 2012). This model also provides a potential framework for linking resonance at low frequencies in bat entorhinal cortex (Heys et al., 2013) to the generation of grid cells in bats (Yartsev et al., 2011;Barry et al., 2012a). As long as resonance strength is sufficiently strong, the model shows that the beat patterns providing spatial periodicity can be generated with resonance at lower frequencies similar to data recorded from the entorhinal cortex of two species of bat (Heys et al., 2013). Note that this model links resonance frequency to the deterministic timing of rebound depolarizations, but does not require the noise-driven subthreshold oscillations of entorhinal stellate cells (Dickson et al., 2000). This is an important advantage for this model, as subthreshold membrane potential oscillations do not have reliable time course (Zilli et al., 2009) and cannot show independent frequency and phase within a neuron (Remme et al., 2010). The use of rebound depolarizations in this model resembles the mechanism of maintaining activity patterns in a model of short-term memory (Lisman and Idiart, 1995) and to encode sequences for modeling theta phase precession (Jensen and Lisman, 1996). The same mechanism was also used to maintain an attractor in a previous paper (Navratilova et al., 2012). The mechanism of attractor maintenance is similar in focusing on rebound depolarization, but that previous paper triggered the intrinsic mechanisms with a spike in an integrate and fire model (Navratilova et al., 2012), rather than due to subthreshold rebound from inhibition. Also, the attractors were more directly shifted by a velocity input, rather than being indirectly shifted by the effects of velocity on the relative phase of interaction of resonant cells with cells receiving MS input. In the model presented here, the loss of MS input prevents grid cell periodicity, as supported by data showing loss of spatial periodicity of grid cells with inactivation of the medial septum (Brandon et al., 2011;Koenig et al., 2011). The previous model does not depend upon MS input, and must also attribute periodic boundary conditions to patterns of synaptic connectivity rather than the cyclical nature of the phase of oscillatory input from the medial septum as used here.
Another previous model generated grid cell firing patterns based on waves of activity moving across the entorhinal cortical sheet (Nadasdy, 2009(Nadasdy, , 2010 to allow coding of sensory input by spike phase. Decoding of this phase signal relative to repeated waves in two dimensions resulted in hexagonal periodicity. That model predicted a hexagonal distribution of grid cell coding across the cortical surface, as supported by a recent paper (Ray et al., 2014). That model did not use intrinsic rebound spiking properties or traveling waves oriented with multiple head directions.
This current model addresses the increase in spacing observed in mice with the forebrain restricted knockout of the HCN1 subunit of the h current  which alters the parameters of the h current so that cells do not show resonance at theta frequency, but instead show resonance at lower frequencies ). Here this difference was simulated by altering the parameter g h in Equation (1) and showing that resonance could still mediate grid cell periodicity. If the lower frequency resonance shows smaller frequency shifts relative to velocity, this would cause broader spacing between grid cell firing fields. Thus, the model is supported by experimental data showing larger spacing between firing fields in more ventral entorhinal cortex (Hafting et al., 2005) and lower resonance frequency of ventral stellate cells , and by data showing that knockout of the HCN1 subunit of the h current causes lower resonance frequency ) and wider spacing of firing fields . The model predicts a total loss of grid cell periodicity with total loss of rebound spiking if h current were completely removed by loss of both HCN1 and HCN2 subunits. The two-dimensional model properties (Figures 12, 13) are supported by recent data showing lack of shared head direction tuning for neurons that fire asynchronously on alternating cycles of theta rhythm (Brandon et al., 2013), and predicts that there should be coupled pairs of asynchronous neurons that fire on opposite cycles but share the same location of their firing fields, while coding different directions of movement and changing their firing pattern on different cycles. In addition to support from data showing loss of grid cell periodicity with medial septum inactivation (Brandon et al., 2011), the model is supported by data (Welday et al., 2011) showing cosine-tuned directional modulation of septal neurons that could be part of ring attractors generating traveling waves. The model also predicts that septal neurons should show speed modulation of firing frequency, and that sets of septal cells should show shared orientation tuning but different temporal theta phase offsets consistent with traveling waves. Note that the model is consistent with evidence of resonance in neocortical and hippocampal circuits (Stark et al., 2013).
The model presented here also provides a potential framework for addressing the effects of acetylcholine on the resonance properties of neurons (Heys et al., 2010;Barry et al., 2012c) and the theorized role of this effect in underlying the experimentally demonstrated change in spacing of grid cells in novel environments (Barry et al., 2012b). A change in acetylcholine levels could cause a shift in overall speed modulation (Newman et al., , 2014, or could be directly involved in modulating intrinsic resonance frequency relative to running speed. Acetylcholine levels have been shown to change in running vs. stationary rats (Mizuno et al., 1991;Marrosu et al., 1995) and have been shown to change the resonance frequency of medial entorhinal neurons (Heys et al., 2010;Heys and Hasselmo, 2012). Alternately, resonance frequency could be shifted by systematic changes in the resonance frequency of neurons due to different levels of depolarization with running speed. Experimental data shows that the resonance frequency decreases with depolarization in MEC (Shay et al., 2012). The model presented here has a potential advantage over attractor models in that it can maintain activity with low firing rates at theta frequency. Traditional attractor models use slow excitatory synaptic currents to maintain activity over long delay intervals between theta cycles. For example, the models by Pastoll et al. (2012b) uses NMDA currents with time constants of 100 ms, which is longer than the time constant of excitatory synaptic potentials demonstrated in the entorhinal cortex (Garden et al., 2008).
The prominence of head direction sensitivity in driving MEC neurons (i.e., both head direction cells and conjunctive cells) suggests the importance of sensory input for updating grid cells, requiring constant monitoring of the orientation of sensory input to the head. The spherical coordinates of each sensory stimulus relative to the animal could influence the phase and amplitude of specific traveling waves, and this could be included in the summation of traveling waves described in the model presented in Figures 14, 15.
In this framework, in addition to integrating body movement, the system may involve learning a link between sensory features and particular phases of traveling waves. This influence on waves could include visual features, but also input from a rodent's whiskers, and auditory input. The influence on waves could involve an interaction of head direction with the relative bilateral sensory input in the two eyes or two whisker systems. A rotational shift in sensory angle in the same direction on both sides of the head corresponds to pure head movement without body movement and should not alter the phase of waves. But movement with the head direction held sideways relative to movement should shift the phase of waves lateral to head direction. If sensory flow is in opposite directions for the two eyes for a given head direction during forward movement then this should shift the phase of waves parallel to head direction. In addition to addressing the predominance of head direction input, a framework driven primarily by sensory input could address the change in firing patterns such as the expansion or compression of grid field spacing with changes in environment size Stensola et al., 2012) and the linearization of grid cell firing in the hairpin maze (Derdikman et al., 2009). To allow sufficient influence of input on the attractor state it may be necessary to change the relative strength of input vs. feedback on different phases of theta rhythm (Hasselmo et al., 2002).
In simulating grid cell firing properties, a modeler does not need to choose between oscillatory interference and attractor dynamics, as these are not incompatible models, as shown by Bush and Burgess (2014). The model presented here uses elements of oscillatory interference, but the effects of excitatory or inhibitory feedback could be interpreted as resulting in dynamical attractors within the model. The model shows that resonance frequency indicates the time course of rebound depolarization that allows generation of beat patterns by interaction of resonance frequency with the frequency of input from the medial septum. The model shows that with sufficient resonance strength, beat patterns can be generated with a range of different resonance frequencies. Figures 12, 13 show that the one-dimensional firing patterns generated by rebound oscillations can be combined into two-dimensional grid cell firing patterns as proposed in previous models of grid cell firing patterns (Burgess, 2008;Navratilova et al., 2012).