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 vpeak 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, xexp(t) and yexp(t), sampled at resolution dtexp = 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 dtexp, so the velocity signals were linearly interpolated, producing velocity signals vx(t) and vy(t).
To adjust the smoothness of the trajectory, vx(t) and vy(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) =
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. 1A, black line, uncoupled noisy cell).
Firing frequency versus input current for individual and networks of simple model cells. A, The firing frequency of an individual cell without noise (light gray) is very similar to the response of an individual noisy cell (black), but even more similar to a network of 250 noisy, gap-junction-coupled neurons (dark gray). For clarity, an inset with an expanded y-axis is shown. B, A network of 250 noisy, synaptically coupled neurons fires at much lower frequencies than any of the cells in A for the same input magnitude. The entire F(I) curve also moves downward as the synaptic coupling strength g is increased (g = 192, light gray; g = 256, dark gray; g = 384, black) (compare with supplemental Fig. S4A, available at www.jneurosci.org as supplemental material). C, Since the firing rate, f, is directly related to velocity, v, v = (f − ωb)/β, the inverse of the F(I) curve (shown for the synaptically coupled network) describes the input an oscillator network requires to encode different velocities along its preferred direction. General parameters: 15 s simulations and see Table 1. Single unit: n = 1, w = 1.2, σ = 0, or σ = 100. Gap-junction coupling: n = 250, w = 0.0048, σ = 100. Synaptic coupling: n = 250, w = 0.0048, σ = 100, ωb = 6.4405 Hz.
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, fi(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 Vi. The frequency of Φi was set on each time step to fi(t), as given above. In theory, if Vi is also at frequency fi(t) at time t and if Vi and Φi are at phase 0 at time 0, then Φi and Vi 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 Vi. In reality, Vi's frequency will not be perfectly controllable, so the difference in phase between Φi and Vi is a measure of the phase error Vi has accumulated. We recorded this error (the difference in phases) each time any cell in Vi 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.
Default parameters
Architecture of the model. The input to the model is a set of nVCO + 1 velocity signals, Veli. These signals go to respective networks Vi that are all-to-all internally coupled over connections Ci, which may be either synapses or gap junctions. The grid cell G receives input over synapses Ci,G from all cells in all Vi.
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 Vi, 0 ≤ i ≤ nVCO, which are identical apart from their inputs. V0 will often be called the baseline network oscillator.
Cells within each network Vi 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 Vi use the same connectivity matrix for reasons of computational efficiency). The cells in each Vi are the only cells that receive external input or noise. All cells in each Vi also project onto G with identical synapses of strength wi. When synapses from all networks have the same strength (i.e., wi = wj for all i, j), the weight is simply referred to as w for conciseness. Importantly, the networks Vi 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 V0 is assumed to have a privileged connection onto G such that inputs from any other oscillator Vi, i ≠ 0, have no effect unless they occur in a tgate second window after any spike from a cell in V0. Thus, a gate that allows inputs from the other Vi 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, cres ± ωres
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 Vi 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 Cm) 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 Vi receive the same scalar input signal IVel,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 Vi. Given a desired frequency fi, defined above, and the frequency response of the network to a given input, F(IVel,i), the input to all cells in network Vi is IVel,i(t) = F−1(fi).
The primary output measure of the model is the firing of cell G. When a single network V0 drives G (using G as an indirect measure of the behavior of V0), 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 2nσ2. The variance is converted from units of seconds2 to radians2 (2nσ2(2π)2/μ2 radians2) and the estimated stability time is the amount of time for the variance to hit a threshold of 2.5 radians2 (for further details, see Zilli et al., 2009).
When two or more networks Vi drive G, the inputs to Vi 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 Vi 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, vr, vt, vpeak, 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 Cm 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 tbase 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 cres was selected similarly, and the resonant frequency was ωres = ωb. Note that when cres ≪ 0, the cell is essentially nonresonant, so we kept cres ≈ 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 Vi 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 V0, which maintained a constant frequency and a number of active oscillators Vi, 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 fi(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 nVCO = 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. 1A) 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.
Two-dimensional grid cell with oscillators comprising single, noise-free simple model neurons. A, The network receives velocity input corresponding to a smooth two-dimensional trajectory (gray) and a leaky integrator postsynaptic cell generates spikes (black) as output. B, The autocorrelogram of the spatial firing in A shows the clear hexagonal pattern characteristic of grid cells. C, D, A simulation of the abstract model produces essentially identical results to A and B apart from spatial displacement and difference in field size. E, F, Phase differences between simple model and abstract model of the two active velocity-controlled oscillators show a linear drift over time, and the fluctuations correspond to the velocities along the preferred directions of the respective VCOs. E, Black line, For comparison, the horizontal component of velocity that is encoded by VCO 1. G, The baseline oscillators accumulate no phase error in this noiseless simulation. Parameters: 320 s simulation, nVCO = 2 (active VCO preferred directions of 0 and 2π/3 radians), n = 1, τ = 40 ms, w0 = 0.8, w1 = 0.14, ωb = 7.8989 Hz; Table 1.
In Figure 3, input from the baseline oscillator neuron V0 is made particularly strong and the inputs from the two active oscillators V1 and V2 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 3E. 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 5A. 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).
Noisy simple model neurons, when uncoupled, are unfit as oscillators. A, C, The network receives velocity input corresponding to a smooth two-dimensional trajectory (gray), and the locations where the simulated grid cell spikes over the entire 320 s (A) and first 80 s (C) are indicated with black dots. B, D, The autocorrelogram of the spatial firing shows no clear regularity, even early in the simulation. E–G, Phase differences between corresponding oscillators show that realistically noisy neurons lose their correct phases very quickly compared with behavioral timescales. Parameters: 320 s simulation, nVCO = 2 [active VCO preferred directions of 0 and 2π/3 radians], n = 1, σ = 100, τ = 40 ms, w0 = 0.8, w1 = 0.14, ωb = 7.8989 Hz; Table 1.
ISI histograms for cells in three conditions. A, ISIs of an uncoupled, noisy simple model neuron are much more variable than when multiple cells are coupled (as in B and C). B, ISIs of one noisy neuron from a network of 250 gap-junction-coupled neurons. C, ISIs of one noisy neuron from a network of 250 synaptically coupled neurons. D, ISIs of a LIF cell driven by the cell in A. The histograms are identical, showing that the postsynaptic cell directly inherits the variability of the cell driving it. E, ISIs of an LIF cell receiving input from all 250 gap-junction-coupled neurons. The variability of individual gap-junction-coupled neurons is very similar to the variability of the postsynaptic cell. F, ISIs of an LIF cell receiving input from all 250 synaptically coupled neurons. The individual cells (C) are noticeably more variable than the postsynaptic cell (and see Fig. 8 for a more extreme case). Parameters: 60 s simulations, n = 1 (w = 1.2) or n = 250 (w = 0.0048), σ = 100; Table 1.
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 Vi 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 V2 is consistent with the phase errors in V1 (Fig. 6E) but not V0 and V2 (Fig. 6F,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. 8A).
Two-dimensional grid cell with network oscillators made of noisy, synaptically coupled simple model neurons. A, The network receives velocity input corresponding to a smooth two-dimensional trajectory (gray) and a gated, leaky integrator postsynaptic cell generates spikes (black). B, The autocorrelogram of the spatial firing in A shows hexagonal periodicity. Parameters: tgate = 0.001 s, τ = 5 ms, w = 0.0016. C, D, The spatial firing and spatial autocorrelogram of the same network but with a resonant postsynaptic cell. Parameters: cres = −0.01, ωres = ωb, w0 = 0.0024, w1 = w2 = 0.002. E–G, Phase differences between network and abstract models of the two active velocity-controlled oscillators. The average of the phase differences remains relatively constant, although individual cells are firing over a much wider range of phases compared with when the cells are instead gap-junction-coupled. General parameters: 280 s simulation. nVCO = 2 [active VCO preferred directions of 0 and 2π/3 radians], n = 250, σ = 100, ωb = 7.2543 Hz; Table 1.
Combining the activity of multiple synaptically coupled oscillators. A, Voltage traces of single neurons from V0 (black), V1 (light gray), and V2 (dark gray) in the first 4 s of the trajectory shown in Figure 6. B, The neurons in A project onto an integrate-and-fire neuron where cells in V0 do not produce activity, but instead open a gate allowing cells in V1 and V2 to produce depolarization and sometimes spikes (black). Parameters: tgate = 0.001 s, τ = 5 ms, w = 0.0016. C, The neurons in A project onto a resonate-and-fire neuron with resonant frequency equal to the baseline frequency. Parameters: cres = −0.01, ωres = ωb, w0 = 0.0024, w1 = w2 = 0.002. General parameters: nVCO = 2 (active VCO preferred directions of 0 and 2π/3 radians), n = 250, σ = 100, ωb = 7.2543 Hz; Table 1.
A large network oscillator with a low connection probability can produce low-variance network oscillations while individual cells remain highly variable. A, Voltage trace over a short period of time of one of 5000 cells in a network where there is a 1% probability that any cell X connects to any other cell Y. The firing is irregular and the effects of noise are clear. B, Activity of postsynaptic integrate-and-fire cell receiving input from all 5000 oscillator cells including the one in A. The inputs to the postsynaptic cell are fairly spread out in time during each period so the cell fires in a burst. Comparing A to B, it can be seen that an individual cell will change its time of firing relative to the rest of the population on a cycle-by-cycle basis. C, A histogram of the ISIs of the postsynaptic cell from B shows that the period between bursts is highly consistent (period SD, 0.0008 s) despite the irregular inputs. This allows the oscillator to be used successfully in full spatial grid simulations (Fig. 11). D–H, ISI histograms for five cells in the network oscillators show the same mean ISI as the postsynaptic cell but much larger variability. The median period SD over all 5000 cells was 0.015 s, ∼20 times that of the postsynaptic cell in B and approximately half of the period SD of the uncoupled cells. Parameters: 10 s simulations, I = 95.8, n = 5000 cells, w = 0.0006, 1% connectivity; Table 1.
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 Vi 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 Vi neurons as gap-junction-coupled inhibitory neurons and G as a simple model neuron driven to fire at the baseline frequency. The inhibition from Vi 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.
Inhibitory noisy simple model neurons as oscillators. A, The network receives velocity input corresponding to a smooth two-dimensional trajectory (gray) and a spontaneously spiking simple model postsynaptic cell generates spikes (black). B, The autocorrelogram of the spatial firing in A demonstrates clear hexagonal periodicity. C, Voltage traces of single cells from gap-junction-coupled networks V0 (light gray), V1 (dark gray), and V2 (black). D, The cells in each Vi project with inhibitory synapses onto a simple model cell firing at ωb. When the Vi are out of phase, the cell is tonically inhibited, but when the Vi move into phase, there is enough time for the cell to fire between volleys of inhibition. Parameters: 320 s simulation, nVCO = 2 (active VCO preferred directions of 0 and 2π/3 radians), n = 250, σ = 100, ωb = 8.9735 Hz, w0 = −3.2, w1 = w2 = −0.76; Table 1.
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 INa, IK, INaP, IH fast, and IH slow). The output of this simulation is at least as good as the previous simulations.
Noisy biophysical neuron networks (synaptically coupled) produce usable velocity controlled oscillators. A, The grid cell is modeled as a gated, leaky integrate-and-fire cell and its spiking (black) is shown along a 2D trajectory (gray). B, The autocorrelogram of the spatial firing in A confirms the clear hexagonal periodicity. C–E, The phase differences for all oscillators remain fairly steady during the simulation, demonstrating that the biophysical model is capable of implementing the oscillatory interference mechanism. Parameters: 240 s simulation. nVCO = 2 (active VCO preferred directions of 0 and 2π/3 radians), n = 250, g = 80, σ = 3.44, tgate = 0.040 s, τ = 50 ms, w = 0.0024, ωb = 8.9519 Hz; Table 1.
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 5A is comparable to the experimental subthreshold oscillation period histograms in Giocomo and Hasselmo's (2008) supplemental Figure S5A (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. 5C,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.
Sparsely coupled simple model neuron networks as oscillators. A, The spatial trajectory (gray) provides the input to three oscillator networks of n = 5000 neurons with a 1% connectivity probability. Though the individual cells in the networks are highly variable (Fig. 8), the network-level variability is low so the integrate-and-fire postsynaptic cell spikes (black) in a stable grid pattern during the 240 s. Parameters: τ = 12 ms, w0 = 0.00019, w1 = w2 = 0.00018. B, Spatial autocorrelogram of the spiking in A. C, The resonate-and-fire postsynaptic cell is also able to combine its inputs to produce a stable grid. Parameters: cres = −0.01, wres = ωb, w = 0.001. D, Spatial autocorrelogram of the spiking in C. E–G, The highly variable individual oscillators spike over most phases, but the average phase of the network shows very little phase drift compared with the abstract model. Every 25th phase difference is plotted due to the large number of spikes. Parameters: 240 s simulation. nVCO = 2 [active VCO preferred directions of 0 and 2π/3 radians], n = 5000, p = 0.01, g = 256, σ = 100, ωb = 8.7942 Hz; Table 1.
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. 7C, supplemental Fig. S5D; but see supplemental Fig. S7C, 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. 1B), 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 Cm), 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 Vi 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 1C. 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
