Abstract
The rodent septohippocampal system contains “theta cells,” which burst rhythmically at 4–12 Hz, but the functional significance of this rhythm remains poorly understood (Buzsáki, 2006). Theta rhythm commonly modulates the spike trains of spatially tuned neurons such as place (O'Keefe and Dostrovsky, 1971), head direction (Tsanov et al., 2011a), grid (Hafting et al., 2005), and border cells (Savelli et al., 2008; Solstad et al., 2008). An “oscillatory interference” theory has hypothesized that some of these spatially tuned neurons may derive their positional firing from phase interference among theta oscillations with frequencies that are modulated by the speed and direction of translational movements (Burgess et al., 2005, 2007). This theory is supported by studies reporting modulation of theta frequency by movement speed (Rivas et al., 1996; Geisler et al., 2007; Jeewajee et al., 2008a), but modulation of theta frequency by movement direction has never been observed. Here we recorded theta cells from hippocampus, medial septum, and anterior thalamus of freely behaving rats. Theta cell burst frequencies varied as the cosine of the rat's movement direction, and this directional tuning was influenced by landmark cues, in agreement with predictions of the oscillatory interference theory. Computer simulations and mathematical analysis demonstrated how a postsynaptic neuron can detect locationdependent synchrony among inputs from such theta cells, and thereby mimic the spatial tuning properties of place, grid, or border cells. These results suggest that theta cells may serve a highlevel computational function by encoding a basis set of oscillatory signals that interfere with one another to synthesize spatial memory representations.
Introduction
The hippocampus and surrounding cortex contain neural circuits that store memories for facts and past experiences (Eichenbaum and Cohen, 1992; Squire and Schacter, 2002). In rodents, these regions contain neurons that fire selectively at preferred locations in space and might thus encode memories of familiar spatial environments (O'Keefe and Nadel, 1978; McNaughton et al., 2006). Several categories of spatially tuned neurons have been identified: place cells fire at one or a few preferred locations (O'Keefe and Dostrovsky, 1971), grid cells fire at multiple locations forming a hexagonal lattice (Hafting et al., 2005), and border cells fire in fixed relationships with environmental boundaries (Savelli et al., 2008; Solstad et al., 2008; Lever et al., 2009). These neurons are believed to participate in computing the animal's location by integrating its movement velocity over time, a process known as path integration (McNaughton et al., 1996; Etienne and Jeffery, 2004).
Spike trains of spatially tuned neurons are often modulated by 4–12 Hz theta oscillations, which have been proposed to play a key role in memory processing (O'Keefe and Recce, 1993; Buzsáki, 2006; Düzel et al., 2010; Hasselmo et al., 2010; Rutishauser et al., 2010; Bissiere et al., 2011). Burgess et al. (2005, 2007) introduced an “oscillatory interference” theory, hypothesizing that theta oscillations are generated by velocitycontrolled oscillators (VCOs), which perform path integration by modulating their frequencies in proportion with the speed and direction of a rat's translational movements. Supporting this idea, theta frequency is indeed modulated by a rat's movement speed (Rivas et al., 1996; Geisler et al., 2007), and oscillatory properties of spatial neurons are correlated with their spatial tuning parameters in accordance with predictions of oscillatory interference models (Burgess et al., 2007; Giocomo et al., 2007; Jeewajee et al., 2008a; Zilli et al., 2009). However, oscillatory interference models explicitly require that VCO frequencies vary as the cosine of an animal's movement direction, and such directional modulation of theta oscillations has never been observed.
Here, rhythmically bursting theta cells were recorded from medial septum, hippocampus, and anterior thalamus in behaving rats. We found that burst frequencies of theta cells were modulated by the rat's movement direction with cosinelike tuning, and directional tuning functions followed rotations of landmark cues, suggesting that theta cells might encode VCO signals predicted by the oscillatory interference theory. Computer simulations revealed that a postsynaptic neuron could exhibit spatially selective firing by detecting locationdependent synchronization among inputs from theta cells with firing properties similar to those observed in our experiments. The postsynaptic neuron could mimic the firing rate maps not only of grid cells, as in prior oscillatory interference models (Burgess et al., 2007; Giocomo et al., 2007; Hasselmo et al., 2007; Burgess, 2008; Zilli and Hasselmo, 2010), but also of place and border cells. Based on these results, we propose how a network of central pattern generator (CPG) circuits composed exclusively from theta cells could provide a basis set of VCO signals for generating diverse populations of spatially tuned neurons such as place, grid, and border cells.
Materials and Methods
All experiments were conducted in accordance with the U.S. National Institute of Health Guide for the Care and Use of Laboratory Animals (NIH Publications No. 8023), and were approved in advance by the animal subjects review committee at the University of California, Los Angeles. Technical descriptions of computer simulations and neurophysiological data analysis (including source code) are available on the ModelDB database (Hines et al., 2004) under accession number 129067 (http://senselab.med.yale.edu/modeldb/ShowModel.asp?model=129067).
Subjects and surgery.
Male Long–Evans rats weighing 350–400 g were housed singly and reduced to 85% of ad libitum weight through limited daily feeding, then trained over 5 d to forage for food pellets in an enclosed environment (see below, Recording sessions and behavior tracking). Under deep isoflurane anesthesia, rats were chronically implanted with tetrode arrays targeting medial septum and dorsal hippocampus (three rats) or anterior thalamus (two rats). Each rat was implanted with16 tetrodes (64 electrode channels), grouped into four independently drivable bundles consisting of four tetrodes each.
Recording sessions.
For data analysis purposes, a “recording session” (also referred to as a “recording” or “session” for short) is defined here as an uninterrupted period of singleunit recording that began when data acquisition was initiated and ended when the experimenter terminated data acquisition just before removing the rat from the experimental environment. Throughout each recording session, rats foraged for 20 mg purified food pellets (Bioserv) in one of three maze environments: (1) a small cylinder (80 cm diameter, 60 cm high, with black walls, a black floor, and a white cue card), (2) a large cylinder (200 cm diameter, 60 cm high, with black walls, a gray floor, and a light blue cue card), or (3) a small square (50 × 50 cm, with white walls, a white floor, and a black cue card). All three mazes were centered within a 2 × 2 m square enclosure surrounded by a black curtain, with light provided by a 40 W bulb mounted on a stand in the corner of the enclosure where the experimenter entered and exited. When a cue card was rotated to assess the influence of visual landmarks on theta cells, the light and entry position were rotated by the same angle to maintain coherence among spatial cues. When the maze environment was swapped between sessions (for example, from cylinder to square), extramaze cues (light and entry position) remained fixed.
Only sessions that met requirements for adequate behavioral sampling in all movement directions were admitted for data analysis (see below, Data inclusion criteria). To help ensure that minimal criteria for behavioral sampling in all movement directions were met, sessions were usually continued for as long as possible, until the experimenter determined by visual observation that the rat was no longer sufficiently motivated to forage for food pellets. A slow rate of food delivery (one pellet per 30–45 s, dropped at pseudorandom intervals) prevented rats from becoming sated too quickly, and thus extended the average duration of the sessions.
Video tracking analysis.
Rats wore a pair of red and green lightemitting diodes (LEDs) spaced 11.25 cm apart from one another, and an overhead video camera sampled LED positions at r = 30 Hz with a resolution of either p = 1.7 pixels per centimeter (hippocampal and medial septum recordings) or p = 4.7 pixels per centimeter (anterior thalamus recordings). Each LED's position was smoothed using a boxcar window 15 samples (0.5 s) wide before computing the midpoint between them. Another iteration of smoothing was performed on the midpoints, using a boxcar window 30 samples (1.0 s) wide. The smoothed midpoints were taken as the rat's estimated position x_{t} = (x_{t}, y_{t}) at each time sample, t, and movement velocity (in centimeters per second) at each sample was estimated by the following: The rat's running speed (in centimeters per second) was estimated by the following: The rat's movement direction was estimated by the following:
Singleunit acquisition.
Theta cells were recorded using a DigitalLynx Sseries acquisition system (Neuralynx). Singleunit waveforms were isolated by manual cluster cutting using Spikesort3D (Neuralynx) software running on a Windows PC, and were required to meet minimum criteria for unit isolation and theta rhythmicity (see below, Data inclusion criteria). Spike trains recorded during different sessions were considered to be from the same theta cell if (1) they were obtained from the same tetrode, (2) the tetrode had been advanced <80 μm between recordings, and (3) cluster boundaries and waveform shapes were visually similar on all tetrode channels for both sessions. Spike trains recorded from the same tetrode during different sessions were considered to be from different theta cells only if they were recorded at coordinates >300 μm apart. In one case, spike data from a session were excluded from analysis because the clusters and waveforms looked similar to a prior recording on the same tetrode, but the tetrode had been advanced by <300 μm since the prior session, so it was unclear whether this was an old or new cell.
Balanced running speed distributions.
To analyze modulation of theta cell bursting by movement direction, the 360° range of movement directions was subdivided into eight 45° wide bins centered at 0, 45, 90, 135, 180, 225, 270, and 315°. Tracking data from each session were parsed to extract all movement epochs containing 12 consecutive position samples (that is, 0.4 s) satisfying two conditions: (1) the rat's movement direction remained within the same bin throughout the epoch, and (2) the rat's running speed was >7.5 cm/s throughout the epoch. Hence, each unidirectional movement epoch was an episode lasting exactly 0.4 s, during which the rat was moving continuously in the same directional bin at a speed of >7.5 cm/s. The rat's mean running speed, s̄, during each unidirectional movement epoch was computed by averaging ‖v_{t}‖ in Equation 2 across all 12 samples in the epoch: where i indexes each position sample. The range of movement speeds from 7.5 to 50 cm/s was evenly subdivided into 17 bins (bin width, 2.5 cm/s). Every epoch was classified into one of these speed bins according its s̄. The total amount of time, z_{d,s}, that the rat spent running in direction bin d and speed bin s during the session was computed as follows: where 0.4 s is the duration of each movement epoch, and N_{d,s} is the total number of epochs in direction bin d and speed bin s. Since there were 17 speed bins, the speed distribution for direction bin d was a 17element vector: To assure adequate sampling of movements in all directions, recording sessions were not included in the data analysis unless the area under z_{d} exceeded 20 s for all directions d.
A uniform distribution of movement speeds in all directions was enforced by sampling movement epochs in such a way that z_{d} = Z for all d, where Z is referred to as the balanced running speed distribution: where Z_{s} is the “balanced” amount of time spent running in speed bin s, defined as where H is the Heaviside function, and ε is an arbitrarily small threshold so that H evaluates to zero whenever min(z_{1,s}, z_{2,s}, …, z_{8,s}) = 0. In Equation 8, multiplication by H[min(z_{1,s}, z_{2,s}, …, z_{8,s}) − ε] sets the value in speed bin Z_{s} to zero if speed bin s is equal to zero in any of the eight directional speed distributions, z_{d} (that is, if z_{d,s} = 0 for any d). Multiplication by max(z_{1,s}, z_{2,s}, …, z_{8,s}) assigns the remaining nonzero speed bins to take on the maximum value for speed bin s that can be found among the eight directional speed distributions, z_{d}. Using Equations 7 and 8 to derive Z (in combination with resampling methods described below, Autocorrelograms and spectral analysis), theta cell burst frequencies could be analyzed in such a way that the speed distribution was uniformly equal to Z for all movement directions. The mean balanced running speed for each session was computed as follows: where c_{s} denotes the center speed of bin s.
Autocorrelegrams and spectral analysis.
To compute directional autocorrelograms, a set of “epoch autocorrelograms” was first created (each from spikes fired during a single 0.4 s epoch of unidirectional movement) by computing a vector of 513 interspike interval counts spanning the time range from −0.4 to +0.4 s (bin width, 1.56 ms). Since each movement epoch's duration was exactly 0.4 s, the autocorrelograms' tails tapered to zero at their boundaries (−0.4 and +0.4 s), so that artificial tapering by a window function was not necessary when later taking the FFT of autocorrelograms to estimate burst frequencies (see below).
To enforce the requirement of constant running speed across all movement direction bins, a resampling method was used in conjunction with the balanced running speed distribution to extract epochs for analysis in such a way that the mean running speed was identical in all directions. To begin, all epoch autocorrelograms from the session were pooled together, and z_{d,s} (see Eq. 5) was initialized to zero for all d and s. Epoch autocorrelograms were then randomly sampled from the pool without replacement. At each sampling, if z_{d,s} < Z_{s}, then z_{d,s} was incremented by 0.4 s, and the epoch autocorrelogram was averaged into a cumulative autocorrelogram for its directional bin, d; otherwise, the epoch autocorrelogram was discarded. When the pool was depleted, it was checked whether z_{d,s} < Z_{s} for any z_{d,s}; if so, then all epoch autocorrelograms were returned to the pool for another round of random sampling without replacement (thus, some epoch autocorrelograms could be averaged into the composite autocorrelogram more than once). When z_{d,s} = Z_{s} for all z_{d,s}, calculation of the speedbalanced directional autocorrelograms was complete.
An FFT of each speedbalanced autocorrelogram was then taken to derive a power spectrum from which burst frequencies could be estimated. Each autocorrelogram was padded with zeros to a length of 2^{19} elements to yield a frequency bin width of 0.0012 Hz for the FFT [as in the study by Jeewajee et al. (2008a)]. After computing the FFT, frequency elements were multiplied by their complex conjugates and then divided by 2^{19} to obtain the power value at each frequency. After obtaining the power spectrum for each directional bin d, the composite autocorrelograms were discarded, and then they were then regenerated by a fresh iteration of resampling using the speedbalancing algorithm described above. The power spectra of the fresh autocorrelograms were then computed, and this cycle of refreshing the autocorrelograms and recomputing their power spectra was repeated for 100 iterations. Multiple iterations were necessary because, as explained above, the speedbalancing algorithm selected movement epochs at random for inclusion in the composite autocorrelogram, and this produced variability in the results of the power spectrum analysis on each iteration. Benchmarking tests confirmed that averaging over 100 iterations was sufficient to yield an accurate estimate of the true mean around which single iteration estimates were varying. For each directional bin, the 100 power spectra obtained from individual iterations were averaged together, and the final averaged power spectrum was smoothed by a boxcar window 14 frequency bins wide.
Power spectra had to meet the criterion of theta rhythmicity in all eight movement directions, or else the session was not admitted for analysis (see below, Data inclusion criteria). To estimate theta cell burst frequencies, the frequency bin of each power spectrum with the highest power in the 5–11 Hz band was located. The power spectrum was then thresholded at 50% of this peak (see Fig. 2E, black shaded regions), and the expected frequency value of the suprathreshold area of the power spectrum was obtained by the following:
where F is the estimated burst frequency; L and U are the frequency bin indices for the lower and upper frequency boundaries, respectively, of the suprathreshold region; f_{i} and p_{i} are the frequency value and suprathreshold power, respectively, of the ith power spectrum bin; and C =
DBFT curves.
To analyze how a theta cell's burst frequency was influenced by the rat's running speed in different directions, the eight points of its DBFT curve were subdivided into three movement direction categories: preferred (the three points nearest to θ̄), antipreferred (the three points 180° opposed from θ̄), and orthogonal (the two remaining points, approximating directions ±90° from θ̄). Within each direction category, epochs were subdivided into five speed categories according to the rat's mean running speed s̄ (Eq. 4) during the epoch <10, 10–15, 15–20, 20–25, and 25–30 cm/s. Movement epoch autocorrelograms were then computed independently for epochs within each speed/direction category, and burst frequencies were estimated from the autocorrelograms by the same methods used for obtaining DBFT curves (see above). Estimated burst frequencies were adjusted to a common baseline across cells by subtracting each cell's estimated burst frequency in the 10–15 cm/s speed bin for the orthogonal direction category from all speed bins in all direction categories. Burst frequencies for each speed bin within each direction category were then averaged together across sessions and then cells. Finally, the populationaveraged data for each speed bin and direction category were all shifted on the yaxis by a common factor that caused the yintercept of a linear fit to the all speed data (averaged within speed bins across directional categories) to pass through y = 0, yielding the graph in Figure 5D. Since each autocorrelogram included only a narrow range of running speeds, it was necessary to protect against estimation errors caused by undersampling. To achieve this, three criteria were enforced for inclusion of data from a particular speed bin of a particular cell in the population average: (1) to assure sufficient behavioral sampling, the session had to yield at least 10 s of data (that is, 25 movement epochs of 0.4 s each) for that speed bin; (2) to assure sufficient spike sampling, the mean firing rate of the cell had to be at least 20 Hz when averaged across all of the movement epochs included in that speed bin; (3) to assure that data for each direction category were not drawn from different subpopulations of theta cells in the population average, data for a given speed bin from a given cell had to meet criteria 1 and 2 in that speed bin for all three direction categories to be included in the analysis.
Each session's DBFT curve was fit to a cosine function, and fitted parameters were used to analyze modulation of theta rhythm by running speed for different directions (see Results). A significance level for the cosine fit was computed by taking the eight data points on the DBFT curve, resorting them into all possible permutation orders, and repeating the cosine fit for each permutation. The r^{2} value measuring the goodness of fit for each permutation was recorded, and the session's p value for cosine directional tuning was defined as the percentage of permutations that yielded better fits (that is, higher r^{2} values) than the permutation observed in the session's DBFT curve.
Positional and directional firing rate analyses.
Spatial firing rate maps for theta cells were generated by parceling the environment into a lattice of 5 × 5 cm spatial bins and then measuring the amount of time, T_{b}, the rat spent in each bin, and counting the number of spikes, C_{b}, fired by the theta cell in each bin. The raw firing rate in each bin was then computed as R_{b} = C_{b}/T_{b}. The spatial information content of the spike train, in bits per spike, was computed from the raw firing rate map by the method of Skaggs et al. (1993). The raw firing rate map was smoothed by a single iteration of adjacent pixel averaging for plots shown in Figures 2⇓–4. For headdirection cells, directional tuning curves were computed using methods described previously (Blair and Sharp, 1995). Like DBFT curves, headdirection cell tuning curves used the convention that degrees increased in the clockwise direction.
Data inclusion criteria.
To insure adequate behavioral sampling of each movement direction, the cumulative duration of the sampled movement epochs from a session had to be ≥20 s for each movement direction, otherwise the session was excluded from the study. In addition, singleunit spike data recorded from theta cells had to meet several criteria during each session. First, to insure good singleunit isolation, theta cell spike waveforms were required to show an amplitude of ≥80 μV (peak to peak) against background noise of <30 μV, and interspike intervals had to show a refractory period ≥1 ms. Second, to insure that movement epochs contained enough spikes from which to create autocorrelograms for burst frequency analysis, theta cells were required to maintain a firing rate of >10 Hz throughout the session. Third, to assure that the neuron was a theta cell and not a thetamodulated spatially tuned neuron, theta cells had to exhibit spatial information content of <0.1 bits per spike. Fourth, for each of the eight movement directions, at least 40% of the area under the autocorrelogram's power spectrum between 4–12 Hz had to lie within a 3 Hz band centered on the power spectrum's peak. This insured that the cell's burst frequency was measurable in all movement directions, so that it would be possible to generate a meaningful DBFT curve for the cell.
Numerical methods for simulations.
Theoretical modeling simulations were performed using Matlab and NEURON (Carnevale and Hines, 2006) running on a Windows PC. Detailed descriptions of the numerical methods and source code for simulations are available on the ModelDB website (Hines et al., 2004) under accession number 129067. Briefly, spatially tuned neurons were simulated in NEURON by a single cylindrical compartment with diameter 10 μm and length 10/π μm, with passive membrane resistance and capacitance of R_{m} = 15 KΩ and C_{m} = 1.0 μF/cm^{2}, and leak reversal potential of E_{leak} = −65 mV. The simulation time step was dt = 0.1 ms. The model neuron received input from N simulated theta cell spike trains (see Results), with burst frequencies modulated by movement velocity data obtained from a recording session with a real rat. Each input spike triggered a synaptic conductance with dynamics governed by AMPA kinetics from Destexhe et al. (1994). For simplicity, excitatory (AMPA) and inhibitory (GABA_{A}) conductances were simulated using the same kinetic parameters, except that the synaptic reversal potential was E_{AMPA} = 0 mV for AMPA and E_{GABAA} = −80 mV for GABA_{A} synapses. In simulations with inhibitory input from theta cells, excitatory drive to the model neuron was provided by a voltagesensitive persistent sodium (Na_{P}) channel, which was implemented using the kinetic scheme of Uebachs et al. (2010). Active Hodgkin–Huxley kinetics were simulated using the standard NEURON mechanism “hh.mod,” with peak conductance parameters for the delayed rectifier and active sodium channels set to ḡ_{K+} = 0.005 mho/cm^{2} and ḡ_{Na+} = 0.05 mho/cm^{2}, respectively. A vector of spike times generated by the model neuron was accumulated by recording all upward crossings of a threshold at −20 mV, from which path plots and firing rate maps for the simulated spatial neurons were generated.
Results
Experimental results
Singleunit activity was recorded from anterior thalamus, hippocampus, or medial septum while rats (n = 5) foraged freely for food pellets in one of three environments (see Materials and Methods): a small cylinder (80 cm diameter), large cylinder (200 cm diameter), or small square (50 × 50 cm). A subset of neurons were classified as theta cells based on robust 6–8 Hz modulation of their spike trains and high firing rates that lacked spatial tuning (see Materials and Methods). Theta cell spike trains were analyzed to test whether their burst frequencies obeyed a “VCO frequency law,” which predicts that the theta frequency should vary as the cosine of the rat's allocentric movement direction (Burgess et al., 2005, 2007; Giocomo et al., 2007; Hasselmo et al., 2007; Burgess, 2008).
VCO frequency law
As a rat navigates across the floor of a 2D environment, its position x = (x, y) at time t can be derived by calculating the path integral of its velocity, v = dx/dt, from an initial starting location, x(0). Just as the rat's position can be obtained by computing the time integral its velocity, an oscillator's phase, ϕ, can be obtained by computing the time integral of its frequency, ω = dϕ/dt. Consequently, an oscillator with a velocitydependent frequency will have a positiondependent phase. Oscillatory interference models exploit this principle to propose that path integration is performed by VCOs with instantaneous angular frequencies, (ω_{1}, ω_{2}, …,ω_{N}), that are modulated by v as follows (Burgess et al., 2005, 2007; Giocomo et al., 2007; Hasselmo et al., 2007; Burgess, 2008): where ϕ_{n} is the nth VCO's phase in radians, Ω is a shared angular base frequency around which all VCOs are modulated, and d_{n} is a fixed “preferred vector” in the horizontal floor plane along which the nth VCO frequency is modulated by v. Here, d_{n} and v are Cartesian vectors in an allocentric coordinate system that is stationary with respect to the lab.
In Equation 11, the dot product term, d_{n} · v(t), implies that a VCO's frequency should vary as the cosine of the rat's allocentric movement direction (that is, the angle between vectors d_{n} and v) when running speed (that is, the length of vector v) is held fixed. We investigated whether theta cell burst frequencies were modulated by the rat's movement direction in this predicted manner to test whether theta cells might function as VCOs.
Directional modulation of theta cell burst frequencies
To analyze how theta cell burst frequencies were modulated by the rat's movement direction, the 360° range of directions was subdivided into eight 45° bins, and tracking data from each recording session were parsed to extract time intervals during which the rat was moving continuously in one of these eight directions (Fig. 1A). It was necessary to analyze modulation of theta frequency by movement direction in isolation from modulation by movement speed, so we devised a “speedbalancing” algorithm (see Materials and Methods, Balanced running speed distributions), which made it possible to probabilistically sample movement intervals from the session in such a way that the distribution of running speeds was rendered identical (that is, balanced) in all directions (Fig. 1B, center graph). Theta cell spike trains (Fig. 1C) from movement intervals sampled in this way were analyzed to generate speedbalanced autocorrelograms for each direction (Fig. 1D), and the power spectrum of each autocorrelogram (Fig. 1E) was taken to estimate theta cell burst frequencies in each direction.
A DBFT curve was plotted to depict how a theta cell's burst frequency varied with the rat's movement direction (Fig. 1F). A leastsquares gradient search was then performed to fit the eight points of each session's DBFT curve to a cosine function, obtained by rewriting Equation 11 as follows: where f is the burst frequency in hertz, and θ is the rat's allocentric movement direction. Before fitting, the running speed parameter, S, was set equal to the mean of the balanced running speed distribution. The fitting algorithm returned estimated values for the three parameters decorated by overbars: F̄ is the estimated base frequency (in hertz) around which the theta cell's burst frequency is modulated, and (r̄, θ̄) are polar coordinates estimating length and orientation, respectively, of the VCO's preferred vector, d_{n}.
The reformulated VCO frequency law of Equation 12 states that the amplitude of the cosine function for directional tuning (that is, the depth of frequency modulation by movement direction) depends on two parameters: r̄ and S. The value of r̄ is assumed to depend on the length of a theta cell's preferred vector, d_{n}, whereas the value of S depends on the rat's running behavior during the recording session. In our experiments, S ranged between 13 and 22 cm/s (Figs. 2⇓–4), with a mean of 16.9 ± 3.4 cm/s. Within this speed range, oscillatory interference models predict that the cosine function's amplitude, r̄S/2π, should be a few tenths of a hertz at most (see derivation below, in Distribution of predicted grid spacings). Hence, to detect the predicted directional modulation of a theta cell's burst frequency, it was necessary to obtain eight independent measurements of the burst frequency (one for each of the analyzed movement directions), all accurate to within ∼0.1 Hz. To help achieve this required precision of measurement accuracy, sessions were included in the analysis only if they met minimal inclusion criteria for behavioral sampling by the rat and theta rhythmicity of singleunit spike trains (see Materials and Methods, Data inclusion criteria). Sessions that failed to meet these criteria were excluded on the grounds that cosine directional tuning of theta burst frequencies would be undetectable even if it were present, and this limited the number of cells we could analyze.
Inclusion criteria were met by a total of 45 recording sessions obtained from 21 theta cells in five rats (data from all of these sessions are summarized in Figs. 2⇑–4). However, satisfying the minimal inclusion criteria did not guarantee that burst frequency measurements would be accurate enough to detect cosine directional tuning during all sessions in which it might have been present. To test for the presence of cosine directional tuning in each session, Equation 12 was fitted to each session's DBFT curve. A p value was then computed for the cosine fit (see Materials and Methods, DBFT curves), providing a confidence measure for the presence of cosine tuning, as well as for the accuracy of the fitted cosine tuning parameters. The degree of overlap between a session's DBFT curve and the fitted cosine function was inversely proportional to the p value computed for the fit. A significance level of p < 0.05 (the standard cutoff for twotailed tests) was beaten by 14 of 45 (31%) of the recording sessions, a much greater proportion than would be expected under the null hypothesis that theta cell burst frequencies were not directionally modulated (binomial test, p < 0.0001). A less stringent significance level of p < 0.1 (the standard cutoff for onetailed tests) was beaten by 31 of 45 (69%) of the recording sessions, which again was a much greater proportion than chance would predict (binomial test, p < 0.0001). These results indicate that among recording sessions that met criterion for inclusion in the analysis, cosinelike directional tuning was much more prevalent than would be expected by chance alone.
Having thus established the presence of directional tuning, we sought to characterize this tuning by further analyzing the data from those sessions in which a robust directional tuning signal was observed. For a session to be included in these analyses, its DBFT curve was required to exhibit a cosine fit of p < 0.1 or better. This criterion was met by at least one session for 19 of 21 (90%) of the recorded theta cells (9 of 9 cells in anterior thalamus, 5 of 6 cells in medial septum, and 5 of 6 cells in hippocampus). Whenever possible, a single theta cell was recorded repeatedly across multiple sessions. In some cases, a cell's DBFT curve beat the p < 0.1 criterion during some sessions while failing to beat criterion during other sessions. In many (but not all) of these cases, the DBFT curve resembled a cosine function even during excluded sessions that failed to beat the p < 0.1 criterion, and tuning parameters were similar to those observed for included sessions (Fig. 2A, cell 1). Such variability in the quality of the cosine fit from different sessions with the same cell may have been caused by variation of uncontrolled factors (such as the rat's movement behavior and the quality of singleunit isolation) that affected the accuracy of theta cell burst frequency estimates from one session to the next.
The proportion of sessions passing the p < 0.1 reliability criterion for cosine tuning differed among the three targeted brain structures. The criterion was beaten by 18 of 21 (86%) of the sessions from anterior thalamus (Fig. 2), 5 of 10 (50%) of the sessions from medial septum (Fig. 3), and 8 of 14 (57%) of the sessions from hippocampus (Fig. 4). This variation in the reliability of cosine tuning across brain structures does not necessarily indicate that the prevalence of cosine tuning differed among structures, because recordings obtained from different structures also varied in their video tracking resolution, session duration, and geometry of the recording environment (Figs. 2⇑–4; also see Materials and Methods), all of which are factors that can influence the accuracy of burst frequency estimates. Hence, variability in the experimental conditions may have affected the proportion of sessions beating the cosine tuning threshold, making it difficult to compare the prevalence of cosine tuning across different brain structures.
Figure 5C shows a population DBFT curve which was generated by averaging curves from all session beating the p < 0.1 threshold, first over sessions (n = 31) and then over cells (n = 19). Before averaging, individual DBFT curves were aligned with respect to their cosine tuning parameters (θ̄ on the xaxis, F̄ on the yaxis). The resulting population curve was well fit by Equation 12 (p < 0.0025), supporting the conclusion that some theta cells exhibited cosinelike directional tuning of their burst frequencies, in accordance with the VCO frequency law.
Cue control over preferred movement directions
Firing rate maps of spatially tuned neurons can follow rotations of familiar landmark cues (Muller and Kubie, 1987; Taube et al., 1990; Knierim et al., 1995, 1998; Hafting et al., 2005; Solstad et al., 2008). If spatial neurons derive their positional tuning from theta cells with directionally tuned burst frequencies, then DBFT curves of theta cells should exhibit a similar tendency to rotate with landmarks. To test this, we recorded theta cells across multiple sessions whenever possible, and sometimes rotated a landmark cue on the wall between sessions (the extramaze light source and experimenter's entry position were also rotated along with the cue card; see Materials and Methods). DBFT rotation errors (henceforth denoted by ε_{DBFT}) were quantified as the difference between rotation of θ̄ and the cue rotation angle. For 12 pairs of consecutive sessions across which a theta cell was held (Table 1), the circular mean of ε_{DBFT} was −9.1 ± 13.4° (Fig. 6A). Circular statistics indicated that ε_{DBFT} values were nonuniformly distributed (Rayleigh's Z = 7.28, p = 0.0002) and significantly clustered near zero (V = 0.77, p = 0.00002), supporting the conclusion that DBFT curves were often controlled by landmark cues.
DBFT curves did not always follow landmark cues, since ε_{DBFT} was large for some session pairs. To further study these cases, we exploited the fact that theta cells in anterior thalamus were sometimes recorded simultaneously with neighboring headdirection cells. In these cases, it was possible to measure a headdirection error (ε_{HD}) that compared rotation of the cue with rotation of the headdirection cell tuning curve (Table 1). For eight pairs of consecutive sessions across which both theta and headdirection cells were held, the circular mean of ε_{HD} was 3.2 ± 22.1°, and ε_{HD} values were clustered near zero (Rayleigh's Z = 6.2, p = 0.0004; V = 0.86; p = 0.00005). These results indicate that headdirection cell tuning curves were controlled by the landmark cue in a manner similar to DBFT curves. In cases where DBFT curves closely followed the cue, headdirection cell tuning also tended to follow the cue (Fig. 6B), but in cases where DBFT curves were poorly controlled by the cue, headdirection cell tuning also tended to be poorly controlled by the cue (Fig. 6C).
To analyze coupling between the directional tuning of theta and headdirection cells, an internal rotation error (ε_{INT}) was computed to compare rotation of theta cell DBFTs against rotation of headdirection cell tuning curves (Table 1). The circular mean of ε_{INT} was 10.9 ± 12.7°, and ε_{INT} values were clustered near zero (Rayleigh's Z = 2.8, p = 0.056; V = 0.591; p = 0.008), indicating that DBFT curves and headdirection cell tuning curves tended to rotate in tandem by similar amounts. Further supporting this, the absolute (that is, unsigned) means of ε_{HD} and ε_{DBFT} (Table 1, bottom) did not differ significantly from one another (Mann–Whitney U = 40, p = 0.92), but the unsigned mean of ε_{INT} was significantly smaller than that for ε_{DBFT} (Mann–Whitney U = 62.5, p = 0.03), implying that DBFT curves were more tightly coupled to headdirection cells than to the landmark cue; that is, theta cell DBFT curves and headdirection cell tuning curves both showed a significant tendency to follow rotations of the landmark cue (as indicated by circular statistics on each population). But both also showed incidences of failure to follow the cue, and in these instances, the directional tuning functions of theta and headdirection cells were better coupled to one another than to the landmark. A likely interpretation of these results is that rotation of the cue sometimes failed to “fool” the rat into rotating its directional reference frame with the cue, but in these cases, preferred movement directions of theta cells remained anchored to the rat's internal directional reference frame (as indicated by better coupling of theta cells to headdirection cells than the cue).
To synthesize spatial tuning functions in two dimensions, different theta cells would need to integrate the rat's movement velocity along vectors with different orientations (see below, Spatially tuned neurons in an open field). To test this prediction, a circular distribution of θ̄ values was plotted for one recording session from each theta cell (Fig. 5A). To maintain the same directional reference frame for all data points, the distribution included only one session from each cell, chosen at random from all sessions during which the cell was recorded with the cue card and extramaze light source in a “standard” position. Three cells were not recorded during a session with cues in the standard position, so for these cells, a nonstandard cue session was randomly chosen, and the cue rotation angle (−90° or +90° from standard) for the chosen session was subtracted from θ̄ before it was entered into the distribution. Values of θ̄ were uniformly distributed (Rayleigh's Z = 2.027, p = 0.13; Rao's spacing test U = 126.0, p > 0.5), supporting the conclusion that each theta cell was tuned to have its own preferred movement direction.
Amplitude of directional modulation
The fitted value of r̄ estimated the amplitude of directional frequency modulation around F̄, for the case where running speed was equal to S. If theta cells function as VCOs, then r̄ estimates the length of the VCO's preferred movement vector. Empirical values of r̄ can be tested against theoretical predictions by recognizing that when a grid cell is formed by interference among VCOs, the vertex spacing of the grid is inversely proportional to VCO vector length (Burgess et al., 2007; Giocomo et al., 2007; Hasselmo et al., 2007). For the case of an equilateral grid, the smallest obtainable vertex spacing is formed by combining three VCOs with preferred movement vectors of the same length, r̄ = d_{1} = d_{2} = d_{3}, and differing orientations that are 120° apart from each other. To compute the spacing of a grid cell formed from such a triplet of VCOs, we may express the phases of the three oscillators (k = 1, 2, 3) as follows:
where Φ is a common “reference phase” shared among all VCOs (Eq. 17). By the VCO envelope formula (Eq. 20), we obtain the following:
where the three vectors given by D_{1} = d_{1} − d_{2}, D_{2} = d_{2} − d_{3}, and D_{3} = d_{3} − d_{1} have identical length D = D_{1} = D_{2} = D_{3}, and their directions are also 120° apart. Thus, the envelope can be regarded as a sum of three straight cosine gratings with the spacing λ = 2π/D. The spacing Λ̄ of the hexagonal grid formed by these three cosine gratings is as follows:
Based on the geometry, we have D =
It should be emphasized once again that these calculations assume grids are formed from triplets of theta cells with preferred direction vectors that are 120° apart. But there are many different ways to form a hexagonal grid from interference among VCOs (Burgess, 2008), and the slope of the inverse relationship between vector length and vertex spacing depends on exactly how VCOs are combined form the grid. For example, adding six VCOs with preferred directions that are 60° apart would lead to a predicted grid spacing that is larger than that in Equation 16 by a factor of
Modulation of burst frequencies by running speed
The fitted value of F̄ estimated the base frequency around which theta cell bursting was modulated. Values of F̄ showed variability among different rats across an approximate range of 6–8 Hz, and showed a lesser degree of variability (up to ∼1 Hz) across different sessions within the same rat (Fig. 4B, xaxis). This pattern of results is consistent with prior studies showing variability of theta frequency across rats and sessions (Jeewajee et al., 2008b). If the base frequency remains rigidly constant at all times and does not vary with running speed, then Equation 11 dictates that a VCO's frequency should increase with running speed when the rat moves in the VCO's preferred direction, and decrease with running speed when the rat moves in the opposite (antipreferred) direction. Averaged across all directions, these opposing influences should cancel each other out, yielding a prediction that the mean VCO frequency would be identical for all running speeds. Contradicting this prediction, prior studies have shown that the theta frequency increases with running speed when movement direction is disregarded (Rivas et al., 1996; Geisler et al., 2007; Jeewajee et al., 2008a). One simple way to account for these results is to assume that the VCO base frequency [Eq. 11, Ω(t)] increases in proportion with running speed (Burgess, 2008). As long as this speeddependent base frequency remains identical for all VCOs, then path integration would not be impaired.
To investigate modulation of the base frequency by running speed, an analysis was conducted in which theta cell burst frequencies were analyzed as a function of the rat's running speed in directions parallel, antiparallel, and orthogonal to each cell's preferred movement direction (Fig. 5D). Averaging results across sessions and then cells, it was found that theta cell burst frequencies increased with running speed in all directions, and the slope of speed modulation was steepest for the preferred and shallowest for the antipreferred direction. The observation that the theta frequency increased with running speed in all directions (rather than increasing in the preferred and decreasing in the antipreferred direction) suggests that the base frequency was not constant, but instead increased with running speed. As long as Ω(t) varies only with running speed, and not with movement direction, then the VCO frequency law still dictates that theta burst frequencies should vary as the cosine of the rat's movement direction, as we have observed.
Analytic simulations of spatially tuned neurons
Experimental results presented above suggest that VCO signals may be encoded by spike trains of rhythmically bursting theta cells. Computer simulations were performed to test whether spatially tuned neurons could be formed by phase interference among theta cells with firing properties similar to those observed in our experiments (source code and input data for simulations are available on ModelDB). Before describing these simulation results, we first introduce an analytic expression for the spatial envelope function that is synthesized from phase interference among any arbitrary set of VCOs.
VCO envelope and carrier equations
Equation 11 dictates that the nth VCO's frequency varies linearly with the rat's velocity along a preferred vector, d_{n}, so that a component of the VCO's phase becomes dependent on the rat's position along that same vector. To isolate this positiondependent phase component, we may separately integrate the two summed terms in Equation 11, and thereby express the instantaneous phase of the VCO as a sum of two terms:
where Φ is a common “reference phase” shared among all VCOs [obtained by integrating the base frequency, Ω, between times 0 and t from an initial starting phase, Φ(0)], and δ_{n} is the nth VCO's unique offset from the reference phase, which depends upon the rat's spatial position, x(t), as follows:
with δ_{n}(x(0)) = ϕ_{n}(0) − Ω(0) denoting the nth VCO's initial offset from the reference phase. It follows from Equation 18 that if a rat navigates to an arbitrary position, x(t), then the nth VCO's phase at that position will be offset from the reference phase by an amount that depends strictly upon the distance of x(t) from the rat's starting position, x(0), along the VCO's preferred vector, d_{n}. Moreover, if all VCOs share the same reference phase, then their phase offsets from each other will also depend strictly on the rat's location. Consequently, every location in space will be encoded by a specific pattern of phase relationships among VCOs. If VCO phase relationships depend strictly on the rat's position, then phase interference among VCOs is guaranteed to synthesize an envelope function that “sits still” in space, regardless of the rat's movement trajectory through the environment (Burgess et al., 2005; 2007). The outputs from N VCOs may be thus summed together to generate an interference pattern as follows (Fig. 7A):
where w_{n} is a weighting coefficient for the nth VCO. The second step in Equation 19 states that the sum of VCO outputs can be rewritten as a product of just two terms: a carrier signal C(t) that varies over time, and an envelope signal E(x) that varies over space. By exploiting a trick based on Euler's formula l (Hartmann, 1998), E(x) can be expressed as a sum of complex exponential terms by the following VCO envelope equation:
where i =
Equation 19 states that the spatial envelope signal is multiplicatively modulated by a temporal carrier signal, assumed here (as in prior oscillatory interference models) to account for theta modulation of spatial neuron spike trains. The carrier signal can be written as follows: where with and
Analytic model of spatial neurons formed from theta cells
The VCO envelope equation implies that any neuron which can detect synchrony among VCOs (by firing selectively when their phases converge) is guaranteed to “pick out” a selected spatial region of the rat's environment (by firing only when the rat is located within that region). Prior interference models have demonstrated this for the special case of grid cells, showing how phase interference among VCOs can synthesize envelopes that form a hexagonal lattice tiling the environment (Burgess et al., 2005, 2007; Giocomo et al., 2007; Hasselmo et al., 2007; Burgess, 2008; Zilli and Hasselmo, 2010). Here, we expand on these results by showing that, more generally, the VCO envelope can mimic almost any spatial function and thereby simulate the firing rate maps of almost any spatially tuned neuron (including place and border cells, in addition to grid cells).
A simplified case is illustrated in Figure 7A, which plots output from eight VCOs (modeled as cosine functions) while a simulated rat runs at a constant speed of 25 cm/s along the length of a 1 m linear track. In this example, all of the VCOs have different fixed frequencies, and they converge in phase at the center of the track. If a postsynaptic neuron sums input from just two of these eight VCOs, then its output can generate a beat envelope that waxes and wanes at fixed intervals along the length of the track, mimicking the spatially periodic firing pattern of a grid cell as in prior models (Fig. 7A, E_{grid}). But if the target neuron sums input from all eight of the VCOs, then its output will grow large only at the center of the track where all VCO phases converge, mimicking the unitary firing field of a place cell (Fig. 7A, E_{place}). Thus, detection of phase convergence among few VCOs (a common event that recurs at fixed intervals along the track) simulates a grid cell, whereas detection of phase convergence among many VCOs (a rare event which occurs at just one location on the track) simulates a place cell.
More generally, a postsynaptic neuron can generate almost any spatial function by detecting phase convergence among a properly chosen set of VCOs. To see why this is so, it is helpful to rewrite Equation 20 in the following form: where A and B are given by Equations 23 and 24, and x is now assumed to be a vector representing the rat's location in a 2D environment instead of a linear track (note that generalization to arbitrary dimensions is possible; see Discussion). Here, the envelope equation becomes a sum of 2D cosine functions (or gratings) formed by interferences between all possible pairs of the N VCOs, indexed by k and l. These cosine gratings define the components of the envelope's 2D spatial frequency spectrum, and they are fully determined by the preferred vectors (d_{1}, d_{2}, …, d_{N}), weighting coefficients (w_{1}, w_{2}, …, w_{N}), and starting phases [ϕ_{1}(0), ϕ_{2}(0), …, ϕ_{N}(0)] of the VCOs that interfere with one another to synthesize the envelope. If these parameters are properly chosen, then a sufficiently large number of VCOs should be able to approximate the spatial frequency spectrum of almost any desired spatial function.
Let us suppose that the rat brain contains a large and diverse population of VCOs, all sharing a common reference phase as assumed by Equations 17 through 25, and that each VCO signal is generated by a rhythmically bursting theta cell (as suggested by our experimental findings above). Such a population of theta cells can encode a basis set from which arbitrary spatial tuning functions could be synthesized. A different spatial envelope would be formed by interference among any unique subset of theta cells in the population, and a few key parameters would determine the frequency spectra (and thus the shapes) of these envelope functions. As noted above, one such parameter is the temporal phase of each VCO signal. Hence, to generate a diversity of different spatially tuned neurons, it would be useful for VCOs to be implemented as multiphase oscillators, which generate multiple copies of the same VCO signal at diverse phase shifts. Supporting this possibility, there appears to be a continuous gradient of theta phase along the septotemporal axis of the hippocampus (Lubenov and Siapas, 2009).
This phase diversity requirement entices us to conceptualize the CPG circuit for theta rhythm as “ring oscillator” composed from a circular layer of theta cells (Blair et al., 2008). A localized “bump” of activity can circulate around the ring at the theta frequency, so that each cell in the ring generates rhythmic theta bursts on a different phase of the same VCO cycle (Fig. 7B). To regulate the VCO frequency in accordance with Equation 11, the ring may receive an external driving input that encodes the rat's translational velocity along the ring's preferred vector, d_{n}. In a population of ring oscillators with different velocity inputs, the nth ring oscillator serves as a multiphase VCO that performs translational path integration along its own preferred vector, d_{n}. This ring oscillator model of translational path integration is quite similar to hypothetical ring attractor models of angular path integration by headdirection cells (Skaggs et al., 1995; Redish et al., 1996; Zhang, 1996; Song and Wang, 2005), raising the possibility that translational and angular path integration might both be performed by similar ring attractor circuits (Blair et al., 2008; Mhatre et al., 2010).
Figure 7C illustrates a hypothetical arrangement of ring oscillators within a structured matrix of VCOs that contains a distributed population of theta cells from which many different spatial tuning functions can be synthesized. The VCO parameters of each theta cell in Figure 7C are determined by the cell's position in the matrix; theta cells residing in the same ring share the same preferred vector, and thus generate the same VCO signal at different phase shifts. Theta cells residing in different rings have different preferred vectors, and thus generate different VCO frequencies. The matrix is organized so that rings in the same row are modulated by running speed with the same slope (and thus have VCO vectors of identical length), whereas rings in the same column prefer movement in the same direction (and thus have VCO vectors of identical orientation).
A target neuron that sums input from theta cells in the VCO matrix can simulate the positional firing properties of almost any spatially tuned neuron, depending on the rows, columns, and ring positions of the theta cells that provide its input. This is illustrated on the right side of Fig. 7C, which shows analytic envelope functions that were simulated by solving Equation 20. To simplify the derivation of spatial envelopes in these simulations, the rat's initial starting position, x(0), was set equal the origin of the spatial coordinate system, which corresponded to the center of the arena. In addition, the starting phases, ϕ_{n}(0), for all VCOs were initialized to zero, and all VCO weighting coefficients, w_{n}, were set uniformly equal to 1. Spatial tuning functions were simulated by assigning preferred vectors, [d_{1}, d_{2}, …, d_{N}], to a set of N VCOs and then solving Equation 20 at each point on a square lattice representing positions in a 1 m^{2} square environment. The amplitude of E(x) was converted to a “firing rate” by filtering through a spike threshold: where R(x) is the firing rate at x, H is the Heaviside function, and κ represents the spike threshold which was set to κ = 0.7 × max[E(x)] for all simulations presented here.
As in prior oscillatory interference models (Burgess et al., 2005, 2007; Giocomo et al., 2007; Hasselmo et al., 2007; Burgess, 2008), a target neuron that receives input from properly chosen theta cells can simulate a grid cell with various vertex spacings (blue and green lines in Fig. 7C). Here, the vertex spacing depends on which matrix row theta cells reside in, and the grid orientation depends on which columns theta cells reside in. But in addition to grid cells, other types of spatially tuned neurons can also be formed. For example, suppose that a target neuron's input comes from theta cells residing in different rows of the same column (Fig. 7C, red lines). Since all of these theta cells have preferred vectors of the same orientation (but with different slopes of speed modulation), the target neuron fires along a linear band oriented perpendicular to the preferred vector orientation. This linear band can mimic the firing rate map of a border cell that fires along a straight wall (Fig. 7C, bottom right) or of a boundary vector cell (Lever et al., 2009) that fires at a fixed distance from the wall (data not shown). Alternatively, if the target neuron's input comes from theta cells in every column of a single row, then it can simulate a place cell (Fig. 7C, top right). Like the vertex spacing of a grid cell, the size of a place cell's firing field would depend on which row of the matrix provides input to the target neuron.
As a general rule, a target neuron that receives input from many theta cells will tend to fire in very few (and possibly zero) subregions within a the circumscribed area of any given environment (since phase convergence among many theta cells is a rare event that occurs at very few locations); the shape and size of this subregion will depend on exactly which theta cells provide the neuron's input. Conversely, a target neuron that receives input from just a few theta cells will tend to fire within multiple subregions of an environment (since phase convergence among a small number of theta cells is a common event that recurs at multiple locations). For example, we have already seen that grid cell firing patterns can emerge from interference among a few systematically chosen theta cells. Alternatively, interference among randomly chosen theta cells can mimic the firing patterns of neurons that fire at multiple locations that are randomly distributed (results not shown), such as those found in the dentate gyrus (Leutgeb et al., 2007).
Biophysical simulations of spatially tuned neurons
In Figure 7, spatial neurons were simulated by analytically solving Equation 20 to generate envelope functions. Hence, the envelopes were obtained by linearly summing a basis set of idealized VCOs, each represented as a velocitymodulated cosine function consisting of a single fundamental harmonic component. But our experimental data suggest that VCO signals may be encoded by noisy theta cell spike trains. Unlike idealized cosine VCOs, these noisy spike trains contain many complex harmonics in addition to the “fundamental” theta frequency of their bursting rate. There is thus no guarantee that spike trains from biological theta cells could interfere with one another to synthesize spatial envelopes, because the additional harmonics introduced into the VCO signals by noise and the spiking process—which are not accounted for by the VCO envelope equations—might interfere with one another in ways that compromise the spatial rigidity and trajectory invariance of the envelopes. Moreover, a biological neuron typically does not compute a perfect linear sum of its inputs, as assumed in the analytic simulations of Figure 7.
Simulating theta cell spike trains
To investigate whether spatially tuned neurons can be formed from spiking theta cells similar to those we observed in our experiments, we simulated theta cell spike trains using a Poisson process. A theta cell spike was generated pseudorandomly at each simulation time step (dt = 2 ms) with an oscillating probability that varied with the amplitude of an ideal cosine VCO (Fig. 8A). Experiments indicated that the base frequency of theta cell VCOs increased with running speed (Fig. 5D). To incorporate this result into simulations of theta cell spike trains, the following version of the VCO frequency law was used: where dϕ/dt is the instantaneous angular frequency of the oscillating spike probability, S(t) is the rat's instantaneous running speed, and Ω_{0} + S(t)/40 defines the speeddependent base frequency, so that Ω(t) increases with running speed at a slope of 1/40 Hz/cm/s. To generate simulated spike trains, movement velocity data from a 60 min segment of an experimental recording session were upsampled to 500 Hz (from the video tracking rate of 30 Hz) to match the simulation time resolution of dt = 2 ms. At each time step, the VCO phase, ϕ(t), was computed by numerically integrating the VCO frequency defined by Equation 27 across the 60 min session, with the values for movement velocity derived from the tracking data. Spike probability at each time step was then computed as follows: where H is the Heaviside function. Thresholding by H gives a spike probability of P(t) = 0 in an interval of width π/2 centered upon the valley of each theta cycle (so that the simulated theta cell falls silent between bursts on each cycle). Outside of this silent interval, P(t) varies as a continuous function of the VCO amplitude, reaching a maximum value of 0.5 spikes per dt at the peak of each theta cycle. When Ω_{0} was set to 7 Hz (as in simulations presented here), Equation 28 produced simulated theta cell spike trains with a mean spiking rate of ∼100 Hz (note that a refractory period of 2 ms between spikes was automatically conferred by the 2 ms time resolution of the simulations). To obtain firing rates lower than 100 Hz, spikes were randomly selected for deletion from the train with a probability inversely proportional to the desired firing rate. When simulated spike trains were analyzed using the same algorithms that were used for experimental data analysis, simulated theta cells exhibited directional tuning of their burst frequencies similar to that observed for real theta cells (Fig. 8).
Place cell synthesized from theta cells on a linear track
Based on existing knowledge of hippocampal and thalamocortical circuits, it is reasonable to assume that many of the theta cells we recorded in medial septum (Fig. 3) and hippocampus (Fig. 4) were inhibitory neurons (Ranck, 1973; Freund and Antal, 1988), whereas theta cells in the anterior thalamus (Fig. 2) may have been excitatory thalamocortical projection neurons (Jones, 1985; Tsanov et al., 2011b). We thus conducted simulations to test whether spatial tuning functions could be formed by inhibitory versus excitatory inputs to a postsynaptic target neuron from theta cells.
To compare interference among excitatory versus inhibitory theta cells, we replicated the place cell simulation of Figure 7A, in which a rat runs at constant speed across a linear track. Each of the eight cosine VCO signals in Figure 7A was converted into a theta cell spike train by the method illustrated in Figure 8. The NEURON simulation environment (Carnevale and Hines, 2006) was then used to create a singlecompartment model of a postsynaptic target cell, which received identically weighted inputs (either excitatory or inhibitory) from the eight simulated theta cell spike trains.
In simulations where theta cell spike trains triggered IPSPs, the target neuron was able to fire selectively in the center of the track by detecting locationdependent synchrony among its inputs from theta cells (Fig. 9A). It may seem counterintuitive that synchronization among inhibitory inputs would excite the postsynaptic neuron, but similar results were also reported by Zilli and Hasselmo (2010) in simulations of grid cells. This paradoxical excitation of the postsynaptic cell by synchronized inhibitory inputs occurs because desynchronized input from inhibitory theta cells produces a constant, uninterrupted barrage of inhibition that prevents the target neuron from firing; when inhibitory theta bursts become synchronized, prolonged relief from inhibition (lasting tens of milliseconds) occurs during gaps between synchronized theta bursts. Under the influence of an excitatory drive, which was provided in our simulations by a persistent sodium (Na_{P}) current (Fig. 9A, blue trace), the target neuron can fire reliably during the disinhibitory gaps. Hence, by detecting coincident gaps between synchronized theta bursts, the target neuron can measure locationdependent burst synchrony among its inputs from inhibitory theta cells, so that its firing rate approximates the spatial envelope defined by Equations 20 and 25.
We next conducted simulations in which the postsynaptic neuron was excited (rather than inhibited) by inputs from theta cells. To detect synchrony among excitatory theta cell inputs, it is necessary for the postsynaptic neuron to sum EPSPs in such a way that its membrane becomes depolarized above spike threshold only when theta inputs are synchronized, but not when they are desynchronized. When EPSPs were triggered by the same eight theta cell spike trains that were used for inhibitory simulations in Figure 9A, the target neuron consistently failed to detect locationdependent synchrony among the theta cells (and thus did not fire selectively in the center of the track), and this failure was replicated over a wide range of parameter values for synaptic and membrane conductances (results not shown).
Synchrony detection by EPSP summation failed because the VCO signals were “undersampled” by noisy theta cell spike trains; that is, since the timing of theta cell spikes was noisy, alignment of burst phases on the slow time scale of the theta period did not necessarily cause individual input spikes to become better aligned on the faster time scale of EPSP summation. One way to overcome this undersampling problem would be to reduce the level of noise in the spike times [for example, as in the study by Hasselmo (2008)], so that alignment of burst phases would cause better alignment of individual spikes (and thus better summation of EPSPs). However, such lownoise spike trains would no longer resemble the noisy spike trains of theta cells recorded in our experiments. Another possibility would be to increase the number of spikes generated per theta cycle (Zilli and Hasselmo, 2010), and thereby increase the density at which each VCO signal is sampled by theta cell spikes. However, the mean firing rate of biological theta cells rarely exceeds ∼80 Hz, and this imposes an upper bound on the number of spikes that a single theta cell can generate on each theta cycle. In our model, simulated theta cells fired at a mean rate of about 40 Hz, similar to the mean firing rate of theta cells we recorded in anterior thalamus (Fig. 2). To increase the number of spikes per theta cycle without exceeding biological limits on the firing rates of theta cells, we conducted simulations in which eight VCOs were each redundantly encoded by four excitatory theta cells with noisy spike trains at mean firing rates of ∼40 Hz (Fig. 9B). Thus, the target neuron received a total of 32 excitatory inputs, from eight VCOs encoded by four theta cells each. In these simulations, the target neuron was able to fire selectively in the center of the track by detecting locationdependent synchrony among excitatory inputs from theta cells (Fig. 9B).
Together, simulation results in Figure 9, A and B, indicate that spatially selective firing can be generated by a postsynaptic neuron that detects locationdependent synchrony among either inhibitory or excitatory theta cells with firing properties similar to those observed in our experiments. However, a larger number of theta cell inputs may be required to achieve robust spatial tuning in the excitatory than in the inhibitory case. This is because in the excitatory case, the postsynaptic neuron must perform the difficult task of detecting coincidences among EPSPs (which are brief events lasting only a few milliseconds), whereas in the inhibitory case, the postsynaptic neuron performs the easier task of detecting coincidences among the gaps between theta bursts (which are more prolonged events lasting tens of milliseconds).
Spatially tuned neurons in an open field
To simulate spatially tuned neurons in an open field environment, theta cell VCOs were modulated by a movement velocity signal derived from the video tracking data for a 60 min segment of an experimental recording session in the 80 cm cylinder. Each simulated theta cell's burst frequency was modulated by the rat's movement velocity along a preferred vector, d_{n}. The singlecompartment target neuron received identically weighted inhibitory inputs from N simulated theta cell spike trains, where N differed depending upon the simulation (Fig. 9). By assigning appropriate values for the phases and preferred vectors of the theta cell VCOs, it was possible for the target neuron to simulate firing rate maps of grid, place, or border cells (Fig. 9C; source code with parameter values available on ModelDB). These results demonstrate that many different kinds of spatially tuned neurons can be formed by detecting burst synchrony among inputs from theta cells with firing properties similar to those observed in our experiments.
Phase precession
In linear track simulations, place cell bursting occurred in phase with synchronized excitatory theta cell inputs (Fig. 9B) or in antiphase with synchronized inhibitory theta cells inputs (Fig. 9A). In both cases, place cell spikes exhibited phase precession against the shared VCO reference phase, Φ(t) (Fig. 9A,B, bottom). This is similar to the manner in which real place cell spikes (O'Keefe and Recce, 1993) and membrane oscillations (Harvey et al., 2009) exhibit phase precession against the hippocampal local field potential (LFP) in behaving rodents.
It is not fully understood what oscillatory signal is actually measured by the LFP in experimental recordings (Blair et al., 2008; Burgess, 2008; Geisler et al., 2010), but we equated the common phase factor Φ(t) with the LFP since it is the only phase component that is shared among all theta oscillators in the model, and might therefore dominate the field EEG in physiological recordings. Under this assumption, spikes generated by spatially tuned neurons can only exhibit phase precession under certain conditions. To define these conditions, let us first suppose that the phase offset term ϕ (Eq. 22) is a constant and does not change in time (that is, dϕ/dt = 0), so that the carrier signal defined by Equation 21 becomes identical to cos Φ(t), up to a constant offset. In this case, the burst frequency of the target cell (such as a place cell) would be identical to the LFP frequency, so there would be neither phase precession nor phase procession. But if ϕ increases in time (that is, dϕ/dt > 0), then the carrier frequency becomes faster than the LFP frequency, and we should see theta phase precession. More generally, it follows from Equation 22 that To simplify this expression without loss of generality, we may assume that the animal is running in the positive x direction in the lab coordinate system. Then Equation 18 reduces to the following scalar form: where d_{n,x} is the x component of the preferred vector d_{n} for theta modulation. Thus we have the following: where dx/dt = v is the animal's running speed.
Now suppose the animal is momentarily located at the peak of the envelope (e.g., place field). The maximum peak value is possible only when the phases of all the oscillators are aligned. Without loss of generality, we write this condition as δ_{n} = 0 for all n. Now Equations 23 and 24 are reduced simply to
Discussion
Here we recorded theta cells with burst frequencies that varied as the cosine of a rat's allocentric movement direction (Figs. 1⇑⇑⇑⇑–6), and thus appeared to behave like VCOs in oscillatory interference models of path integration (Burgess et al., 2005, 2007; Giocomo et al., 2007; Hasselmo et al., 2007; Blair et al., 2008; Burgess, 2008; Hasselmo, 2008; Zilli and Hasselmo, 2010). Computer simulations demonstrated how grid, place, and border cells could all be formed by a postsynaptic neuron that detects locationdependent synchrony among inputs from theta cells with firing properties similar to those observed in our experiments (Figs. 7⇑–9). These results suggest that populations of theta cells might generate VCO signals from which spatial memories can be synthesized via oscillatory interference.
Theta cells as VCOs
It has been proposed that theta oscillations encoding VCO signals might arise from intrinsic membrane currents within the dendrites of spatially tuned neurons (Kamondi et al., 1998; Lengyel et al., 2003; Burgess et al., 2007; Giocomo et al., 2007; Hasselmo et al., 2007; Losonczy et al. 2010) (but see Remme et al., 2010) or from populations of neurons that fire rhythmically at velocitymodulated frequencies (Blair et al., 2008; Burgess, 2008; Hasselmo, 2008; Zilli et al., 2009; Zilli and Hasselmo, 2010). In the latter case, theory predicts that rhythmic firing should be modulated by a rat's movement direction. It has been shown previously that anterior thalamus contains neurons with firing rates that are modulated by both head direction and theta rhythm (Tsanov et al., 2011a) and that medial septum contains neurons that exhibit maximal theta rhythmicity when a rat runs in a specific direction on a linear track (King et al., 1998). However, our present findings provide the first evidence for theta oscillations that with frequencies that vary as the cosine of a rat's allocentric movement direction, which is absolutely necessary for path integration by oscillatory interference.
Cosine directional tuning was observed in a higher proportion of recording sessions from anterior thalamus (Fig. 2) than medial septum (Fig. 3) or hippocampus (Fig. 4), which may have been a consequence of variable experimental conditions that affected the accuracy of burst frequency measurements, rather than an indicator of functional differences among theta cells recorded in different brain regions (see Results). At the running speeds we observed here, theory predicts that VCO frequencies should vary with movement direction by only fractions of a hertz. To measure theta cell burst frequencies with this level of accuracy, we applied strict criteria for inclusion of session data in the analysis, which unavoidably biased the sample of theta cells we analyzed. Thus, from our present findings, it is difficult to estimate the overall proportion of thetamodulated cells that exhibits directional tuning of their burst frequencies. However, directional modulation of theta cell burst frequencies appeared to be quite prevalent in the structures we targeted, since our reliability threshold for cosine directional tuning was beaten during at least one session by 90% (19 of 21) of the cells admitted for analysis. Moreover, directional tuning followed landmark rotations (Fig. 6), as would be expected for neurons involved in path integration.
Spatial envelopes can only be synthesized from a set of VCOs that are modulated around exactly the same base frequency, and thereby coupled to one another through a shared reference phase (Eq. 17). To test whether different theta cells share a common reference phase, it would be necessary to simultaneously record pairs of theta cells (which was not done here) and show that the offset between their burst phases varies with the rat's position in a manner that forms a 2D cosine grating (Burgess, 2008). Such experiments could provide a critical test of whether theta cells encode VCO signals that share a common reference phase. However, they would be challenging to perform, because only theta cells with sufficiently different preferred vectors (separated by sufficiently long difference vectors) would form gratings with spatial frequencies high enough to observe within the confines of an experimental environment.
Spatial memory coding by theta cells
Computational models have assumed that place cells form their firing fields by summing inputs from grid cells (McNaughton et al., 2006; Solstad et al., 2006; Blair et al., 2007; Gorchetchnikov and Grossberg, 2007; Savelli and Knierim, 2010; Monaco and Abbott, 2011) or boundary vector cells (Hartley et al., 2000). Here, we showed how place cells might instead form by summing direct excitatory or inhibitory inputs from theta cell VCOs (Figs. 7, 9). This does not necessarily conflict with prior models, since grid and boundary vector cells can themselves be synthesized from theta cells (Fig. 9C). Hence, place cells could either receive direct inputs from theta cell VCOs (as in our simulations) or indirect inputs from grid, border, or boundary vector cells (as in prior models) that are each formed by “bundling” several VCOs together. In either case, a place cell which fires over a small region of space would ultimately derive its positional selectivity by detecting rare occurrences of locationdependent synchrony among many VCOs.
Recent studies reported that hexagonal firing patterns of entorhinal grid cells were severely impaired when theta rhythm was disrupted by inactivation of the medial septum (Brandon et al., 2011; Koenig et al., 2011), supporting the possibility that grid cells may derive their positional firing from theta oscillations. However, hippocampal place cells can continue to exhibit spatial firing (and retain some theta modulation of their spike trains) after disruption of ascending theta pathways (Mizumori et al., 1989; Sharp and Koester, 2008; Koenig et al., 2011), after loss of inputs from entorhinal cortex (Van Cauter et al., 2008; Brun et al., 2008a), and before the development of entorhinal grid cells in rat pups (Wills et al., 2010). One explanation for this pattern of findings could be that place cells receive diverse VCO inputs from both direct (e.g., theta cells) and indirect (e.g., grid or boundary vector cells) sources, and any of these inputs may be sufficient to support some degree of spatial firing by place cells after other inputs are lost. Moreover, some inputs may be less susceptible than others to disruptions of theta rhythm. For example, boundary vector cells might derive some of their spatial tuning from sensory signals generated by encounters with landmarks or boundaries, rather than by detecting synchrony among VCOs that perform path integration (O'Keefe and Burgess, 2005; Burgess et al., 2007; Wills et al. 2010).
Anatomical origins of VCO signals
To synthesize diverse populations of spatially tuned neurons from VCOs, identical copies of each VCO signal must be generated at different temporal phase shifts. To satisfy this computational requirement, we conceptualized the VCO circuit as a “ring oscillator” composed from a circular layer of theta cells, each bursting on its own phase of the VCO cycle (Fig. 7B) (Blair et al., 2008). In simulations, we showed how place, grid, and border cells could all be formed by target neurons that detect synchrony among inputs from different theta cells residing in a matrix of such ring oscillators (Fig. 7C).
Evidence suggests that theta oscillations in hippocampus and cortex may be driven and synchronized by subcortical CPGs, which send ascending projections to spatial memory circuits via the medial septum and limbic thalamus (Petsche et al., 1965; Gogolak et al., 1967; Vertes and Kocsis, 1997). These ascending theta pathways exhibit considerable anatomical overlap with ascending headdirection pathways (Vann and Aggleton, 2004; Taube, 2007; Blair et al., 2008; Sharp and Koester, 2008). Interestingly, our ring oscillator model of the theta CPG circuit performs translational path integration in almost exactly the same way that angular path integration is performed by ring attractor models of headdirection cells (Skaggs et al., 1995; Redish et al., 1996; Zhang, 1996; Song and Wang, 2005; Blair et al., 2008). This suggests the intriguing possibility that headdirection signals and velocitymodulated theta rhythm might both emerge from similar subcortical ring attractor circuits and then follow parallel ascending pathways to reach spatial memory circuits in the hippocampus and cortex. To test this, future studies should investigate whether theta cells with directionally tuned burst frequencies can be identified in subcortical structures from which theta rhythm originates.
Neural coding by envelope synthesis
The nervous system generates a multitude of oscillatory rhythms spanning a wide range of frequency bands (Buzsáki, 2006), and phase relationships among such oscillators have been proposed to encode information not only about space, but also about time (Miall, 1989; Hopfield and Brody, 2001; Matell and Meck 2004) or complex sensory stimuli (Laurent et al., 2001; Tiesinga et al., 2008; Koepsell et al., 2010). If the phase of a neural oscillator shifts smoothly in conjunction with a timevarying state variable, x, then that oscillator's frequency would inherently be modulated by dx/dt, and thereby exhibit firing properties analogous to the velocitymodulated theta cells we observed here. In principle, the existence of such oscillators in the nervous system would make it possible to synthesize almost any desired envelope function (as specified by Eqs. 20, 25) in the state space of x, and thereby provide a powerful and flexible mechanism for constructing mental representations in biological neural networks.
Footnotes

This work was supported by National Science Foundation/NIMH Collaborative Research in Computational Neuroscience Program Grant MH079511 awarded to H.T.B. and K.Z. We thank Neil Burgess and James Knierim for helpful discussions.
 Correspondence should be addressed to Dr. Hugh T. Blair, Psychology Department, 1285 Franz Hall, University of California, Los Angeles, Los Angeles, CA 900951563. blair{at}psych.ucla.edu