Controlling Phase Noise in Oscillatory Interference Models of Grid Cell Firing

Oscillatory interference models account for the spatial firing properties of grid cells in terms of neuronal oscillators with frequencies modulated by the animal's movement velocity. The phase of such a “velocity-controlled oscillator” (VCO) relative to a baseline (theta-band) oscillation tracks displacement along a preferred direction. Input from multiple VCOs with appropriate preferred directions causes a grid cell's grid-like firing pattern. However, accumulating phase noise causes the firing pattern to drift and become corrupted. Here we show how multiple redundant VCOs can automatically compensate for phase noise. By entraining the baseline frequency to the mean VCO frequency, VCO phases remain consistent, ensuring a coherent grid pattern and reducing its spatial drift. We show how the spatial stability of grid firing depends on the variability in VCO phases, e.g., a phase SD of 3 ms per 125 ms cycle results in stable grids for 1 min. Finally, coupling N VCOs with similar preferred directions as a ring attractor, so that their relative phases remain constant, produces grid cells with consistently offset grids, and reduces VCO phase variability of the order square root of N. The results suggest a viable functional organization of the grid cell network, and highlight the benefit of integrating displacement along multiple redundant directions for the purpose of path integration.


Introduction
Grid cells recorded in the medial entorhinal cortex (MEC) of freely moving rodents fire whenever the animal occupies any one of an array of locations arranged in a regular triangular grid across the environment (Hafting et al., 2005). These cells are thought to support "path integration," providing a representation of location updated on the basis of movement-related information (O'Keefe and Fuhs and Touretzky, 2006;McNaughton et al., 2006). The oscillatory interference model  explains grid cell firing patterns in terms of input from neuronal oscillators whose frequency varies as a function of the animal's movement velocity ("velocity-controlled oscillators" or VCOs; Blair et al., 2008;Hasselmo, 2008).
Predictions for the theta-band modulation of grid cell firing rates as a function of running speed and grid scale  have been verified (Jeewajee et al., 2008), stellate cells in MEC exhibit membrane potential resonance at frequencies that vary consistently with changes in grid spacing along the dorsoventral axis of MEC (Giocomo et al., 2007), and septohippocampal "theta cells" appear to function as the predicted VCOs (Welday et al., 2011). However, a key question concerns whether the level of noise in neuronal oscillations is too high to allow the spatial stability observed in grid fields (Zilli et al., 2009;Remme et al., 2010). The observation of theta-phase precession in grid cells (Hafting et al., 2008;Reifenstein et al., 2012) implies that grid cell populations do maintain reliable spatially demarcated phase relationships in vivo. These consistent phase relationships highlight the fact that the interfering oscillations are not independent. For example, the baseline (local field potential) oscillation against which individual neuronal firing precesses may be the mean of the population activity (Burgess et al., 1993;Geisler et al., 2010), and oscillatory activity may be stabilized relative to the environment by location-specific phase-resetting inputs Monaco et al., 2011;Sreenivasan and Fiete, 2011;. Zilli and Hasselmo (2010) showed how synaptically or electrically coupling VCOs with the same preferred direction reduces their overall frequency variability. Here we show how coupling multiple VCOs carrying information about displacement along different directions (via entrainment of the baseline oscillation), or with different spatial offsets (via ring-attractor dynamics; Blair et al., 2008), could produce the phase reliability and consistency needed to generate stable and coherent grid-like firing patterns.

Materials and Methods
We use the neuronal-VCO grid cell model described previously  see Results for an overview). VCOs are simulated as neurons with a sinusoidal membrane potential oscillation (MPO) with velocitydependent frequency described by Equation 15. Three VCOs are used, with preferred directions spaced by 120°. This configuration satisfies the requirement of 60°multiple spacing for producing hexagonal grids, and enables a straightforward mechanism for maintaining VCO consistency (see Results). Each VCO spikes periodically when its phase is zero if the animal is traveling within 90°of its preferred direction (Fig. 1A). The grid cell's VCO inputs are modeled as exponentially decaying EPSPs with a time constant of 25 ms (Fig. 1C;Garden et al., 2008).
The zero-speed VCO frequency f 0 was 8 Hz, with velocity-dependent modulation scaled by ␤ ϭ 0.026 cm Ϫ1 (which is in the range of typical theta frequency modulation reported in Jeewajee et al., 2008). Phase noise is Normally distributed (SD 0.006 rad) and added to each VCO phase i every 1 ms time step, producing a minimum phase SD of 0.067 rad (i.e., 1.34 ms per 125 ms cycle of an 8 Hz oscillator). The baseline oscillation frequency f b increases from f 0 proportional to running speed (Eq. 18) and equals the mean frequency of the VCOs (see Results).

Location coding by multiple VCOs
Following Equation 16, the ith VCO's phase i relative to the baseline phase b encodes distance traveled (d i ) along its preferred direction i . If distance is represented independently along three or more directions (i.e., with three or more VCOs), the 2D location represented can be inconsistent (i.e., each pair can represent a different location; see Results; Fig. 2). Therefore, when analyzing the spatial drift of a grid due to phase noise, we assume an ongoing mechanism for realigning all phases to indicate the mean location represented by each pair. With a frequent (or continuous) realignment, such that accumulated phase noise between realignments is small (compared with their period), this is a linear transformation. Thus the mapping from phase noise offsets to mean location offset is continuous, i.e., each small phase perturbation produces a small perturbation in grid cell firing locations after realignment. The linear transformation for this realignment is described here, and an example mechanism implementing this is detailed below and in Results.
We define the state of the animal at time t as a vector v គ͑t͒ comprising its location coordinates v គ xy ͑t͒, and the phase of the baseline oscillation b ͑t͒. Then (from Eq. 16) the phases of all the oscillations obey the following: where the first n components of ͑t͒ are the VCO phases and the last component is the baseline phase, e គ͑t͒ is the phase noise, and A is the (n ϩ 1) by 3 transformation matrix: A ϭ ΄ 2␤ cos 1 2␤ sin 1 1 Ӈ 2␤ cos n 2␤ sin n 1 0 0 1 ΅ . (2) Considering the VCO phases relative to baseline (i.e., with the baseline phase subtracted), A can be viewed as projecting the location of the animal onto a plane P in the n dimensional space of VCO phases relative to baseline. The points on P correspond exactly to the relative phases in which the location represented by each VCO pair is consistent, and consequently will be inconsistent if phase drift (due to noise) takes the relative phases out of P. The closest consistent point in phase space can be found by projecting onto P along its normal vector. The location corresponding to this point v͑t͒ can be found by applying a linear transformation B, such that: which is the Moore-Penrose pseudo-inverse of A. If the phase noise in each oscillation has zero mean, applying B to noisy phases is equivalent to finding the location estimate that minimizes the mean squared distance from the correct (noise-free) location. The consistent set of phases Ј͑t͒, corresponding to this location, can be found by applying AB, i.e., Ј͑t͒ ϭ AB͑t͒.
Since VCOs encode distance relative to the baseline oscillation, common noise across all oscillations can be ignored. We assume that each oscillation has independent, identically distributed Gaussian noise. Therefore, the covariance matrix of is given by ⌺ ϭ 2 I, where 2 is the phase variance of each oscillation and I is the identity matrix of order n ϩ 1. Consequently, v is Normally distributed with covariance matrix: where the first 2-by-2 element submatrix in ⌺ v គ is the covariance matrix of v xy . We define the point of unrecoverable error as when 50% of the distribution of v xy (which is contained within Mahalanobis distance r Ϸ 1.1774, see the cumulative 2 ͑2͒ distribution of interval r 2 ) has drifted into adjacent fields (Fig. 4A).

Coupling phase-offset VCOs as ring attractors
To examine noise reduction from coupling phase-offset VCOs, we simulated an oscillating ring-attractor network of N cells (following Blair et al., 2008;for review, see Zhang, 1996). We used a 1D adaption of the 2D continuous attractor model (Burak and Fiete, 2009) with asymmetric recurrent connections making the activity bump rotate around the ring, and rate-coded neurons obeying: where i is neuron i's synaptic activation, f͑x͒ ϭ x for x Ͼ 0, 0 otherwise, and t is the noise added at time t. Neurons have positions x i evenly spaced around the ring (x i , increasing clockwise). There are two classes of neurons: clockwise or anticlockwise, in alternating positions around the ring. The recurrent connection to cell i from cell j has weight: where [ ] N indicates modulo N. W 0 is a function of the circular distance between the cells ␦, shifted by k j l and . For clockwise cells k i ϭ 1, for anticlockwise k i ϭ Ϫ1, so that the distribution of weights is shifted clockwise or anticlockwise, biases the otherwise symmetric distribution of connections so that the bump of activity propagates around the ring at a fixed frequency (see below). The weighting function is a difference of Gaussians: Recurrent inputs are inhibitory, creating a stable bump of activity due to local disinhibition and distal inhibition (i.e., ␤ Ͻ ). The feedforward input to neuron i is as follows: providing differential excitation to clockwise versus anticlockwise cells (due to k i ), which modulates the speed of propagation of the activity bump according to running velocity (m͑t͒ ϭ ␤s͑t͒cos͑͑t͒ Ϫ d ͒; see Eq. 15). Parameter values were as follows: a ϭ 25, ␤ ϭ 1/͑0.88N͒ 2 , ϭ 1.05␤, l ϭ 0.075N, ϭ 0.075N, ϭ 10 ms. t was drawn from the Normal distribution with zero mean, and SD , at each 0.5 ms time step. Values of the membrane noise level used were 0.025, 0.05, and 0.1. The activity bump had oscillation frequency ϳ8 Hz when running speed was zero.

Entraining a baseline oscillator to the mean VCO frequency
To demonstrate that activity in a baseline ring-attractor network could oscillate at a frequency entrained to the mean frequency of three VCO ring attractors, we developed a spiking network whose activity bump can be entrained to rotate at a frequency defined by its inputs. We use the Neural Engineering Framework (NEF; Eliasmith and Anderson, 2003) to derive feedback connection weights facilitating attractor dynamics for a relatively narrow Gaussian bump of activity that naturally oscillates around the ring. The NEF defines neural representations in terms of nonlinear encoding of the variable to be represented (in this case phase) by populations of neurons with tuning curve responses, and weighted linear decoding of their activities to retrieve the represented value.
We assume that neurons in the baseline ring attractor have Gaussian tuning curves in the phase value x represented by their location on the ring, such that the population is able to represent a rotating Gaussian of the form f͑x͒ ϭ exp(Ϫ (x Ϫ ͑t͒) 2 /2 2 ) with width , where ͑t͒ is the phase of the ring oscillator at time t (i.e., the phase indicated by the position of the bump of activity on the ring). Consider estimating the Gaussian activity profile in Fourier space by the vector ␥ of Fourier coefficients (where ␥ 1,n and ␥ 2,n are the nth frequency cosine and sine coefficients, respectively), using the components F0 -F6 (with small enough to make F3 significant, as this harmonic will be the target of entrainment by the VCO inputs). The oscillation of the bump around the ring, and the effect of external inputs, can be expressed in terms of the dynamics of these coefficients (Orchard et al., 2013). Noting that rotation of the activity bump corresponds to rotation of the Fourier coefficients, these dynamics are given by the following: where R is a 14-by-14 matrix composed of seven 2-by-2 rotation submatrices along the diagonal, each corresponding to coefficient pair n and with rotation angle 2nf 0 (where f 0 is the zero-speed VCO frequency, see above). The vector u គ͑t͒ is the frequency domain representation of the external input (i.e., its first seven sine and cosine Fourier coefficients).
Equation 9 can be implemented in a neural population by applying the NEF (Eliasmith and Anderson, 2003;Conklin and Eliasmith, 2005 for derivation of a similar network). Using linear integrate and fire neurons whose subthreshold voltage dynamics are given by the following: where i RC is the membrane time constant, ␦͑t͒ is the delta function, and t in is the time of the nth spike of neuron i. The cell spikes when V i Ͼ 1 and subsequently its voltage is held at zero for a refractory period, ref .
Each cell i has scalar parameters ␣ i and J i bias , the input gain and bias, respectively. The convolution ‫)ء(‬ of each input spike with an exponential decay filter, h(t), with time constant , captures postsynaptic current dynamics. The feedback weights from neuron j to neuron i required to implement Equation 9 are given by the following: where ⌬t is the simulation time step. The encoding vector for neuron i, i , specifies its tuning for coefficients in the frequency domain, i.e., the first 14 Fourier coefficients approximating a Gaussian that defines its tuning curve. The vectors i are optimal linear decoders, such that the population estimate of ␥ is as follows: These decoders are found by minimizing the mean square error of the estimate, E ϭ 1 2 ͳ ͩ ␥ Ϫ ␥ ͪ 2 ʹ with ͗ ⅐ ͘ denoting integration over the desired range of ␥ (in this case, the Fourier coefficients of Gaussians tiling phase space x).
In our simulations of this network, 200 neurons were used with ϭ 0.45rad, i RC ϭ 20 ms, ref ϭ 2 ms, ϭ 5 ms, and ⌬t ϭ 0.5 ms. Parameters ␣ i and J i bias were chosen to give each cell a random peak firing rate between 100 and 200 Hz and a random tuning width (both uniformly distributed), to reflect heterogeneity in real neuronal populations (for review, see Eliasmith and Anderson, 2003). The recurrent connectivity of the network is determined as described above. After an initial input to start it spiking, this network settles down to an approximately Gaussian bump of activity that oscillates around the ring at f 0 .
The activity profile of the external inputs to the ring is estimated by its Fourier components u គ͑t͒. Each frequency component of the external input can be used to influence the corresponding component of the activity profile of the baseline network. Thus, to entrain the F3 component of the network, the two coefficients, u 1,3 and u 2,3 , should be the cosine and sine amplitudes of the desired control oscillation, with all other components of u set to zero. The desired control oscillation should have a frequency that is three times the mean frequency of the three VCOs, i.e., its phase should be the sum of the phases of the three VCO rings (see Results). The cosine and sine Fourier components of this oscillation are simply derived from consideration of the real and imaginary components of the complex representation of the phases of each VCO ring (the product of these representations represents the desired oscillation: summing their phases), giving: where v 1 i and v 2 i are, respectively, the cosine and sine amplitudes of the ith VCO phase. The VCO amplitudes v 1 i and v 2 i could readily be decoded from VCO ring networks as inputs to the baseline network with appropriate synaptic weights, similarly as outlined above. Computation of the products of these inputs as in Equation 14 could be provided by multiplicative synapses onto neurons in the entrained network (e.g., the Sigma-pi units of Mel, 1993), or by directly (linearly) decoding the transformation in Equation 14 from a single attractor network representing the phases of all three VCO directions (e.g., by applying the NEF).

Results
The spatial firing pattern of grid cells is produced as an interference pattern between a set of VCOs with different directional tuning ( Fig. 1; see Materials and Methods; . Each VCO has a preferred direction d and an "active" frequency f a . This frequency varies according to the animal's running speed s(t) and direction (t), relative to a common "baseline" frequency f b (t): The displacement of the animal over a period of time is simply the integral of its velocity during that period, while the phase of an oscillation is the integral of its frequency. Accordingly, the phase of an oscillation whose frequency varies according to velocity along some direction will be proportional to displacement along that direction. Specifically following Equation 15, VCO phase a relative to baseline phase b encodes the animal's displacement along the VCO's preferred direction: where h(t) and (t) are the net displacement distance and angle from the origin, respectively (Fig. 1A, VCOs). Position can be tracked in 2D using two or more VCOs with differing preferred directions (Figs. 1, 2). Furthermore, the interference envelope between such oscillators will generate a periodic pattern across 2D space. With a particular spacing between their preferred directions (60°or noncollinear multiples thereof), this pattern has the same tessellated hexagonal structure as grid cell firing fields ; Fig. 1). The development of this structure, from periodic representations of displacement along specific directions, may arise from Hebbian-like unsupervised learning, which favors coincident activity from representations separated by 60° Mhatre et al., 2012), and which may also explain the presence of less regular patterns formed from superposition of these representations (Krupic et al., 2012). VCOs are neuronal oscillations with frequencies modulated by velocity-dependent synaptic inputs. They were proposed as MPOs occurring in individual neurons that synapse onto the grid cell, or within dendritic subunits of the grid cell itself ; see ), Hasselmo (2008, and Welday et al. (2011) for the former implementation, and Remme et al. (2010) for limitations of the latter. Here, we assume that the VCOs and baseline oscillations are external to the grid cell, providing synaptic inputs. The VCOs spike periodically at a particular phase of their oscillation, and might correspond to the velocity-controlled "theta cells" found along the septohippocampal circuit (Welday et al., 2011). Unlike the VCOs, the baseline oscillation does not depend on running direction and is consistent with the "somatic input" described previously . See  for discussion of the relationship of the baseline oscillation to the local field potential.
A grid cell receives inputs from multiple VCOs with preferred directions spaced at multiples of 60°and fires when input spikes arrive at similar phases ( Fig. 1; Hasselmo, 2008; i.e., the grid cell performs coincidence detection). Thus a grid cell fires at a grid node because the VCO phases align there, and will fire again when the VCO phases have all changed by a similar amount (modulo 2), i.e., when the animal has moved similar distances along the VCOs' respective preferred directions (Fig.  2). This arrangement effectively modulates the grid cell's firing rate according to the interference pattern between the VCOs, thus generating their characteristic spatial firing pattern.
We follow  in which (1) the velocity-dependent inputs to VCOs can only increase their frequency from some minimum, f 0 ; (2) the baseline oscillation takes the mean frequency of the local VCOs; and (3) a VCO only contributes to grid cell firing when the animal is running within 90°of its preferred direction. The first two constraints imply: (without affecting Eq. 15), while the final constraint allows grid cells to show thetaphase precession relative to the baseline (which is identified with the local field potential), as observed in MEC layer II (Hafting et al., 2008).

Phase noise and spatial coherence of grid cell firing
Since VCO phases encode the animal's distance from grid nodes along their preferred directions, spatially coherent grid cell firing requires that VCO phases continue to coincide at the nodes of the grid. In the presence of noise, this requires an ongoing phase-correction mechanism. One solution is for feedback from place cells, driven by location-specific environmental information, to reset the VCO phases for a particular grid cell at a grid node Monaco et al., 2011). Here we consider local mechanisms for maintaining the phase consistency of VCOs.
If VCO phases are aligned with the baseline oscillation at the origin, then the relative phase of the ith VCO encodes displacement along its preferred direction i with distance d i ͑t͒ given by the following (from Eq. 16): A B C Figure 1. The oscillatory interference model of grid cell firing. A, VCOs (red, green, and blue circles) spike periodically such that their spiking phases relative to a baseline oscillation (colored sinusoid) reflect translation along their preferred directions (black arrows) as the animal navigates its environment. Spike phase maps for a simulated run in a circular environment are shown above their respective cells (spike locations indicated by dots, color coded by baseline phase; see phase color bar). A grid cell (purple circle) is driven by VCOs with different preferred directions, and acts as a coincidence detector for spikes arriving from those inputs, by "leaky" integration (C; see Materials and Methods). The grid cell is additionally subject to modulation by the baseline oscillation (black sinusoidal arrow). These combined inputs cause the grid cell to fire when the VCOs all spike during similar positive baseline phases (which occurs within the ringed regions on each spike phase map), giving the grid cell its spatial firing fields (grid cell, spike phase map is shown below the cell, and corresponding firing rate map shown in B). Multiples of 60°spacing (here 120°) between VCO preferred directions causes the firing fields to be arranged in a regular triangular grid. Environment diameter is 1 m, grid scale is 0.6 m. The animal's trajectory is shown as a gray line. For details see Results and Burgess (2008).
Two (noncollinear) VCOs encode an unambiguous 2D location at the intersection of orthogonal lines from the endpoints of their displacement vectors ( Fig. 2A). However, the orthogonals corresponding to any additional VCOs must intersect the same point to provide a consistent representation (Fig. 2B). Any errors in their phases will cause each pair to intersect at different points (Fig. 2C). Such inconsistencies will change their interference pattern, producing smaller sub-peaks in firing around each (errorfree) grid node of the downstream grid cell. Therefore, variable inconsistency will lead to unstable and unpredictable firing fields (Fig. 3, top row; Burgess et al. , 2007Giocomo and Zilli, 2007). If three VCOs with 120°preferred direction spacing together represent a consistent location, their encoded distances will sum to zero, i.e., d 1 ϩ d 2 ϩ d 3 ϭ 0, and (Eq. 16) their relative phases will too, modulo 2. This can be seen by applying any three 120°s paced preferred directions to Equation 20, e.g., the preferred directions 1 ϭ 0Њ, 2 ϭ 120Њ and 3 ϭ 240Њ shown in Figure  2B. An error e in one of the VCO phases can abolish this zero sum property, causing the three orthogonals to intersect at the vertices of an equilateral triangle with altitude e/2␤ (Fig. 2C). If the source of error is unknown, the triangle's center best estimates location (Fig. 2D). Subtracting e/3 from all VCO phases will make their orthogonals coincide at this point, restoring consistency (and their zero sum). The situation for independent errors in all three VCOs is similar: the orthogonals intersect at the vertices of an equilateral triangle whose center best estimates location, and subtracting one-third of the total phase error from each VCO will realign them there.
Since VCO phases encode distances relative to baseline, the above error correction (specifically, subtracting one-third of the total phase error from each VCO) will happen automatically if the baseline oscillation is coupled to the VCOs to include their mean phase noise, i.e., if the baseline frequency takes the mean of the VCO frequencies (Fig. 3, bottom row). This arrangement was originally proposed in , to enable the model to function with only positive VCO frequency modulation, rather than for any contribution to noise reduction. Following this, the baseline phase at time t is as follows: where e i (t) is the ith VCO's phase error. This maintains the zero sum of relative phases (and therefore their consistency) despite the presence of noise: Tracking the mean VCO frequency could be implemented by entrainment of the baseline oscillation by the VCOs. A network mechanism to implement this is described below.
The ability to bring VCO phases into consistency by a common phase shift (e.g., by altering the baseline phase) is unique to the 3 ϫ 120°configuration. This can be seen by inspecting Figure  2D: any additional VCOs (in a multiple of 60°configuration) would be oriented in the opposite direction to one of the 120°t riples. It would thus require a correctional phase shift of the same magnitude but opposite sign, leading to a nonzero sum of phases. Furthermore, arrangements of this configuration can have a zero sum of phases with an inconsistent representation. For example, adding equal but opposite offsets to the phases of oppositely oriented VCOs makes them inconsistent with the remaining VCOs, but leaves the sum of phases unchanged. Therefore, no common baseline shift can by itself make six 60°spaced VCOs consistent. However, the principle can be extended to this configuration by combining the baseline shift (as above) with a subtraction of the mean phase between each opposing VCO from their respective phases.
For any symmetrical configuration of VCO directions, the single best estimate to use for maintaining coherent grid firing is Figure 2. Phase noise and spatial coherence. A, Each VCO represents displacement (red lines) along its preferred direction (arrows) by its phase i relative to a baseline oscillation. Two VCOs can together represent a position in 2D space, as visualized by the intersection (rat location) of perpendicular lines from the ends of the displacement vectors (red dashed lines). B, When additional VCOs are used, all perpendiculars should intersect. Note the negative displacement along the 240°preferred direction. C, If a single VCO is perturbed by phase noise (e; blue arrow), the pairwise intersections describe three separate points (blurred rat locations), which form the vertices of an equilateral triangle. The size of e is the altitude of the triangle (but would be undetectable with only two VCOs). D, If the source of the error is unknown, the best estimate of the correct location is the center of the triangle. The VCO phases can be aligned so that all three consistently represent this point by subtracting e/3 from each.  Figure 3. Effect of phase noise in VCOs on grid cell firing. Firing rate maps are shown (as in Fig. 1) of a simulated grid cell with phase noise for different trials with the same trajectory. Each trial was run in two configurations: with the background oscillation set at the noiseless mean VCO frequency (f 0 ϩ ␤s͑t͒, first row), and with the actual mean VCO frequency, including their noise (second row). The latter automatically prevents phase inconsistency forming so that grid fields remain well defined, despite the drifting of the field centers. Thus preserving representational consistency maintains spatially coherent grid firing and mitigates error accumulation. Model parameters as in Figure 1. VCO preferred directions shown by black arrows. Simulations had cumulative VCO phase noise with an ISI SD of 1.34 ms per 125 ms cycle (see Materials and Methods).
the average location (i.e., the mean vector) indicated by all unique pairwise combinations of VCOs. Realigning every VCO to indicate this location is equivalent to setting their position in relative phase space to the nearest point on the plane of consistent phases. We discuss this further in Materials and Methods, where we derive the linear transformation for computing this estimate from the phases of an arbitrary configuration of VCOs. Note that with any configuration, noise common to all VCOs and the baseline will not affect the represented location.

Phase noise and locational drift of grid cell firing
A mechanism for maintaining VCO phase consistency (e.g., following the principle described above) will ensure coherent grid cell firing. However, accumulating phase noise in any number of VCOs will still cause grid fields to drift (Fig. 3). We now characterize the drift rate for different configurations of VCOs: two VCOs with 60°spaced preferred directions, three 120°VCOs, and six 60°VCOs. In the latter two cases assuming that VCO phases are continuously aligned to the vector mean of the locations in-dicated by all pairwise combinations of VCOs to maintain coherent grid firing, as discussed above. The distribution of represented locations can be calculated from the distributions of oscillation phases (see Materials and Methods). Two main factors determine its spread: the number of VCOs and the configuration of preferred directions (Fig. 4A). Increasing numbers of VCOs result in more independent estimates of location to average over, tightening the distribution of estimated location. Increasing from two VCOs with 60°s paced preferred directions to three 120°V COs decreases the area containing 50% of the distribution by ϳ66%, and increasing from two to six 60°VCOs decreases it by ϳ83%.
We also calculated the relationship between VCO phase variance and the duration of spatially stable firing ( Fig. 4B; see Materials and Methods for derivation). Using three 120°VCOs extends the stable duration ϳ200% compared with two 60°V COs, and increasing from three to six extends this by a further ϳ100%.

Coupling VCOs with offset phases
Finally, we considered the noise-reducing effects of coupling the VCOs required to drive a population of grid cells with different spatial phases. VCOs with the same preferred directions can be configured as a ring-attractor network (Blair et al., 2008), in which a localized bump of activity oscillates around the ring at theta frequency. The bump should oscillate faster when driven by velocity-dependent inputs according to Equation 15. Such a network enforces the phase offsets between VCOs and, given a common baseline oscillation as described above, would maintain consistent spatial offsets between grids. The maintenance of fixed phase differences between each VCO in a ring means that the activity bump of each ring effectively encodes a single oscillation (with a single phase error).
We simulated such a ring-attractor network of N cells, with noisy synaptic inputs, to examine phase stability. Networks were simulated with three different levels of noise, and with N ranging from 20 to 200 (see Materials and Methods). We found that the phase variability of the network decreased with increasing N, with SD approximately of the order 1ͱN (Fig. 4C). Increasing neuronal noise in the network (i.e., increasing the membrane noise SD parameter, ) caused upward shifts in the relationship between N and phase SD.
A population of grid cells can be driven by two or more ringattractor networks containing phase-offset VCOs: one ring for each preferred direction (Blair et al., 2008). Grid cells with gridlike firing patterns of different spatial offsets can be generated by combining inputs from VCOs in each ring with different phase offsets. Furthermore, the baseline oscillation modulating a grid igure 4. Effects of VCO configuration and phase noise. A, Location estimate distribution density plots with 2 ϫ 60°VCOs (left), and 3 ϫ 120°VCOs (right). Distribution centers are aligned to a grid vertex (the correct location), with adjacent grid vertices indicated by red crosses. Blue ellipses delineate contours of equal density containing 50% of the distribution (Mahalanobis distance r Ϸ 1.1774). Both plots show the estimated location distribution once phase noise in each VCO has accumulated to reach an overall phase SD of 2 rad (ϳ39 ms per 125 ms cycle) compared with the correct phase. With two VCOs, the area enclosing 50% of the estimated location distribution is equal to the area of a hexagon at half-grid scale (dashed black lines), which we define as the point of unrecoverable error (see Materials and Methods). However, with three VCOs, the 50% area is only ϳ33% of the two VCO case. B, Stability time limits for grid representations using different configurations of VCO preferred directions as a function of oscillator phase noise (grid scale G is 0.45 m). C, Phase noise as a function of number of cells N in a ring-attractor model of VCOs with different values for the membrane noise parameter, (see Materials and Methods). SD is of the order 1 ͱN (line fits show a ͱN Ϫ b ϩ c). Data points (colored circles) show phase SD calculated from 1500 simulations of VCO phase after 5 s runs with zero running speed, i.e., after oscillating at ϳ8 Hz. Error bars show the SD of these phase SDs over five blocks of 300 trials each.
cell's firing can also be implemented as a ring-attractor network. This ring could be entrained by the VCO rings to follow their mean frequency. We developed a spiking network model to demonstrate this entrainment using the NEF (Eliasmith and Anderson, 2003), see Materials and Methods and below.

Entrainment of oscillator networks
To demonstrate the feasibility of the error correction scheme outlined above, we demonstrate that the frequency of a baseline ring attractor can follow fluctuations in the mean frequency of three VCO ring attractors, which might be caused by phase noise in the VCOs. This requires using the mean of the three phases of each VCO ring as an ongoing target to entrain the phase of the activity bump in the baseline ring. However, there are three possible solutions for the mean of the three phases, spaced 120°apart (the phase of each VCO ring is only determined modulo 360°, so the mean phase is only determined modulo 120°). Since modulation of grid cell firing by the baseline oscillation restricts firing to positive baseline phases, each solution corresponds to three different (shifted) grid patterns   Fig. 7a). To ensure stability, the baseline phase should track a single solution and not switch between different solutions. This will be achieved if the three VCO rings start with coherent phases (e.g., following reset by place cell feedback), and baseline entrainment occurs faster than the accumulation of phase error. This entrainment can be implemented by the appropriate connectivity between the VCO rings and the baseline ring attractor, as described below (see Materials and Methods for details).
We specifically targeted entrainment of the third harmonic of the baseline oscillation. As long as the oscillation of the activity bump around the ring remains coherent (i.e., the phases of the other significant harmonics are locked with the third), this achieves the desired result. Moreover, since this component is three times the base frequency, it allows us to use the sum of the phases of the three VCO rings for entraining the baseline ring, which is easier to implement than using the mean of the phases of the three VCO rings. To do this, we constructed an oscillating ring-attractor network with dynamics able to stabilize a relatively narrow Gaussian bump of activity, by calculating the appropriate feedback connectivity. Since narrower Gaussians have larger components at high frequencies, a narrow bump of activity is easier to entrain via its third harmonic. We show that, with the appropriate connection weights, the baseline oscillatory frequency can be entrained by an external input representing the sum of the phases of the three VCO rings, and indicate how this sum of phases could be derived from the VCO ring attractors (see Materials and Methods). We applied the NEF (Eliasmith and Anderson, 2003) to derive the weights necessary for feedback connections to stabilize an appropriately sized Gaussian bump in the baseline ring attractor, and for external inputs to modulate its third harmonic (Fig. 5A, schematic;Materials and Methods;Conklin and Eliasmith, 2005).
We ran simulations with the network modulated by inputs oscillating with the sum of the phases of three VCO rings (Fig. 5 B,C). For this demonstration, we used VCO ring frequencies starting at 5, 6, and 7 Hz and all increasing by 3.5 Hz over 1 s. This network was very reliably entrained to the mean phase of the VCO rings nearest to its activity bump (given the 120°ambiguity in mean phase), and thus tightly followed the mean frequency of the VCO rings.

Discussion
According to the oscillatory interference model, grid cell firing reflects interference between multiple VCOs whose phases, relative to a baseline oscillation, represent the animal's distance traveled along the VCOs' "preferred directions." Consequently, the model is sensitive to accumulating phase noise in the VCOs. When using three or more VCOs with different preferred directions, phase noise additionally corrupts grid cell firing due to inconsistency in the locations represented by pairs of VCOs (Fig. 3).
Phase noise must be controlled sufficiently for the spatial representation to be stable over short distances and durations between resets by environmental information (Etienne et al., 1996). We addressed the control of VCO phase noise in two ways: (1) coupling the baseline oscillation to VCOs with three different preferred directions to eliminate phase inconsistency between The phases of each VCO ring-attractor oscillation are represented by their sine and cosine components (dashed and solid blue curves, respectively). Multiplication of these components yields an oscillation with a phase equal to the sum of the phases of the VCO rings (green curves; see Materials and Methods for details). Cells in the baseline ring attractor (black circles) have feedback connections that stabilize a rotating Gaussian bump of activity. They additionally receive inputs carrying the oscillation corresponding to the sum of the phases of the three VCO rings weighted so that they entrain the baseline oscillation's third harmonic (see Materials and Methods and Results). B, The baseline ring attractor is entrained to oscillate at the mean frequency of the VCO rings, specifically following their mean phase (see Results). In simulations, baseline activity (shown as the linearly decoded F1, black curve) first initializes to, and then follows, the mean VCO oscillation (blue curve) as it increases from 6 to 9.5 Hz during the 1 s simulation. The entraining input, oscillating at the sum of the phases of the VCO rings, is shown in green. C, Individual cells in the baseline ring attractor are phase locked to the entrained signal. Spike rasters are shown for cells in the baseline ring during entrainment, ordered by their preferred firing phase. them and (2) coupling VCOs with the same preferred direction to reduce the effects of phase noise in individual VCOs. We then showed how both approaches could work together.

Coupling VCOs with three different preferred directions
Under the configuration of three VCOs with 120°spacing of preferred directions, coupling a modulatory baseline oscillation  to the VCOs so that it follows their mean frequency enforces consistent phases corresponding to the mean location represented by pairs of VCOs. With three VCOs, phase noise with SD ϳ2-3 ms per 125 ms cycle allows location to be stably represented for 1-2 min (Fig. 4B).
This approach follows from the observation that, as the animal moves around, changes in the phases of VCOs with symmetrically distributed preferred directions should always sum to zero. Conversely, phase changes that sum to zero are always coherent: corresponding to a grid cell representation of a shifted location. In the face of independent identically distributed phase noise in three VCOs with 120°spacing, we can optimally mitigate the effect of noise and ensure continued phase coherence, by subtracting one-third of the summed phase change in the VCOs from each (Fig. 2). This is achieved automatically by making the phase of the baseline oscillation track the mean phase of the VCOs, so that its frequency tracks their mean frequency. The zero sum of phase changes can be understood intuitively, reflecting the symmetry of the VCO's preferred directions, which ensures that the net translation along these directions sum to zero. It is also consistent with the Fourier shift theorem, if the VCO rings are identified as the Fourier components of the 2D grid firing pattern (Orchard et al., 2013).
Coupling VCOs with the same preferred direction Oscillatory MPOs measured in MEC slices have a phase SD of ϳ15 ms per cycle (Zilli et al., 2009), which would give stable grid firing for only ϳ3 s. However, the slice preparation may well disconnect grid cells from the ascending theta system that may contain coherent VCO networks (Welday et al., 2011). Furthermore, coupling neuronal oscillators to entrain one another to the same phase , or arranging them as a ring attractor such that they fire with sequentially offset phases (Blair et al., 2008;Welday et al., 2011), would achieve much lower variance and might correspond to the extended theta network operating in vivo. We show that in the latter case, a ring-attractor network of N VCOs maintains fixed offsets between their phases and increases their phase reliability by the order of ͌N (Fig. 4C) for a range of neuronal noise levels.

Combining both approaches
A population of grid cells with spatially offset grids can be driven by ring attractors containing phase-offset VCOs (Blair et al., 2008;Welday et al., 2011;. The baseline oscillation could also be implemented in a ring attractor, entrained to track the mean frequency of the VCO rings. We showed one way this entrainment could be implemented using connections between different ring-attractor networks (Fig. 5). The baseline ring attractor will precisely track the mean frequency of the VCO rings if its third harmonic is entrained to the sum of the VCO phases. We demonstrated network connectivity between ring attractors that would be able to implement this.
We also discussed how configurations with additional VCO directions could be implemented, such as the 6 ϫ 60°configuration in which each VCO ring could be phase coupled to the VCO ring with the opposing preferred direction, without compromising the baseline entrainment mechanism. Note that Figure 4 shows increases in performance relating to both phase noise cancellation by an entrained baseline (e.g., when going from the 2 ϫ 60°to the 3 ϫ 120°c onfiguration) and from increasing the number of coupled VCOs (e.g., the ͌N increase in stability in a ring attractor when going from the 3 ϫ 120°to the 6 ϫ 60°configuration).

General implications
The oscillatory interference model can be seen as mapping location in 2D space onto the phase differences between VCO rings . From this perspective, our configuration of three VCO rings with 120°spacing represents location by the phase differences between the VCO rings and the baseline. The dimensionality of the phase difference code (3) is reduced to match that of 2D space by the constraint that the sum of phase differences is constant, ensuring a continuous mapping between phase and location. The baseline in this configuration can be thought of as a fourth VCO ring, which could also be velocity modulated similarly to other VCO rings. As long as the other three VCOs are symmetrically arranged around the "baseline" in terms of their modulation by the direction and speed of movement, then entraining the baseline to the mean phase of the other VCO rings will have the desired effect. We note that our analyses concern mitigating the effects of phase noise within a "module" of grid cells and VCOs with the same grid scale. The coherence of representations across modules of different scales is not considered here, but does have interesting implications for decoding error (for review, see Sreenivasan and Fiete, 2011;Towse et al., 2014).
Our implementation of baseline phase entrainment suggests a functional role for harmonics of theta frequency oscillations, which may, for example, manifest as theta-gamma frequency coupling (Belluscio et al. 2012;Lisman and Jensen, 2013). This may reflect an important aspect of oscillatory interference models: that the proposed oscillations are tightly coupled, rather than independent oscillations, which would rapidly decohere. Such coupling suggests a role for local circuits, and for phase coordination across the septohippocampal theta network (Blair et al., 2008;Welday et al., 2011). Failure to detect theta rhythmicity in the firing of grid cells in bats (Yartsev et al., 2011) argues against oscillatory models; however, rhythmicity may be undetectable due to the very low firing rates and movement speeds in this study (Barry et al., 2012). Additionally, the absolute frequency band of oscillations in the model is unimportant, with some suggesting that delta-band rhythmicity plays the role of rodent theta in bats (Heys et al., 2013) and humans (Watrous et al., 2013), while a baseline frequency of zero corresponds to a nonoscillatory network that integrates self-motion along multiple preferred directions (Mhatre et al., 2012).
Without our scheme to maintain phase consistency between VCOs, the amount of inconsistency between three or more VCOs (e.g., the area of the triangle of pairwise intersections; Fig. 2C) will increase as noise accumulates, potentially signaling the uncertainty in represented location. This measure of uncertainty is another advantage of tracking distance redundantly along more than two directions . Similarly, sailors obtain a position fix using multiple "lines of position" from landmarks, enabling them to gauge uncertainty: using three lines produces a triangular region (a "cocked-hat") whose size reflects measurement inconsistency (Cort, 2009). Common noise in VCO phases will not affect the spatial firing pattern, although common errors in self-motion information will cause drift in any path integration system, and can only be corrected by environmental input (Cheung et al., 2008).
In summary, the oscillatory interference model  produces coherent grid cell firing patterns that remain stable for seconds in the face of reasonable assumptions regarding phase noise in an isolated system. Extending this scheme to the in vivo situation, in which location-specific information (perhaps mediated by place cells) might reset VCO phases Monaco et al., 2011;Sreenivasan and Fiete, 2011; and VCOs with constant phase offsets could be arranged as ring attractors (Blair et al., 2008), potentially provides unlimited additional stability. This situation is more consistent with the spatially stable theta-phase of firing observed in grid cells (Hafting et al., 2008;Reifenstein et al., 2012) and largescale coherence of the septohippocampal theta system (Buzsáki, 2002) than analyses of MPOs in vitro (Zilli et al., 2009;Remme et al., 2010). The provision of a local baseline oscillation might be a role for local circuits in MEC, while the septohippocampal theta system may contain the proposed coherent sets of VCOs (Welday et al., 2011).