## Abstract

One of the two primary classes of models of grid cell spatial firing uses interference between oscillators at dynamically modulated frequencies. Generally, these models are presented in terms of idealized oscillators (modeled as sinusoids), which differ from biological oscillators in multiple important ways. Here we show that two more realistic, noisy neural models (Izhikevich's simple model and a biophysical model of an entorhinal cortex stellate cell) can be successfully used as oscillators in a model of this type. When additive noise is included in the models such that uncoupled or sparsely coupled cells show realistic interspike interval variance, both synaptic and gap-junction coupling can synchronize networks of cells to produce comparatively less variable network-level oscillations. We show that the frequency of these oscillatory networks can be controlled sufficiently well to produce stable grid cell spatial firing on the order of at least 2–5 min, despite the high noise level. Our results suggest that the basic principles of oscillatory interference models work with more realistic models of noisy neurons. Nevertheless, a number of simplifications were still made and future work should examine increasingly realistic models.

## Introduction

Following the initial discovery and characterization of grid cells in rat entorhinal cortex (Fyhn et al., 2004; Hafting et al., 2005; Sargolini et al., 2006), a number of models of the spatial firing properties of these cells were offered. Generally, the mechanisms of the models fall into two categories: continuous attractors (Fuhs and Touretzky, 2006; McNaughton et al., 2006; Guanella et al., 2007; Burak and Fiete, 2009) and oscillatory interference (O'Keefe and Burgess, 2005; Blair et al., 2007, 2008; Burgess et al., 2007, 2008; Hasselmo et al., 2007; Hasselmo, 2008) (but see Gaussier et al., 2007; Kropff and Treves, 2008).

Oscillatory interference models were first used to explain hippocampal place cell phase precession (O'Keefe and Recce, 1993; Lengyel et al., 2003; Huhn et al., 2005) by combining two oscillators to produce an interference pattern of activity in the spiking of a neuron. Grid cell models simply add one or more additional oscillators to produce a two-dimensional (2D) interference pattern. One oscillator maintains an arbitrary baseline frequency and the remaining, active oscillators are driven to various frequencies above or below the baseline frequency. The specific frequencies used at each time are based on the animal's velocity in different directions through a simple transformation that results in the phase difference between an active oscillator and the baseline oscillator encoding one-dimensional (1D) positional information. Active oscillators maintaining 1D positional information along directions 60° apart produce a regular hexagonal interference pattern in space.

One prominent criticism of these models (Giocomo and Hasselmo, 2008; Hasselmo, 2008; Welinder et al., 2008; Burak and Fiete, 2009; Zilli et al., 2009) is that they are presented in terms of abstract, perfect oscillators, whereas oscillators in the brain are noisy and have more complicated dynamics. Zilli et al. (2009) showed that experimentally measured variability in the examined biological oscillators was large enough that an oscillatory interference model's spatial firing would be expected to remain stable for a few seconds at best. This strongly supported the noise criticism, but there are at least two ways for these models to overcome this problem. One possibility is that sensory cues can frequently or even constantly reset the grid network (Redish and Touretzky, 1997; Redish, 1999; Burgess et al., 2007; Samu et al., 2009). A second possibility, as suggested by Zilli et al. (2009), is that these individually noisy oscillators may be coupled *in vivo*, and through synchronization the network's oscillations may be less variable (Manor et al., 1997; Needleman et al., 2001; Ly and Ermentrout, 2010; Tabareau et al., 2010).

To test the latter solution, we use numerical simulations of two different neural models and verify that coupled, noisy neural oscillators can be used successfully in the oscillatory interference framework. We close with discussion of our assumptions as well as implications and directions for future study.

## Materials and Methods

##### Computational methods.

Simulations were performed in MATLAB 7.0.0.19920 (R14). MATLAB source code to reproduce all figures and to support many claims labeled preliminary simulations or unpublished observations is available on ModelDB. Equations of the models were solved using the forward Euler method with a time step of *dt* = 0.1 ms (simple model) or *dt* = 0.01 ms (biophysical model). Of these two models, only one was used to model the velocity-controlled oscillators (VCOs) in any particular simulation. Preliminary simulations showed these step sizes were sufficiently small. As our simulations often comprised hundreds of cells (15,000 in the largest simulation) and were up to 320 s long, it was important for processing time that the step size not be too small. The grid cell itself was modeled as a leaky integrate-and-fire (LIF), resonate-and-fire, or simple model neuron, and it was simulated alongside either the simple model or the biophysical model VCOs, using the same size time step the oscillator model used.

Spike times were determined by comparing the voltage variable with a fixed value (thresholds of 1 for the LIF model, 1 for the resonate-and-fire model, θ_{spike} for the biophysical model, and a peak value *v*_{peak} for the simple model). For all models, a spike was recorded on the time step where the voltage crossed the threshold from below (and for the LIF, resonate-and-fire, and simple models, the voltage was then immediately set to the reset voltage, as specified below).

To provide two-dimensional trajectories as input for our simulations, we used experimentally collected rat trajectory data from Hafting et al. (2005) (available for download at http://www.ntnu.no/cbm/moser/gridcell). The trajectory data are a set of coordinates, *x*_{exp}(*t*) and *y*_{exp}(*t*), sampled at resolution *dt*_{exp} = 0.02 s (note that some trajectory files seem to contain multiple, concatenated trajectories separated by a discontinuity). The difference between adjacent position samples was used as the velocity input to the simulation. Simulations were performed at a finer temporal resolution than *dt*_{exp}, so the velocity signals were linearly interpolated, producing velocity signals *v*_{x}(*t*) and *v*_{y}(*t*).

To adjust the smoothness of the trajectory, *v*_{x}(*t*) and *v*_{y}(*t*) were bidirectionally (in time) low-pass filtered with a third-order Butterworth filter at 0.4 Hz (an arbitrary value controlling the smoothness of the trajectory). Most of our 2D simulations used the filtered velocity input, as the unfiltered velocity results were often poorer in comparison. Which of the two cases is more realistic is unclear, as the experimental trajectory recording likely contains high-frequency jitter from the video tracking system and the vestibular and proprioceptive systems of the animal themselves may also act to low-pass filter their acceleration or velocity signals to prevent small movements of the head from being path integrated. From the trajectory, we can calculate two spatial functions that are of frequent interest: speed, *s*(*t*) = *t*) = atan[*v _{y}*(

*t*)/

*v*(

_{x}*t*)], where atan is the four quadrant arctangent. Filtering decreases the velocity of the trajectory. For example, in 320 s of trajectory prepared in this manner, the mean (peak) instantaneous speed is 22 cm/s (89 cm/s) in the unfiltered trajectory compared with 17 cm/s (49 cm/s) in the filtered trajectory.

Our simulations required knowing the frequency response, *F*(*I*), of the oscillator cell or network to current injections. In most cases, *F*(*I*) was calculated numerically using input magnitudes at fixed intervals to produce a resolution of ∼0.02 Hz over ranges of at least 4 Hz and other values were linearly interpolated. For the *n* = 5000 case, a resolution of ∼0.15 Hz was used. A range of 4 Hz was needed because we used β = 2 Hz/(m/s) (see below) and allowed for a maximum instantaneous velocity of 1 m/s, thus requiring 2 Hz above and below the baseline frequency. When the networks comprised noisy neurons, the measured *F*(*I*) curve became noisy, which imposed a maximum level of accuracy (Fig. 1*A*, black line, uncoupled noisy cell).

We also simulated abstract VCOs Φ_{i} and grid cells using a common form of the oscillatory interference model (Eq. 1). The abstract model was simulated using the forward Euler method using the same time resolution as our network model (which depended on the neural model in use). Each abstract VCO's state was characterized by its phase ϕ_{i} evolving at a time-varying frequency, *f*_{i}(*t*) = ω_{b} + β*s*(*t*)cos[ϕ* _{i}* − ϕ(

*t*)], which was a function of speed,

*s*(

*t*); body direction ϕ(

*t*); and each VCO's preferred direction ϕ

_{i}.

A quantitative measure of the phase of the network oscillators was desired, but translating the state of a neuron into the corresponding phase is not a trivial task and translating the state of a network of coupled neurons into a population phase is harder still. Generally, each point in the neuron's state space can be identified with an asymptotic phase (Izhikevich, 2007), but finding the phase of an arbitrary point is computationally intensive. Instead of identifying phases in this manner, we simulated an abstract VCO Φ_{i} alongside each network VCO *V*_{i}. The frequency of Φ_{i} was set on each time step to *f*_{i}(*t*), as given above. In theory, if *V*_{i} is also at frequency *f*_{i}(*t*) at time *t* and if *V*_{i} and Φ_{i} are at phase 0 at time 0, then Φ_{i} and *V*_{i} will always be at the same phase at any time and we can use the phase of Φ_{i} as a measure of the phase of *V*_{i}. In reality, *V*_{i}'s frequency will not be perfectly controllable, so the difference in phase between Φ_{i} and *V*_{i} is a measure of the phase error *V*_{i} has accumulated. We recorded this error (the difference in phases) each time any cell in *V*_{i} emitted a spike.

Inaccuracies in *F*(*I*) can be particularly noticeable in the phase differences between baseline oscillators because these both maintain a constant frequency, which should be identical. If a network oscillator's frequency is slightly different from the abstract VCO's frequency due to *F*(*I*) interpolation, the phase difference will show an apparent linear drifting error that does not actually affect the network model. We avoided this cosmetic problem by specifically selecting the baseline frequency as one of the points of the measured *F*(*I*) curve so that no interpolation was needed. This ensured the abstract and network oscillators were truly operating at identical frequencies at baseline.

##### Model description.

Our model is summarized in Tables S1–S7 (parameters in Table 1) following the good model description practices suggested by Nordlie et al. (2009). The architecture of the model is summarized in Figure 2.

Our network oscillatory interference model is composed of a single cell *G* (the grid cell itself), which receives input from one or more (generally three) oscillatory networks *V*_{i}, 0 ≤ *i* ≤ *n*_{VCO}, which are identical apart from their inputs. *V*_{0} will often be called the baseline network oscillator.

Cells within each network *V*_{i} are recurrently coupled all-to-all (no self-connections) by identical synapses or gap junctions of strength *g* (except in one simulation where the connectivity probability is *p* = 0.01 and all *V*_{i} use the same connectivity matrix for reasons of computational efficiency). The cells in each *V*_{i} are the only cells that receive external input or noise. All cells in each *V*_{i} also project onto *G* with identical synapses of strength *w*_{i}. When synapses from all networks have the same strength (i.e., *w*_{i} = *w*_{j} for all *i*, *j*), the weight is simply referred to as *w* for conciseness. Importantly, the networks *V*_{i} are not coupled to each other, to prevent them from synchronizing with each other (a problem for single-cell oscillatory interference models) (Remme et al., 2009).

*G* is generally modeled as a leaky integrate-and-fire neuron characterized by its membrane time constant τ. The threshold of the cell is 1, and after a spike its voltage is reset to 0. We also considered a modification where oscillator *V*_{0} is assumed to have a privileged connection onto *G* such that inputs from any other oscillator *V*_{i}, *i* ≠ 0, have no effect unless they occur in a *t*_{gate} second window after any spike from a cell in *V*_{0}. Thus, a gate that allows inputs from the other *V*_{i} networks to reach *G*'s soma is opened only briefly after the baseline oscillator fires (Ang et al., 2005; Jarsky et al., 2005).

Alternatively, we sometimes modeled *G* as a resonate-and-fire neuron (Izhikevich, 2001). The model is the simplest model of resonance: the subthreshold dynamics are a two-dimensional linear time-invariant system with a pair of complex conjugate eigenvalues, *c*_{res} ± ω_{res} *c*_{res} < 0) and with resonant frequency ω_{res}. The model can be interpreted as a linearization of a resonant biophysical model (Izhikevich, 2001). The 2D state of the neuron can be written as a single complex number. The real part of the neuron's state is thresholded as a spiking mechanism. When the real part of the state crosses from <1 to >1, a spike is emitted and the state of the neuron is reset to

Finally, for the simulation using inhibitory synapses onto *G*, we modeled *G* as a simple model (Izhikevich, 2003, 2007) neuron that receives a current injection to fire continuously. The dynamics of the simple model describe a neuron's action potential upstroke and its subthreshold and postspike behavior. The action potential downstroke is modeled explicitly by resetting the membrane potential to a specified value when it reaches a peak value (the peak of an action potential, not a firing threshold). After each spike, the recovery variable is incremented by a fixed amount. This increment represents the build-up of slow currents during the action potential. Our implementation of this model was based on the MATLAB code given in Izhikevich (2007).

Networks *V*_{i} are usually modeled as sets of identical simple model neurons with an extra additive voltage noise term. Since discontinuities (like the action potential downstroke) can affect synchronization properties (Teramae and Tanaka, 2004) and other approximations made in the simple model may affect our results, we repeated some of the simulations with a biophysical model.

The biophysical model of an entorhinal cortical layer II stellate cell uses parameters from Acker et al. (2003) but with an additive noise term (instead of their stochastic persistent sodium channel model) and delta-current synapses (instead of their AMPA synapse model). This model of entorhinal cortex layer II stellate cells contains standard spiking sodium and potassium channels as well as a noninactivating persistent sodium current and a two-component (fast and slow) H current (Fransén et al., 2004) thought to play a large role in the stellate cell subthreshold dynamics (Magistretti and Alonso, 1999, 2002; Dickson et al., 2000). We decreased the conductance densities in the model (by increasing *C*_{m}) to slow the firing frequencies such that the model could fire at 5 Hz, as needed to match biological noise levels as described below.

All cells in network *V*_{i} receive the same scalar input signal *I*_{Vel,i}(*t*) as a current injection (each network has its own input signal). The magnitude of the input at any time is calculated in terms of the desired population frequency of *V*_{i}. Given a desired frequency *f*_{i}, defined above, and the frequency response of the network to a given input, *F*(*I*_{Vel,i}), the input to all cells in network *V*_{i} is *I*_{Vel,i}(*t*) = *F*^{−1}(*f*_{i}).

The primary output measure of the model is the firing of cell *G*. When a single network *V*_{0} drives *G* (using *G* as an indirect measure of the behavior of *V*_{0}), the mean μ and SD σ of the interspike intervals (ISIs) of *G* are measured. In case *G* fires a burst of spikes in response to sustained inputs, ISIs <50 ms long are collapsed into one long ISI, which represents the interval between first spikes of each burst (e.g., consecutive ISIs of length 100, 7, and 90 ms become a pair of ISIs of 100 and 97 ms). This easily allows synchronous population activity in the VCO to be measured. Henceforth, when the mean or variance of the periods of a VCO zzare referred to, we mean the corresponding statistic of the spikes emitted by *G*. We can use these statistics to evaluate how well oscillators will perform in the oscillatory interference model. Zilli et al. (2009) derived an expression for the estimated stability of a grid pattern as a function of the mean and variance of the periods of a noisy oscillator. Using the approximation that noise results in oscillator period durations independently drawn from a normal distribution with mean μ and SD σ (in seconds), the estimated stability time formula *t* = 5μ^{3}/(4πσ)^{2} gives the number of seconds until a grid cell is no longer likely encoding the correct spatial location. This is derived by noting that the difference in lengths of periods of two identical but noisy oscillators at the same frequency has twice the variance 2σ^{2} of the length of a period of either oscillator, and the difference in lengths of *n* consecutive (independent) periods has the variance 2*n*σ^{2}. The variance is converted from units of seconds^{2} to radians^{2} (2*n*σ^{2}(2π)^{2}/μ^{2} radians^{2}) and the estimated stability time is the amount of time for the variance to hit a threshold of 2.5 radians^{2} (for further details, see Zilli et al., 2009).

When two or more networks *V*_{i} drive *G*, the inputs to *V*_{i} represent a spatial trajectory, and a primary measure is the set of spatial positions where *G* emits spikes. The spatial autocorrelation of the firing positions was calculated numerically in MATLAB by binning the spike positions using the hist3 function, then calculating the autocorrelation by taking the two-dimensional convolution of the histogram with itself (using the conv2 function with the same option). A final measure commonly used is the phase difference between a network oscillator *V*_{i} and the corresponding abstract, idealized oscillator Φ_{i}, as described above.

The parameters of the model are summarized in Table 1, where default values are given (where a simulation uses a different value, it is given in the text or figure caption). The parameters of the simple model neurons (*C*, *v*_{r}, *v*_{t}, *v*_{peak}, *a*, *b*, *c*, *d*, *k*) correspond to a resonant regular spiking neuron (e.g., a cortical pyramidal cell), as given in Izhikevich (2007) (the original model had *b* < 0, which we negated to get a resonant model). The biophysical model parameters were taken from Acker et al. (2003), with *C*_{m} increased as shorthand for proportionally decreasing all ionic conductances. The noise level σ was selected to phenomenologically match experimental data (supplemental material). As the noise σ increases, the number of cells *n* may need to increase and coupling strength *g* may need to change to maintain synchronization. The baseline frequency ω_{b} was set to a frequency near the middle of the measured *F*(*I*) curve at a current value where the exact frequency is known. For the LIF model, the time constant τ and synaptic strength *w* (and *t*_{base} when used) were adjusted to ensure all VCOs fired at nearly the same time, which is need in order for *G* to fire, and there was a range of τ and *w* that produced similar results (possibly changing the size of the fields themselves, up to the point of firing occurring at all or no locations). For the resonate-and-fire model, the decay rate *c*_{res} was selected similarly, and the resonant frequency was ω_{res} = ω_{b}. Note that when *c*_{res} ≪ 0, the cell is essentially nonresonant, so we kept *c*_{res} ≈ 0 to focus on the case where resonance may play a large role. The velocity-to-frequency slope β controls grid spacing as well as the range of frequencies of the *F*(*I*) curve that must be calculated. We used β = 2 Hz/(m/s), the mean value of β for theta frequency oscillations in entorhinal cortex layer II stellate cells (Jeewajee et al., 2008).

Initial suggestions linked grid cells to the membrane resonance of entorhinal cortex layer II stellate cells (Giocomo et al., 2007). Both our biophysical and simple model parameterizations are also resonant; however, as no current evidence specifically implicates entorhinal cortex layer II stellate cells (or resonance in general) in the generation of the grid pattern, preliminary simulations of many of these results were also tested with the networks *V*_{i} comprising simple model neurons parameterized as integrators instead of resonators (setting *b* = −2) and all tested results were reproducible with that model. Generally, both types of neurons are found in all regions. Resonant cells in the medial temporal lobe also include, for instance, CA1 pyramidals (Leung and Yu, 1998; Hu et al., 2002, 2009) and interneurons (Chapman and Lacaille, 1999), cells in lateral and basolateral amygdala (Pape et al.,1998), and >80% of layer II parasubicular neurons (including principal stellate and pyramidal neurons in addition to interneurons) (Glasgow and Chapman, 2007). Cells reported to be nonresonant (integrators) include nonstellate cells in entorhinal cortex layer II (Alonso and Klink, 1993; Haas and White, 2002; Erchova et al., 2004) and pyramidal cells in entorhinal cortex layer III (Erchova et al., 2004). However, conclusions regarding resonance drawn from *in vitro* experiments may not always correspond to *in vivo* behavior (Fernandez and White, 2008; Prescott et al., 2008).

## Results

### Oscillatory interference

Our oscillatory interference model was composed of one baseline oscillator *V*_{0}, which maintained a constant frequency and a number of active oscillators *V*_{i}, i ≥ 1, which changed their phase relative to the baseline oscillator to encode positional information. The simulated animal's body velocity input changed the oscillators' frequencies slightly above or below the baseline frequency such that the set of locations where any single active oscillator was in phase with the baseline oscillator formed a set of parallel bands in the environment (though these bands do not appear as spiking in this model).

The function of the model is easily demonstrated using the abstract version (Burgess et al., 2007), where the VCO phases evolve according to
and
We define *f*_{i}(*t*) = ω_{b} + β*s*(*t*)cos[ϕ_{i} − ϕ(*t*)], the desired frequency of oscillator *i* at time *t*, as a function of velocity at time *t* (*s* and ϕ) and preferred direction ϕ_{i}. Notice in Equations 1 and 2 that the phase difference ϕ_{i} − ϕ_{0} changes at frequency β*s*(*t*)cos[ϕ_{i} − ϕ(*t*)] Hz. If the animal moves in direction ϕ_{i} at *s*(*t*) = 1 m/s and if β = 1 Hz/(m/s), the phase difference will oscillate at frequency β*s*(*t*) = 1 Hz. If the phase difference is initially 0, the phase difference will return to zero each 1 m the animal walks (like a digit in an odometer rolling over). If the animal were to move at *s*(*t*) = 2 m/s then the phase difference would oscillate at β*s*(*t*) = 2 Hz and so the phase difference will still return to zero each 1 m the animal walks. In fact, this is true for any 1D trajectory *s*(*t*). If β = 2 Hz/(m/s), then the phase difference will return to zero each 1/2 m (generally each 1/β m), so β acts to set the distance in space over which the phase difference repeats.

The abstract grid cell is then said to fire if Σ_{1≤i≤nVCO}(cosϕ_{0} + cosϕ_{i}) exceeds a threshold θ. Each term in the sum represents the readout via an interference mechanism of the respective phase differences, and the model fires when all of these phase differences are near 0. When the directions ϕ_{i} are separated by 60° increments (and are not all colinear), the phase differences all equal 0 at positions that are hexagonally arrayed.

Three neighboring fields form an equilateral triangle where the distance from the base of the triangle to the opposite vertex is 1/β m. By trigonometry, the distance between neighboring fields is then λ = 2/(√3β) m. This model is thoroughly described in a number of previous publications (Burgess et al., 2007; Giocomo et al., 2007; Burgess, 2008; Hasselmo, 2008) and further details can be found there.

### Single cell oscillators

In Figure 3, noiseless simple model neurons are used as single-cell VCOs (Burgess, 2008). For computational efficiency, this manuscript focuses on the case of two active VCOs at angles 0° and 120°, although preliminary simulations suggest our model is successful when *n*_{VCO} = 4 (three active VCOs at angles 0°, 120°, and 240°). On each time step of the simulation, the instantaneous velocity inputed along the directions 0° and 120° set the frequency of the active oscillators (via the *F*(*I*) curve) (Fig. 1*A*) so that the phase difference changed to represent the distance traveled during that time step. At the positions where the oscillators were all in phase, the postsynaptic cell *G* was receiving maximal input.

In Figure 3, input from the baseline oscillator neuron *V*_{0} is made particularly strong and the inputs from the two active oscillators *V*_{1} and *V*_{2} are scaled so that a single input is subthreshold and both inputs must occur shortly after the baseline input to push *G* over threshold. If all synaptic weights were equal, *G* would produce either tiny or asterisk-shaped fields (Burgess, 2008, his Figs. 7 and 8) with lines extending along directions where subsets of VCOs are in phase. Burgess (2008) used a multiplicative baseline sinusoid to solve this problem, and our comparatively stronger baseline input is an alternative solution. With the strong baseline, the firing fields of the grids become triangular in shape, as seen in Figure 3, *A* and *B*.

Despite the unusual shape, the hexagonal spatial pattern compares well to the output of the abstract version of the model shown in Figure 3, *C* and *D*. The phase error between respective neural and abstract oscillators, a quantitative measure of performance, is plotted in Figure 3, *E–G*. The near-flatness of the lines shows that the neural model here is performing almost perfectly. The slight fluctuations in phase errors of a VCO are closely related to its velocity input, as shown by the superimposed velocity signal in Figure 3*E*. These fluctuations are much larger when the input is an unsmoothed spatial trajectory (supplemental Fig. S1, available at www.jneurosci.org as supplemental material).

The similarity of the phase difference error to the velocity input suggests a relationship between changes in the input and changes in the phase error. In fact, it is a general property of neurons that changes in input produce transient, history-dependent changes in the *F*(*I*) relation (supplemental text and supplemental Figs. S2, S3), which we refer to as adaptation. For example, changing *I* to a level at which the neuron should fire at 7 Hz actually results in a brief period of slightly higher or lower frequency, which will introduce a small phase error. This error can be minimized by reducing β to reduce the range of VCO frequencies. However, decreasing β will also increase the relative impact of noise-related VCO frequency fluctuations. Additionally, the grid spacing λ = 2/(√3β) is a function only of β, as shown above, so the spacing cannot be held constant as β is varied, meaning these opposing pressures on β directly correspond to pressures on grid spacing.

On the single cell level, however, the effects of noise are far greater than the effects of adaptation. A simulation using single, noisy cells as oscillators (Fig. 4) completely failed to produce the grid pattern. The intrinsic noise level used in this and subsequent simulations was selected (see supplemental material) so that an uncoupled simulated cell had the same ISI mean and variance as the persistent spiking cell with the median estimated stability time from Zilli et al. (2009). An example ISI histogram of a simple model cell with this noise level is shown in Figure 5*A*. From the phase errors in Figure 4, *E–G*, it is clear that the oscillators lose their correct phase relative to the baseline almost immediately. This result was expected on the basis of previous noisy simulations (Burgess et al., 2007; Giocomo and Hasselmo, 2008) and is consistent with the estimated time scale for the level of noise used (Zilli et al., 2009).

### Network oscillators

Because coupling can reduce the period variance of noisy oscillators as a network (Needleman et al., 2001) and individually (Ly and Ermentrout, 2010), we investigated whether networks of coupled, noisy oscillators could play the role of VCOs. For both synaptic and gap-junction coupling, we found that *n* = 250 all-to-all connected, noisy, simple model cells (coupling strengths *g* = 256 and *g* = 0.1, respectively) produced synchronized populations with acceptably low period variance of the network as a whole and the individual units (for the effects of varying important parameters, see supplemental material and supplemental Fig. S4, available at www.jneurosci.org as supplemental material). We can measure the *F*(*I*) curves of these networks (Fig. 1), then drive the network to fire at any desired frequency, just as with the single cell VCOs.

Three synaptically coupled networks, one acting as a baseline oscillator, together produce the simulation output shown in Figure 6. A voltage trace of one neuron from each *V*_{i} is shown in Figure 7, along with the activity of the gated LIF cell simulated in Figure 6, *A* and *B*, and a resonate-and-fire cell shown in Figure 6, *C* and *D*. Both types of postsynaptic cell produced reasonable hexagonally arrayed firing. The blurring directed along the preferred direction of *V*_{2} is consistent with the phase errors in *V*_{1} (Fig. 6*E*) but not *V*_{0} and *V*_{2} (Fig. 6*F*,*G*) as expected from previous analysis (Zilli et al., 2009). This simulation using synaptically coupled cells is similar to that of Hasselmo (2008), where the persistent spiking input was assumed to come in population bursts, shown here in Figure 7. Here the bursts are caused by strong synaptic coupling, as the neural model does not burst when uncoupled (supplemental Fig. S5, available at www.jneurosci.org as supplemental material) or when less tightly synchronized (Fig. 8*A*).

Instead of using synaptic connections, successful performance also occurs if the cells in the VCOs are coupled through gap junctions. Supplemental Figure S6 (available at www.jneurosci.org as supplemental material) shows the same simulation but with gap-junction coupling instead of synaptic coupling, and a voltage trace of one cell from each *V*_{i} and *G* are plotted in supplemental Figure S7 (available at www.jneurosci.org as supplemental material). Achieving large fields proved difficult with the resonate-and-fire model with a low damping constant, but was successful with the gated, integrate-and-fire cell. When gap-junction coupled, the cells did not burst and the smaller range of spiking phases (compared with synaptically coupled cells) may have a role in the increased difficulty we found in achieving larger fields.

The lack of excitatory recurrent connections in entorhinal cortex layer II (Dhillon and Jones, 2000; Couey and Witter, 2010) suggests that the oscillators may be inhibitory interneurons (if the oscillators are located there). To test this theory, we modeled the *V*_{i} neurons as gap-junction-coupled inhibitory neurons and *G* as a simple model neuron driven to fire at the baseline frequency. The inhibition from *V*_{i} onto *G* is modeled as an exponentially decaying input. As shown in Figure 9, when the VCOs are out of phase, *G* is tonically inhibited, but when the VCOs are in phase (in a grid field) there is a long enough window between the inhibitory inputs that the cell can fire, and it produces the characteristic field spacing as a result.

Though the simple model is an excellent approximation of the subthreshold dynamics of biophysical models, it is worth verifying our results with a relatively more realistic model. In Figure 10, we show the results of a simulation where the VCO networks are now composed of noisy biophysical neurons (with currents *I*_{Na}, *I*_{K}, *I*_{NaP}, *I*_{H fast}, and *I*_{H slow}). The output of this simulation is at least as good as the previous simulations.

In physiological experiments, it is usually only feasible to record from very few cells at any time, so it is worth considering how the ISI statistics of individual cells in a synchronized VCO network relate to the ISI statistics of the network itself in these simulations. In Figure 5, ISI histograms are shown for an uncoupled cell compared with cells from coupled networks and the variability-reducing effect of synchronization is particularly clear (Fig. 8). The simulated histogram in Figure 5*A* is comparable to the experimental subthreshold oscillation period histograms in Giocomo and Hasselmo's (2008) supplemental Figure S5*A* (available at www.jneurosci.org as supplemental material). For synchronized gap-junction-coupled neurons, Figure 5, *B* and *E*, individual neurons had essentially identical statistics to the network as a whole. When neurons were synaptically connected, individual neurons also tended to have similar periods to the network as a whole, but individual cells had a higher variability (larger ISI SD) (Fig. 5*C*,*F*). However, if the level of connectivity in the network is much lower, the network as a whole can produce low-variance oscillations even though the individual cells in the network show considerable ISI variance (Fig. 8). These sparsely coupled oscillators can still be used successfully in the model (Fig. 11). Thus, the observation of irregular firing *in vivo* does not rule out highly regular population-level oscillations.

However, coupling is not helpful against all forms of noise. With as little as 5% of the total noise variance correlated across all cells in one network, preliminary simulations of a network of 250 gap-junction-coupled simple model neurons showed an order of magnitude decrease in stability times compared with the case of each cell receiving the same total level of fully independent noise. In fact, a network of 2500 gap-junction-coupled simple model neurons with 5% of the total noise variance correlated (and 95% independent) has approximately the same variability as a single cell with a noise variance at only 5% of σ^{2}, the value matched to the experimental data, clearly demonstrating the inability of coupling to correct this type of noise. We used independent noise sources in our simulations to avoid this problem, although realistically there is likely to be some unknown degree of correlation. This problem is an important point for future work, as it could potentially be a fundamental flaw with this approach.

Overall, our simulations suggest that a variety of realistic mechanisms can implement an oscillatory interference model, despite large amounts of independent noise in individual oscillators.

## Discussion

We demonstrated that neural oscillators with realistic levels of noise can be coupled to produce network oscillations with less variability than the individual, uncoupled oscillators and that the individual oscillators in a network may or may not show this reduction in variance. We then showed that the frequency of such an oscillatory network can be controlled to an acceptable but not perfect level of accuracy. Finally, multiple such oscillatory networks can be combined to produce grid cells with spatial firing that is relatively stable on behavioral timescales. In particular, the oscillator networks may be either gap-junction or synaptically coupled and both excitatory and inhibitory oscillator neurons can be used. We also showed that stable two-dimensional grids can be generated using biophysical model neurons, suggesting a wide range of neural models might successfully play the role of the oscillators in this model.

Previous oscillatory interference models that divided the oscillators among multiple cells (Burgess et al., 2007; Burgess, 2008; Hasselmo, 2008) assumed the postsynaptic cell *G* integrated its inputs, despite *G* often being identified with a resonant cell type: entorhinal cortex layer II stellate cells (Alonso and Klink, 1993; Giocomo et al., 2007). The distinction between resonators and integrators is a fundamental one that relates to the dynamics of the excitability of the neuron and affects the behavior of the cell in many ways (Izhikevich, 2007). Our present results suggest the oscillatory interference models might also work when the interference occurs in a resonant cell, supporting that previously untested assumption. It is interesting that the spiking of the resonate-and-fire cell in our simulations often occurred on alternate cycles (Fig. 7*C*, supplemental Fig. S5*D*; but see supplemental Fig. S7*C*, available at www.jneurosci.org as supplemental material), similarly to reports in ventrally located grid cells (Deshmukh et al., 2009).

We have made no explicit claims about the anatomical location or cell type of the oscillator networks or the grid cell proper. In principle, both parts of the model could be anywhere upstream of entorhinal cortex layer II (ECII), where the strongest grid signal to date has been found. In ECII, there is a dorsoventral gradient of grid field spacing and it is worth considering how this may relate to the model. In the present model, the intrinsic properties of the grid cell *G* do not affect the spacing of the grid fields. Instead, the spacing depends on the velocity-to-frequency transformation performed by both the velocity inputs and the oscillator networks. One possibility is that there are separate velocity signals for each VCO spatial scale. Alternatively, a single output signal might go to all VCO networks and a difference in spacing would come from changes in the *F*(*I*) slope (shallower ventral) over the range of velocity signal levels. One way to decrease the slope (increase spacing) consistent with the dorsoventral gradient along ECII (Hafting et al., 2005) is to increase the synaptic coupling in the network (Fig. 1*B*), which is consistent with increased temporal summation of EPSPs in ventral stellate cells (Garden et al., 2008). Another successful way to change the *F*(*I*) slope is to proportionally scale down all the conductance densities in the model (e.g., increase *C*_{m}), though other ways of changing the slope are known (Mehaffey et al., 2005; Fernandez and White, 2010).

Simplifications made in this work sidestepped many remaining issues in creating completely realistic oscillatory interference models. Cell-to-cell variations of membrane and synaptic conductances can be very high (Goaillard et al., 2009), and increasing oscillator heterogeneity decreases synchrony (Winfree, 1967; Cumin and Unsworth, 2007; Galán et al., 2007). Similarly, the inputs to individual VCO neurons themselves may be nonuniform, which can also interfere with synchronization (Golomb and Rinzel, 1993; Tsodyks et al., 1993; Wang and Buzsáki, 1996; White et al., 1998). Individual cells may receive inputs from different numbers of velocity-coding cells with different firing rates and with different synaptic strengths, and the result of these nonuniformities on synchrony and the frequency response of the network is difficult to predict. In contrast, homeostatic mechanisms might cause each cell to adjust its membrane properties to produce the same *F*(*I*) curve as other cells, despite dissimilar profiles of membrane conductances.

Another relevant issue is that correlated noise can synchronize oscillators (stochastic synchrony), particularly when their frequencies are similar (Teramae and Tanaka, 2004; Galán et al., 2007; Nakao et al., 2007; Marella and Ermentrout, 2008; Ly and Ermentrout, 2009). Noise in the body velocity signal might apply pressure toward synchronization of the VCO networks with each other and away from the correct phase differences (in addition to the spatial phase errors that will arise from the velocity noise itself). Further, correlated noise to neurons within networks *V*_{i} cannot be corrected through coupling, so even a small amount may have a noticeable effect on the stability of the system.

An as yet unanswered question about the grid network in animals is how many distinct spatial scales are represented (but see Barry et al., 2007). In single-cell oscillatory interference models, each grid cell can have a unique spatial scale. Continuous attractor models, however, require networks of grid cells that have identical spatial scales. With the present model, arguments of parsimony (using as few cells as possible) suggest that because large networks are needed to provide stable oscillations, there will likely be fewer VCOs than grid cells. Many grid cells would then share the same VCOs and so would be of the same spatial scale.

Our simulations provide a number of novel predictions that might be tested experimentally. First, we predict that an oscillator's input should come from cells with activity that is the inverse of the VCO's *F*(*I*) curve, with a velocity relation like that shown in Figure 1*C*. Second, when an animal makes a sudden reversal of direction (e.g., on a linear track), the firing rate adaptation may be detectable as slight changes in grid position in proportion to the animal's speed and angle of motion compared with the grid cell orientation. Also, if VCOs are shared among grid cells, drifts in the spatial firing of individual cells may be partially correlated (in proportion to the number of VCOs shared and along the directions of the shared inputs). This is in contrast to a continuous attractor network where drift among all neighboring cells would be highly correlated (Burak and Fiete, 2009) or other oscillatory interference models where the drift of each cell is independent (Zilli et al.,2009). Finally, this model may be directly tested through the use of single-cell intracellular *in vivo* recording, which is increasingly feasible (Harvey et al., 2009; Lee et al., 2009). Our model predicts that synaptic potentials recorded in the cell generating the grid pattern should resemble the postsynaptic traces we have shown in Figure 7 and supplemental Figures S5 and S7. Specifically, each grid cell should always be receiving at least three strong rhythmic inputs that change their phase relative to each other according to the animal's spatial position. Of all the model predictions, this is the most important as it offers a direct, falsifiable test of this model.

We have provided a successful biophysical simulation of an oscillatory interference grid cell model and demonstrated its robustness against realistic levels of noise. This provides a proof of concept that suggests biological systems can implement the oscillatory interference model despite the high levels of intrinsic neuronal noise, but only through experimental tests designed to distinguish between models will the true mechanism underlying the grid pattern be discovered.

## Footnotes

This work was supported by National Institute of Mental Health Silvio O. Conte Center Grant MH71702, and R01 Grants MH61492 and MH60013; National Science Foundation Science of Learning Center Grant SBE 0354378 (CELEST); and the Office of Naval Research. We thank Lisa Giocomo and Jim Heys for helpful comments.

- Correspondence should be addressed to Eric A. Zilli, Center for Memory and Brain, 2 Cummington Street, Boston University, Boston, MA 02215. zilli{at}bu.edu