## Abstract

Hippocampal place cells convey spatial information through a combination of spatially selective firing and theta phase precession. The way in which this information influences regions like the subiculum that receive input from the hippocampus remains unclear. The subiculum receives direct inputs from area CA1 of the hippocampus and sends divergent output projections to many other parts of the brain, so we examined the firing patterns of rat subicular neurons. We found a substantial transformation in the subicular code for space from sparse to dense firing rate representations along a proximal-distal anatomical gradient: neurons in the proximal subiculum are more similar to canonical, sparsely firing hippocampal place cells, whereas neurons in the distal subiculum have higher firing rates and more distributed spatial firing patterns. Using information theory, we found that the more distributed spatial representation in the subiculum carries, on average, more information about spatial location and context than the sparse spatial representation in CA1. Remarkably, despite the disparate firing rate properties of subicular neurons, we found that neurons at all proximal-distal locations exhibit robust theta phase precession, with similar spiking oscillation frequencies as neurons in area CA1. Our findings suggest that the subiculum is specialized to compress sparse hippocampal spatial codes into highly informative distributed codes suitable for efficient communication to other brain regions. Moreover, despite this substantial compression, the subiculum maintains finer scale temporal properties that may allow it to participate in oscillatory phase coding and spike timing-dependent plasticity in coordination with other regions of the hippocampal circuit.

## Introduction

Hippocampal place cells represent an animal's location in the world through a combination of rate and phase coding. These neurons fire selectively in certain regions of an environment (“place fields”) and are otherwise nearly silent. The resulting sparse pattern of ensemble activity among hippocampal place cells can be used to decode an animal's path through an environment (Wilson and McNaughton, 1993; Brown et al., 1998; Huxter et al., 2008) and distinguishes between different environments (Leutgeb et al., 2005; Karlsson and Frank, 2009). The phase of firing relative to ongoing theta oscillations (5–10 Hz) in the local field potential (LFP) also conveys spatial information. A hippocampal place cell will fire at progressively earlier phases of theta as the animal passes through its place field, a phenomenon known as theta phase precession (O'Keefe and Recce, 1993; Skaggs et al., 1996).

Hippocampal place cells are part of a larger circuit for spatial representation that includes the entorhinal cortex (EC) and subiculum (Sharp, 1999; Barry et al., 2006; Moser et al., 2008). While recent investigations have characterized rate and phase coding in the EC (Hafting et al., 2008; Mizuseki et al., 2009), the subiculum has received relatively little scrutiny. The subiculum receives exceptionally dense projections from area CA1 of the hippocampus (Amaral et al., 1991; Cenquizca and Swanson, 2007; Jinno et al., 2007; Fuentealba et al., 2008), and like area CA1, has topographically organized reciprocal connections with the EC (Naber et al., 2001; Kloosterman et al., 2004). The subiculum is also a major output structure of the hippocampal circuit, sending projections to diverse targets including the prefrontal cortex, amygdala, nucleus accumbens, and hypothalamus (Witter, 2006).

Previous studies reported that neurons in the subiculum tend to fire in multiple irregularly spaced firing fields and are less spatially selective than neurons in the hippocampus (Barnes et al., 1990; Sharp and Green, 1994). Some evidence suggests that subicular neurons maintain invariant spatial firing fields across distinct environments (Sharp, 2006), unlike hippocampal neurons which tend to “remap” when the environment is changed (Muller and Kubie, 1987; Best et al., 2001). Similarly, a recent study has proposed that some neurons in the subiculum are “boundary-vector cells” that are insensitive to environmental context (Lever et al., 2009). These findings suggest, quite surprisingly, that the hippocampal representation of specific locations and environments may not be well maintained in the subiculum, just one synapse downstream from CA1 place cells.

We recorded spikes and LFPs in the subiculum of freely running rats to determine (1) whether the distributed spatial representation in the subiculum might be advantageous for conveying information about spatial location and environmental context and (2) whether subicular neuronal activity exhibits a phase code as well as a rate code. We discovered that neurons in the subiculum exhibit robust theta phase precession and that the distributed code found in the subiculum is well suited to transmit information about both the animal's current spatial location and environmental context using many fewer neurons than the upstream CA1 region.

## Materials and Methods

##### Extracellular recordings in behaving rats.

The data presented here come from five adult (400–600 g) male Long–Evans rats. All procedures were approved by the Institutional Animal Care and Use Committee at the University of California, San Francisco. Rats were handled and trained to run back and forth on a U-shaped elevated running track (6 cm wide, 350 cm long) for drops of saccharin-sweetened milk. During training and later recording, food intake was restricted to maintain each rat at 85–90% of its baseline free-feeding body mass. The daily training schedule consisted of three 15 min task sessions on the U track, with intervening 15 min rest sessions in a nest enclosure. Rats were trained for at least 5 d before undergoing surgery for electrode implantation. After recovery from surgery, rats were retrained for at least five additional days to ensure that the load on the head and the cable tether did not impair running.

Rats were implanted with microdrive arrays that bilaterally targeted the intermediate subiculum (Table 1). In each hemisphere, six tetrodes and one single-wire electrode made parallel penetrations into the brain tissue, with ∼0.3 mm spacing between adjacent electrodes. The electrodes were constructed from polyimide-insulated 12.5 μm nichrome wire, and the recording tips were plated with gold to achieve impedance in the range of 200–250 kΩ at 1 kHz. The single-wire electrode in each hemisphere was used to record a local reference signal in the white matter overlying the subiculum. Spike waveforms (600–6000 Hz) were recorded relative to the local reference. Local field potentials (1–400 Hz) were recorded relative to a stainless steel screw in the skull over the midline cerebellum which served as a common ground.

Recordings began on the last day of training on the familiar U-track (“day 0”). On the next day of recording (“day 1”) and on all subsequent days, the second 15 min run session took place in a geometrically identical U track with the same surface color and texture, located on the other side of the room which the rat had never before visited. The two sides of the room were separated by a tall barrier and contained different visual landmarks. We refer to the training U track as “environment 1” and the second, novel track as “environment 2.” The two environments were separated by a wall so that they were not visible from each other. For one rat, environment 1 and environment 2 were oriented anti-parallel to each other; for the other four rats, the two environments were oriented parallel to each other, so that the rat faced the same wall when approaching the food wells. On each day after the conclusion of recordings, tetrodes were adjusted to optimize recording yield for the following day. Recordings continued for as many days as well isolated single units could be recorded.

Details of data acquisition and off-line processing have been published previously (Karlsson and Frank, 2008). Only well isolated putative single units that had an absolute refractory period of at least 2 ms (1 ms for fast-spiking neurons with high firing rates) were included in the final analyses. For units that were cleanly isolated in only a subset of the recording sessions (e.g., as a result of slow tetrode drift), we included spike data from only those sessions, provided that the waveform cluster was identifiable across all sessions without the possibility of confusion with other units. We did not attempt to match putative single units across days, so it is possible that the some of the same neurons may have been recorded in multiple days. We included data from “silent cells” that fired few spikes during the task sessions (Thompson and Best, 1989), which we identified by their stable spike waveforms across intervening rest sessions.

We recorded video of the rats' running behavior in the two environments with an overhead video camera. To ensure visibility, infrared light-emitting diodes were fixed on the animals' heads. In off-line video processing, we determined the pixel coordinates of the light-emitting diodes in every video frame with a semiautomated object-tracking algorithm. We then transformed these pixel coordinates to physical distances along the length of the environment, using the known geometry of the environment to calibrate image scale and camera distortion. Next, we smoothed the linearized position data in two steps. In the first step, we used robust locally weighted regression to suppress outliers and pixel jitter (Cleveland and Loader, 1996). In the second step, we fitted a smoothing spline with a data-adaptive roughness penalty (de Boor, 2001). The roughness penalty of the smoothing spline was based on the third derivative of the fit, which guaranteed the smoothness of the first derivative (running speed), and it was weighted in proportion to the local mean squared relative acceleration (acceleration divided by speed) at each time point. The adaptive weighting of the roughness penalty suppressed spurious high-speed transients due to head-bobbing when the rat was standing in place, while preserving veridical accelerations and decelerations that occurred when the rat was running. We estimated running speed by analytical differentiation of the fitted smoothing spline.

##### Histological reconstruction of recording sites.

At the end of the experiment, each rat was anesthetized with isoflurane, and direct current (10 μA for 3 s) was passed through the electrodes to make small electrolytic lesions at the recording tips. One day later, the rat was lethally overdosed with pentobarbital and transcardially perfused with isotonic sucrose followed by formaldehyde fixative. After perfusion, the brain was further soaked in fixative for 24 h with the electrodes remaining *in situ*; this postfixation step made the electrode penetration tracks stand out clearly in the fixed tissue. The fixed brain was cryoprotected by infiltration in 20% glycerol/2% DMSO, cut in 50 μm coronal sections on a cryostat microtome, and stained with cresyl violet. Cytoarchitectonic boundaries of area CA1, subiculum, and presubiculum were resolved according to established criteria (Witter and Amaral, 2004). We subdivided the subiculum into longitudinal strips corresponding to proximal (closer to CA1), middle, and distal (closer to presubiculum) thirds along the transverse dimension. We also distinguished, separate from the subiculum proper, a CA1/subiculum transition zone and a subiculum/presubiculum transition zone where the principal cell layers overlapped.

##### LFP theta oscillations.

We recorded all LFPs relative to ground. To estimate the instantaneous phase of theta oscillations, we bandpass filtered the LFP and decomposed the filtered signal into amplitude and phase components using the Hilbert transform (Harris et al., 2003; Siapas et al., 2005). We used a phase-preserving acausal filter with a 5–10 Hz bandpass; this bandpass was chosen to accommodate fluctuations in the instantaneous frequency of theta (Montgomery et al., 2008). We defined 0° phase to be the positive peak of the theta oscillation.

Instantaneous theta phase could not of course be reliably estimated whenever theta oscillations became indistinct. To ensure that we always worked with a well defined theta phase, we restricted our phase-dependent analyses to times when there was a clear spectral peak in the 5–10 Hz frequency band of the LFP power spectrum. We determined when this peak was present by first prewhitening the LFP to compensate for the 1/*f* shape of the power spectrum, so that our estimate of power in the 5–10 Hz frequency band would not be overwhelmed by leakage of background power from lower frequencies (Percival and Walden, 1993). After prewhitening the LFP, we estimated the power spectrum in overlapping 2 s time windows using the multitaper method (*K* = 3 Slepian tapers, *W* = 1 Hz; Mitra and Pesaran, 1999). The resulting whitened power spectrum was approximately flat over the frequency intervals 1–4 Hz and 18–400 Hz. We defined a “theta power ratio” to identify time windows during which theta power was elevated against this flat background. The theta power ratio was computed as the integrated power of the whitened spectrum over 5–10 Hz (theta), divided by the integrated power over 10–25 Hz (supratheta). We specifically avoided using the delta frequency band in our calculation of the theta power ratio, because low-frequency (<4 Hz) movement artifacts occurred when rats were running. Contiguous overlapping time windows in which the theta power ratio exceeded 5 dB were determined to be times when theta oscillations were clearly present in the LFP. We checked that our results were robust to reasonable changes in the choice of threshold, as well as choice of parameters in the multitaper method.

We used circular statistics to measure the relationship between single-unit spiking and the phase of ongoing theta oscillations in the LFP (Siapas et al., 2005). We determined an instantaneous phase value for every spike by interpolating the theta phase component of a reference LFP signal at the spike times. In practice, we found that the choice of reference LFP signal did not qualitatively alter any of our results, because LFPs were highly coherent in the 5–10 Hz frequency band across all recording sites. However, for the sake of definition, we referenced each neuron's spikes to the LFP signal recorded in the ipsilateral subiculum that showed the largest-amplitude theta oscillations.

We tested the unimodality of each neuron's spike phase distribution with the Rayleigh test of circular uniformity. Because theta oscillations in the LFP have a characteristic sawtooth asymmetry, the LFP theta phases that we estimated using the Hilbert transform were slightly nonuniformly distributed. To adjust for this baseline nonuniformity, we transformed the spike phases to circular ranks in the empirical distribution of LFP theta phases, and then performed the Rayleigh test on these circular ranks (Siapas et al., 2005). We also fit a von Mises distribution to the spike phases via maximum likelihood. The estimated mean parameter of the von Mises fit was taken as the neuron's preferred phase of firing, and the estimated concentration parameter was taken as a measure of phase locking to theta oscillations in the LFP. A larger value of the von Mises concentration parameter indicates a smaller circular variance, i.e., greater selectivity for the preferred phase of theta.

##### Burst index.

Previous work has suggested the presence of different subicular cell types based on burstiness (Sharp and Green, 1994), so we calculated a burst index to describe subicular neurons. A conventional measure of burstiness is the proportion of interspike intervals that are shorter than some threshold value, typically 6–10 ms (Frank et al., 2001; Harris et al., 2001; Anderson and O'Mara, 2003). This conventional measure is reasonable for comparing neurons that have low mean firing rates, such as principal neurons in the hippocampus, but it gives misleadingly inflated values for neurons with high mean firing rates. We observed a wide range of mean firing rates among neurons in the subiculum, so we needed a measure of burstiness that was not confounded by firing-rate differences. We therefore defined a new “burst index” that was based on a neuron's spike autocorrelogram, rather than the distribution of interspike intervals. This burst index was computed as the integrated power of the spike autocorrelogram over 1–6 ms, divided by the integrated power over 1–20 ms. Larger values indicate more burstiness. Our burst index is similar to the “first moment” of the (truncated) autocorrelogram (Csicsvari et al., 1999), but unlike the first moment, it is selectively weighted to measure bursts.

##### Firing-rate maps.

To characterize rate and phase coding by neurons in the subiculum and area CA1, we estimated each neuron's instantaneous firing rate as a function of the rat's linearized position in the environment and the phase of theta oscillations in the LFP (Mehta et al., 2002; Maurer et al., 2006). Separate position-phase firing-rate maps were estimated for the right-bound and left-bound directions of travel in each environment. We required data from a minimum of five passes in a given environment/direction to estimate a corresponding firing-rate map. We defined a pass to be a single traversal of the environment from one food well to the other food well. We excluded passes in which the rat backtracked or failed to attain a mean running speed of at least 10 cm/s. We also excluded time intervals when the theta power ratio fell below the minimum 5 dB criterion, because theta phase was ill-defined in the absence of clear theta oscillations in the LFP; furthermore, if these intervals of low theta power comprised >10% of the duration of a pass, we excluded the entire pass.

To accurately estimate position-phase firing-rate maps, we needed a statistically efficient estimator that could cope well with noisy and undersampled data. Because the rats behaved freely in our experiments, the number of passes in any given environment/direction was modest and there was pass-by-pass variability in the sampling of linearized positions and LFP theta phases; consequently, the position-phase space (that is, the set of all possible combinations of position and phase) was unevenly and sparsely sampled. To overcome this sampling variability and the intrinsic spiking variability of neurons, we used local polynomial kernel regression (Fan et al., 1995) to smooth the data, under the natural assumption that a neuron's instantaneous firing rate should be a smooth function of linearized position and theta phase. This smoothness assumption reflects our prior knowledge about the shapes and sizes of the neurons' firing fields (Barnes et al., 1990; Sharp and Green, 1994). Details of the statistical estimation process and how it makes use of finite amounts of data without overfitting rate maps, are given in a section at the end of these methods.

##### Spatial activity fraction.

We quantified the sparseness of spatial coding in the hippocampus and subiculum by the spatial activity fraction, which can be interpreted as the proportion of the environment over which a neuron fires (Rolls et al., 1998; Battaglia et al., 2004; see also Olypher et al., 2003; Ego-Stengel and Wilson, 2007; Wilent and Nitz, 2007 for information theoretic approaches). This measure is identical to standard measures of spatial sparseness (Skaggs et al., 1996) adapted to be appropriate for the continuous theta phase variable.

Given a position-phase firing-rate map *f*(*x*,φ), where *x* is linearized position on the track and φ is the instantaneous phase of the LFP theta oscillation, the spatial activity fraction was computed as follows:
where
is the spatial firing-rate profile averaged over all theta phases. The spatial activity fraction ranges between zero and one. A value close to zero indicates that the neuron fires in only a small region of the environment and is silent everywhere else, whereas a value close to unity indicates that the neuron fires uniformly at all locations throughout the environment.

##### Comparison of firing-rate maps across environments.

We used two complementary measures to quantify the similarity between the firing-rate maps of a given neuron in environment 1 versus environment 2. The first measure, which we call cosine similarity, was computed as follows:
Here, *f*_{1} and *f*_{2} are the neuron's position-phase firing-rate maps in the two environments. The cosine similarity is the mean, taken over all theta phases, of the normalized dot product between spatial firing-rate profiles conditioned on theta phase. Because firing rates are strictly positive, this measure can assume values between zero and one; a value of unity indicates that the two firing-rate maps are identical up to a proportionality constant. The second measure, which we call normalized overlap, was computed as follows:
where
is the firing-rate map normalized at each theta phase. The normalized overlap is a generalization of place-field overlap (Battaglia et al., 2004) to account for LFP theta phase. This measure also ranges between zero and one. For the purposes of computing cosine similarity and normalized overlap, we concatenated firing-rate maps for the right-bound and left-bound directions of travel in each environment.

##### Mutual information.

Using the experimentally derived firing-rate maps of the neurons that we recorded, we constructed model populations of neurons and simulated their ensemble spiking responses. From these simulations, we estimated how much spatial information is carried by neurons in the subiculum and in area CA1, while making as few assumptions as possible about how this information might be decoded by downstream circuits. Our combined use of experimental data, modeling, and information theory is similar to the approach used by Osborne et al. (2008). The inferences that we derived are valid to the extent that the firing-rate maps of the neurons in our dataset were representative of the true coding statistics in the larger neuronal populations.

We had three reasons to rely on an experimentally derived model instead of estimating information-theoretic quantities directly from the experimental data. First, estimates of mutual information from limited samples are biased, and corrections for this sampling bias require considerable effort (Panzeri et al., 2007). By using a model, we could perform exact calculations that were equivalent to simulating infinite amounts of data. Second, estimates of mutual information depend on both the neural response and the stimulus distribution (in this case, the occupancy distribution over spatial locations and environments). In our model, we fixed the stimulus distribution to be exactly the same for all neurons, so that we could fairly compare mutual information between different neurons and regions without having to worry about confounding differences in spatial behavior. Third, the model allowed us to assess combinatorial coding by diverse ensembles of neurons. In our experiments, we were never able to simultaneously record more than a few well isolated single units, so it was necessary to model pseudo-simultaneous ensemble responses from data recorded over multiple days in multiple animals.

While our model-based information-theoretic approach has a number of advantages, correct interpretation of the results requires an understanding of its assumptions and limitations. In our model, we assumed that each neuron's spike train could be described as a conditionally independent inhomogeneous Poisson process, given the rat's current position and the instantaneous phase of theta oscillations. In other words, we did not model bursting or other history dependence in the spike trains, nor did we model “excess” correlations between neurons that might be present in the data. Thus, the number of bits of information that we calculated cannot be interpreted literally. That said, our approach does provide a principled way to compare the information content of different firing-rate maps, and that is how the results should be understood.

We constructed populations of model neurons for each of the following anatomical subregions: distal subiculum, middle subiculum, proximal subiculum, and distal area CA1 (including the CA1/subiculum transition zone). The responses of each model neuron were completely described by position-phase firing-rate maps of a corresponding real neuron in our dataset. For simplicity, we ignored directionality in our model; that is, the virtual rat ran in the same direction on every pass. This unidirectionality allowed us to duplicate each real neuron in our model, by mirror-reflecting the estimated firing-rate map in one of the two directions and treating this as an additional model neuron. We always coupled the same real neuron's firing-rate maps in environment 1 and environment 2, so that the corresponding model neurons properly reflected the observed degree of similarity/remapping between the two environments.

From these model populations, we drew ensembles of *K* neurons at a time, for *K* = 1,…, 7. Because it was computationally intractable to simulate every possible combination of neurons for *K* > 2, we randomly sampled at most 500 unique ensembles from each population for each value of *K*. We simulated the responses of these *K* neurons, indexed by *k* = 1,…, *K*, using the experimentally derived position-phase firing-rate maps, *f _{k}*(

*x*,φ,

*e*).

*x*is the linearized position of a (virtual) rat, φ is the phase of theta oscillations in the LFP, and

*e*∈{1,2} is a variable that indicates whether the rat is running in environment 1 or environment 2. We simulated realistic distributions of

*x*and φ by randomly selecting 100 passes from our dataset in which the animal's mean running speed exceeded 10 cm/s and the theta power ratio was above the minimum 5 dB criterion for at least 90% of the pass duration. These passes were identically replicated in both environment 1 and environment 2 to simplify our calculations and to enable us to estimate how much information was conveyed about the identity of the environment.

We divided each pass into short (but not infinitesimally small) nonoverlapping time bins of duration Δ*t*. We assumed that information was conveyed in the patterns of simultaneous spiking and silence among the ensemble of neurons within each time bin. For simplicity, each neuron's response in a single time bin was denoted as either “0” (silence) or “1” (one or more spikes). With an ensemble of *K* neurons, 2* ^{K}* possible binary ensemble responses (“words”) were possible. Patterns of spiking across time were not considered. We set the bin size Δ

*t*to be 20 ms, because this is consistent with the known temporal resolution of spiking in the hippocampus (Harris et al., 2003; Lisman, 2005). However, we obtained qualitatively similar results with 10 and 30 ms bins.

The neurons in our model fired Poisson spike trains and were conditionally independent given their respective firing-rate maps. For a single neuron, the probability of spiking/silence in a given time bin was calculated as follows:
where *n _{k}* is the spike response of the neuron and

*f*(

_{k}*x*,φ,

*e*) is the neuron's firing-rate map. The time-varying continuous variables

*x*(

*t*) and φ(

*t*) were taken directly from experimental data, and the environment

*e*was imposed by simulation. Assuming conditional independence among

*K*neurons in an ensemble, the probability of the ensemble word

*n⃗*≡ [

*n*

_{1}, …,

*n*] was then simply Because the virtual rat performed identical passes in the two environments, we had the following simple marginal probabilities: Note that

_{k}*p*(

*x*) was not uniform because rats sped up and slowed down in a stereotyped manner over the environment across the multiple passes. The mutual information between the ensemble responses

*n⃗*and the spatial variables (

*x*,

*e*) is defined to be the difference between the total entropy of spike responses and the noise entropy of the spike responses: The first term is the total entropy of the response, and the second term is the noise entropy, wherein we average the conditional distribution of response over all linearized positions and both environments. Notice that although the LFP theta phase φ contributes to the spiking probabilities, it does not appear directly in the mutual information. φ is effectively a nuisance variable, which we included to realistically model theta modulation and phase precession. The simulated passes by the virtual rat contained enough random variation that φ was uncorrelated with

*x*.

We also estimated the mutual information between ensemble responses *s⃗* and the environment *e*, conditioned on the current linearized position *x*:
This conditional mutual information quantifies how much information the neurons convey to distinguish between locations in environment 1 versus environment 2. We used this measure to quantify the remapping of spatial representations between the two environments.

We then performed additional simulations to determine whether theta phase precession was essential for conveying spatial information. We modified the position-phase firing-rate maps of the model neurons to remove theta phase precession while preserving spatial modulations of firing rate. For each experimentally derived firing-rate map *f*_{experiment}, we computed a rank-1 approximation *f*_{approx}, using the singular value decomposition (Linden et al., 2003). The rank-1 approximation *f*_{approx} is separable, meaning that it can be factored as the product of a position-only component and a phase-only component:
Thus, by construction, a model neuron whose firing-rate map is given by *f*_{approx} will not exhibit theta phase precession. *f*_{approx} is an optimal approximation to the experimentally derived firing-rate map *f*_{experiment} in the following sense:
That is, the spatial profile of the firing-rate map and the phase profile match the original as closely as possible, subject to the separability condition.

Our estimate of mutual information is related, but not identical, to the “spatial information rate” that has been widely used to characterize hippocampal spatial representations (Skaggs et al., 1993; see also Olypher et al., 2003; Ego-Stengel and Wilson, 2007 for alternative measures of spatial information). The spatial information rate is also a model-based estimate of mutual information. However, it is limited to a single neuron's spiking response, and it is computed from the neuron's firing rate as a function of position alone without considering theta modulation. In contrast, we modeled the dependence of firing rate on LFP theta oscillations to simulate realistic spiking responses, and we enumerated ensemble patterns of spiking and silence to account for combinatorial coding synergies, which cannot be captured by simply summing the spatial information rates of individual neurons.

##### Phase precession in unitary place fields.

Theta phase precession is conventionally quantified by computing the correlation between position and the preferred theta phase of spiking within a neuron's place field (Mehta et al., 2002; Huxter et al., 2003). If a neuron has multiple place fields in an environment, these unitary place fields are segmented and a separate correlation is computed for each field. Previous investigators simply segmented place fields by the local minima and maxima of spatial firing rate as a function of linearized position alone, without considering theta phase (Hafting et al., 2008; Mizuseki et al., 2009). However, segmentation in the position dimension alone is prone to error, because two adjacent place fields can be partially overlapping in position but clearly separated in theta phase (Maurer et al., 2006). We therefore developed a new approach using morphological grayscale algorithms (Vincent, 1993) to segment unitary place fields in position-phase firing-rate maps. Our procedure is fully automated, objective, and involves only two parameters: a minimum peak firing rate (*f*_{min}), and a parameter (ε) that specifies the contrast of the segmented place field against the “background” bumpiness of the firing-rate map.

Given a position-phase firing-rate map *f*(*x*,φ), we first identified all peaks (local maxima) that exceeded *f*_{min}. We set *f*_{min} to be 10 spikes/s, although qualitatively similar results were obtained with a more stringent threshold of 15 spikes/s or a more permissive threshold of 5 spikes/s. All identified peaks were checked sequentially in decreasing order of firing rate. Given a peak with firing rate *f*_{peak} > *f*_{min}, we applied the H-maxima transform to suppress all local maxima in the firing-rate map whose height was smaller than ε · *f*_{peak}. This smoothing step was necessary to make the procedure robust to noisy bumps in the firing-rate map. We set ε to be 0.3 (30%). Next, we used a flood-fill algorithm to trace the ε · *f*_{peak} contour that surrounded the base of the firing-rate peak, and we also used a watershed algorithm to identify watershed lines between neighboring convex domes in the firing-rate map. The traced contour was accepted as a valid place-field segmentation only if it did not intersect the contours of any other place fields that had already been segmented, did not intersect any watershed lines, did not intersect itself, and did not extend to the minimum or maximum linearized position of the firing-rate map. These criteria distinguished well isolated, unimodal, completely segmented fields from multiple overlapping fields or truncated fields. The contrast parameter ε tunes the algorithm's tolerance for departures from strict unimodality. Setting ε to a more stringent value of 0.2 made the algorithm too sensitive to small bumps in the firing-rate map and excluded many place fields that appeared to be effectively unimodal by visual inspection, such as fields with closely spaced twin peaks or small side lobes on their flanks. It should be noted that our morphological segmentation procedure is completely symmetric and invariant with respect to periodic theta phase; we did not impose any prior bounds on the peak phase of the segmented place field or its orientation in position-phase space.

After segmenting the unitary place fields, we then identified their principal axes. Place fields that exhibit phase precession have a characteristic oblique orientation in position-phase space. To capture this orientation, we fit a two-dimensional Gaussian function to the firing-rate surface within the segmented baseline contour of each place field. We did not impose any prior bounds on this Gaussian fit; “backwards” place fields with positive-valued phase/position slopes could be fit just as well as place fields with negative slopes. From the covariance matrix of the Gaussian fit, we computed the slope (in degrees/cm) and the Pearson correlation coefficient between linearized position and LFP theta phase. We repeated the Gaussian fitting procedure over all jackknife estimates of the firing-rate map to estimate a Studentized confidence interval on the Fisher-transformed correlation coefficient. We excluded place fields whose correlation coefficients were not significantly different from zero, because the phase/position slope of a field is ill-defined in the absence of correlation.

##### Spike phase spectra.

The preceding analysis of theta phase precession on a per-field basis was limited to place fields that could be automatically segmented. To confirm the generality of our results, we performed an alternative analysis that did not rely on segmentation of unitary place fields. We measured theta phase precession by estimating the peak of the spike phase spectrum (Mizuseki et al., 2009; Geisler et al., 2010). The spike phase spectrum is the power spectral density of a spike train after it has been transformed from a sequence of spike times to a sequence of instantaneous spike phases relative to the ongoing LFP theta oscillation. The frequency domain of the spike phase spectrum has units of relative frequency, where the instantaneous frequency of theta oscillations in the LFP is defined to be 1 cycle^{−1}. A neuron that exhibits theta phase precession has a peak in its spike phase spectrum at a frequency greater than unity. A higher peak frequency indicates a faster rate of phase precession.

We introduced three methodological improvements to the estimation and interpretation of the spike phase spectrum. First, we estimated the spike phase spectrum using only data from valid running passes when LFP theta power was high, deliberately excluding times when rats were stopped at the food wells. Second, we used the derivative of the spike phase spectrum to accurately resolve the frequencies of peaks in the spectrum. Third, we performed shuffle tests to identify statistically significant spectral peaks. We describe these methods in further detail in a section below.

##### Estimating position-phase firing rate maps.

To motivate the use of local polynomial kernel regression and to introduce necessary notation, we first consider a simple estimator that is known in statistics as the Nadaraya-Watson kernel smoother (Hastie et al., 2009). The Nadaraya-Watson kernel smoother has been used in previous studies to estimate the firing-rate maps of hippocampal place cells (Harris et al., 2001, 2003). Conceptually, it is straightforward: simply divide a smoothed spike density map by a smoothed position-phase occupancy map to obtain a smoothed firing-rate map. Formally, the computation proceeds as follows: We discretize the data in small time steps (Δ*t* = 2 ms) that are indexed by *t* = 1,…, *T*. In this discrete representation, the number of spikes fired by the neuron in each time step (typically zero or one) is given by the time series *N*_{1},…, *N _{T}*; the linearized position of the rat is given by

*X*

_{1},…,

*X*; and the phase of the LFP theta oscillation is given by Φ

_{T}_{1},…, Φ

*. Then the Nadaraya-Watson estimator of the neuron's instantaneous firing rate*

_{T}*f*at a given linearized position

*x*and LFP theta phase φ is given by the following formula: where

*K*is a kernel function supported in [−1, +1]⊗[−1, +1]. The phase difference (Φ − φ) in the argument of the kernel function must be taken in a way that respects the periodicity of phase.

The parameters *h _{x}* and

*h*

_{φ}specify the smoothing bandwidth of the Nadaraya-Watson kernel smoother. The kernel-weighted contributions of nearby data points that fall within the position bandwidth [

*x*−

*h*,

_{x}*x*+

*h*] and the phase bandwidth [φ −

_{x}*h*

_{φ},φ +

*h*

_{φ}] are averaged together to produce the estimate

*f̂*(

_{N-W}*x*,φ). If

*h*and

_{x}*h*

_{φ}are too small, then the Nadaraya-Watson kernel smoother will exhibit unacceptably high variance, because few data points will be included in the smoothing bandwidth around each point in position-phase space. On the other hand, if

*h*and

_{x}*h*

_{φ}are too large, then the Nadaraya-Watson kernel smoother will systematically underestimate the heights (depths) of peaks (valleys) in the firing-rate map, because the local structure of the data within the smoothing bandwidth will be averaged away. In other words, the Nadaraya-Watson kernel smoother controls sampling variance at the cost of increasing bias in regions of position-phase space where the true firing-rate map

*f*is curved. Unfortunately, these regions of curvature—i.e., place fields—are precisely the regions where we want to accurately quantify a neuron's response.

We can achieve a better tradeoff between variance and curvature bias by adding parameters to explicitly model the local curvature of the firing-rate map. This is the basic principle of local polynomial kernel regression (Fan et al., 1995). Assuming that *f* is everywhere smooth, we can locally approximate log *f* with a low-degree polynomial. (In fact, the Nadaraya-Watson kernel smoother happens to be a special case of local polynomial kernel regression with a 0-degree polynomial, i.e., a constant function.) We fit the logarithm of the firing rate with a polynomial to guarantee that the estimated firing rate is always positive. Furthermore, the logarithm is the canonical or natural link function for Poisson-distributed count data, with attractive mathematical properties for efficient maximum-likelihood estimation (McCullagh and Nelder, 1989).

For example, the quadratic Taylor approximation to log *f* in the neighborhood of the fitting point (*x*,φ) is the following:
where
is a vector of coefficients. β* _{x}*(

*x*,φ) and β

_{φ}(

*x*,φ) approximate the first-order partial derivatives of log

*f*at (

*x*,φ); likewise, β

*(*

_{xx}*x*,φ), β

_{φφ}(

*x*,φ), and β

_{x}_{φ}(

*x*,φ) approximate the second-order partial derivatives. We can use nearby data points that are within a smoothing bandwidth around the fitting point to estimate these derivatives. Just as the Nadaraya-Watson kernel smoother averages the kernel-weighted contributions of data points around each fitting point, local polynomial kernel regression maximizes the kernel-weighted log-likelihood of the data (Fan et al., 1995). Assuming that the spike counts are Poisson-distributed, maximization of the kernel-weighted log-likelihood takes the following form: where the kernel function

*K*and the bandwidth parameters

*h*and

_{x}*h*

_{φ}are the same as in the Nadaraya-Watson estimator. This log-likelihood maximization problem is equivalent to fitting a generalized linear model with weighted observations, for which there are well developed numerical recipes (McCullagh and Nelder, 1989). Exponentiating the constant coefficient of the fit gives the “local quadratic” estimator of the firing-rate map: Although the local quadratic estimator exhibits less curvature bias than the Nadaraya-Watson kernel smoother, it is susceptible to overfitting in regions of position-phase space where the firing rate is very close to zero. The Nadaraya-Watson kernel smoother, in contrast, is immune to overfitting and performs well in regions where the firing rate is flat and close to zero. Combining the two estimators yields a robust shrinkage estimator that adapts to the local curvature of the data without overfitting: The shrinkage parameter α(

*x*,φ) takes values between zero and one, and can be selected in a data-adaptive manner by cross-validation, as described below.

We estimated firing-rate maps on a uniform grid of points that were spaced 2 cm apart in the linearized position dimension and 6° apart in the theta phase dimension. This grid spacing was fine enough to support smooth interpolation between grid points. We excluded linearized positions within 20 cm of either food well, because theta oscillations in the LFP often were indistinct when rats stopped to feed, and also because linearized position was not well defined during turning behavior. At each grid point, we computed both the Nadaraya-Watson kernel smoother and the local quadratic estimator. We used a tensor-product Epanechnikov kernel:
which downweights data points at the outlying extremes of the smoothing bandwidth. We chose bandwidth parameters *h _{x}* = 20 cm and

*h*

_{φ}= 90° to match the scale of features in the firing-rate maps. To confirm that this smoothing bandwidth was reasonable, we also tested smaller and larger values (

*h*= 15 cm, 30 cm and

_{x}*h*

_{φ}= 60°, 120°) on data from a randomly selected subsample of neurons. As assessed by visual inspection of the firing-rate maps and by cross-validation scores (see below), these smaller and larger values of the bandwidth parameters resulted in worse fits for a majority of the neurons in the random subsample.

To determine the optimal degree of shrinkage between the Nadaraya-Watson estimator and the local quadratic estimator, we partitioned the data into passes and performed cross-validation. Let *p* = 1,…, *P* be the pass number, and let *S _{p}* be the set of all time bins that belong to pass

*p*. We excluded each pass from the dataset and estimated firing-rate maps from the remaining

*P*− 1 passes; these leave-one-out estimates are denoted by

*f̂*

_{N-W}

^{(−p)}and

*f̂*

_{quadratic}

^{(−p)}. We then selected the shrinkage parameter α(

*x*,φ) that minimized the kernel-weighted cross-validation deviance: where is the expected spike count in time bin

*t*that was predicted from the cross-validation dataset

*t*∉

*S*. The deviance is a natural error measure for Poisson count data that is analogous to the sum of squared errors for normally distributed data. The cross-validation procedure also conveniently provided jackknife estimates of the firing-rate map, which we used to test the statistical significance of theta phase precession as described below.

_{p}##### Estimating spike-phase spectra.

For each neuron, separate spike phase spectra were estimated for the right-bound and left-bound directions of travel in each environment. We required data from a minimum of five valid passes in a given environment/direction to estimate a corresponding spike phase spectrum. To ensure that our estimates were not contaminated by spikes emitted at times when theta phase was ill-defined, we excluded passes in which the theta power ratio fell below the minimum 5 dB criterion for >10% of the pass duration. We also excluded time intervals at the start and end of each pass when the rat was within 20 cm of either food well, because rats were stopped and/or feeding at these reward locations, and theta phase precession is not observed in these behavioral states (Skaggs and McNaughton, 1996). We computed a multitaper estimate of the spike phase spectrum by taking a weighted sum of eigenspectra over passes and tapers (Jarvis and Mitra, 2001). Because different passes contained different numbers of theta cycles, we constructed a separate set of Slepian tapers for each pass to maintain a consistent smoothing bandwidth of ±0.05 cycle^{−1} (Walden et al., 1995).

The multitaper method reduces the variance of the spectral estimate by averaging the true spectrum over the smoothing bandwidth. This averaging flattens out local structure in the spectrum whose frequency scale is smaller than the smoothing bandwidth. As a result, peaks in the estimated spectrum have distorted shapes and are not sharply resolved. Simply selecting the numerical maximum of a multitaper spectral estimate is not a consistent estimator of the true peak. However, it is possible to resolve spectral peaks by examining the covariance structure of the multitaper eigenspectra. We used this method to estimate the peak frequency of the spike phase spectrum, which was the parameter of interest for measuring theta phase precession. A detailed explanation of the theory behind this method is given by Prieto et al. (2007).

Assuming that the spectral density varies smoothly with the bandwidth of the multitaper estimate, it can locally approximated by a second-degree Chebyshev polynomial expansion:
where *W* is the half-bandwidth, *T _{n}* are Chebyshev basis polynomials, and

*a*(

_{n}*f*) are coefficients. The coefficient of the linear term is approximately proportional to the first derivative of the spectral density at

*f*: Using this Chebyshev polynomial approximation, the local covariance of the eigenspectra within the bandwidth [

*f*−

*W*,

*f*+

*W*] can be expanded as where and

*V*

_{1},…,

*V*are the

_{K}*K*orthonormal Slepian tapers (in frequency space) that were used to compute the multitaper estimate.

*V*

_{1},…,

*V*are Hermitian functions, meaning that their real components are even functions and their imaginary components are odd functions.

_{K}*T*

_{0}and

*T*

_{2}are even functions, and

*T*

_{1}is an odd function. Because the domain of integration [−

*W*,+

*W*] is symmetric about zero,

*a*

_{1}(

*f*) can be estimated from the imaginary component of the complex covariance matrix

*C*(

_{ij}*f*).

Around each frequency *f*, we constructed a global covariance matrix with block-diagonal structure, in which each submatrix along the diagonal was the multitaper covariance matrix for a single pass. We estimated the first derivative of the spectral density, *a*_{1}(*f*) · *W*, by least-squares fit to the imaginary component of this global covariance matrix. We then found all local maxima in the spike phase spectrum (i.e., positive-to-negative zero-crossings of the first derivative) within the 0.6–1.4 cycle^{−1} frequency band. We deliberately chose this wide, inclusive frequency band to avoid selection bias.

After finding all candidate local maxima, we defined the peak of the spike phase spectrum to be the local maximum with the largest integrated spectral power within ±0.05 cycle^{−1} (the multitaper smoothing bandwidth). We wanted to be confident that this spectral peak reflected a significant oscillation in the neuron's spike train, rather than a spurious fluctuation due to finite sampling. To test statistical significance, we computed spike phase spectra for 500 surrogate spike trains in which each spike phase was independently jittered by a random offset drawn from a uniform distribution on the interval [−π, +π]. This jittering disrupted oscillatory structure on the theta timescale while preserving slower spatial modulations of firing rate. The test statistic was the ratio of integrated power within [*f*_{peak} −*W*, *f*_{peak} +*W*] divided by the total power within the 0.6–1.4 cycle^{−1} frequency band. We deemed an observed peak in the spike phase spectrum to be significant at the 0.05 level if the test statistic for the experimentally observed spike train was greater than the value attained by 95% or more of the phase-jittered spike trains.

## Results

We collected data from 5 rats while they ran on two geometrically identical running tracks that were located in separate partitions of a room. The rats experienced many training sessions on one of the tracks (environment 1) before recording, whereas the other track (environment 2) was relatively novel because it was introduced during recording. This experimental design allowed us to characterize neural representations of spatial context between different environments as well as spatial location within an environment. We did not record enough neurons every day to be able to analyze novelty-triggered dynamics as a function of the number of exposures to environment 2, so instead we combined data from all days of recording. Overall we found that the rats spent on average 49% of the total recording time (25 h across all 5 rats in the two tracks) not at a food well. Seventy percent of this time the rats ran toward a food well at velocities >5 cm/s, and during these running periods we observed a theta power ratio exceeding 5 dB 99% of the time.

### Diversity of firing patterns in the subiculum

Using tetrodes, we recorded spikes from 91 single units in the subiculum at intermediate locations along the septal-temporal axis and extending along the proximal-distal axis of the subiculum (Fig. 1*A*). Previous studies proposed various schemes to classify subicular neurons on the basis of firing rates, bursting patterns, and extracellular spike waveforms (Sharp and Green, 1994; Anderson and O'Mara, 2003). We therefore examined these parameters to identify salient differences in the firing properties of the neurons in our dataset (Fig. 1*B*,*C*).

We were able to distinguish separate fast-spiking and non-fast spiking classes of neurons, but we did not find any clear basis for classifying neurons according to burstiness. Three neurons in our dataset had narrow, symmetric spike waveforms and fired at high mean rates (>30 spikes/s) with no appreciable bursting tendency. Non-bursting, fast-spiking neurons in the subiculum have been previously identified as putative inhibitory interneurons (Greene and Totterdell, 1997; Menendez de la Prida et al., 2003). The other 88 neurons had broad, asymmetric spike waveforms and were quite diverse in their firing properties. There was no clear category boundary between “bursting” and “non-bursting” neurons, and burstiness did not significantly correlate with any features of the spike waveform that we measured. We classified these neurons as putative pyramidal cells, in agreement with previous studies that found substantial variability in burstiness among subicular neurons with indistinguishable pyramidal morphologies (Taube, 1993; Staff et al., 2000).

Neurons at different proximal-distal locations within the subiculum receive different inputs and send output projections to different targets (Witter, 2006). Sharp and Green (1994) reported that neurons in the proximal subiculum (closer to area CA1) have slightly lower mean firing rates and fire over a somewhat smaller proportion of the environment than neurons in the distal subiculum (closer to presubiculum). We therefore examined functional differentiations along the transverse axis of the subiculum.

We found a pronounced proximal-distal gradient in the distribution of firing rates among neurons in the subiculum. We grouped neurons by anatomical location in transverse thirds of the subiculum: proximal (closer to CA1), middle, or distal (closer to presubiculum). Plotting the mean firing rates of subicular neurons grouped by transverse location revealed a striking pattern (Fig. 1*C*). The median (interquartile range) firing rates of principal neurons in the proximal, middle, and distal groups were, respectively, 2.5 (1.4–3.7), 6.3 (4.1–11.4), and 10.3 (6.6–14.1) spikes/s. The effect of proximal-distal location on firing rate was statistically significant (*p* < 10^{−7}, Kruskal–Wallis test). We did not find a statistically significant effect of proximal-distal location on spike waveforms (Fig. 1*D*).

For reference, we compared the firing properties of neurons in the subiculum to those of simultaneously recorded neurons in area CA1 at nearby locations along the septal-temporal axis. In 3 of the 5 rats, some tetrodes penetrated the distal part of area CA1 and the transition zone where the principal cell layer of area CA1 superficially overlaps the principal cell layer of the subiculum. The neurons that we recorded at these locations seemed to constitute a homogeneous sample; none of our analyses revealed significant differences between neurons in the CA1/subiculum transition zone versus neurons in CA1 proper. Therefore, to improve the statistical power of comparisons with neurons in the subiculum proper, we assigned all of these neurons to a single “distal CA1” group, with the caveat that this may be an oversimplification of the neuronal population in the CA1/subiculum transition zone.

As expected (Barnes et al., 1990; Sharp and Green, 1994), we found that the median (interquartile range) firing rate of principal neurons in the distal CA1 group was 1.4 (1.1–1.7) spikes/s, which was significantly lower than the overall median firing rate of principal neurons in the subiculum (*p* < 10^{−9}, Kruskal–Wallis test). The firing rates of principal neurons in area CA1 cluster near zero; in the proximal subiculum, firing rates are slightly higher; and in the middle and distal subiculum, the firing-rate distributions are shifted even higher.

### Spatial representation in the subiculum

We also found significant proximal-distal differences in the spatial coding properties of neurons in the subiculum. We estimated firing-rate maps as a function of the animal's linearized position along the track and the instantaneous phase of theta oscillations in the LFP (see Materials and Methods). For each neuron, separate position-phase firing-rate maps were estimated in each environment and direction of running. Figure 2 shows examples of firing-rate maps for representative neurons in the proximal, middle, and distal subiculum. These firing-rate maps contain well defined regions of high firing rate in position-phase space, which we refer to as unitary place fields in line with previous terminology (Maurer et al., 2006). Note that these place fields have oblique orientations that are indicative of theta phase precession; this aspect of spatial coding will be addressed later in the paper.

We found that a majority of neurons in the subiculum have multiple, irregularly spaced unitary place fields, consistent with previous descriptions of “patchy” spatial firing (Barnes et al., 1990; Sharp and Green, 1994). Remarkably, the firing-rate maps seemed to vary from those in proximal subiculum (Fig. 2*A*) that closely resembled maps from CA1 neurons to those in the middle and distal parts of the subiculum (Fig. 2*B*,*C*) that resembled the multipeaked maps from entorhinal grid cells in linear environments (Hafting et al., 2008; Derdikman et al., 2009; Mizuseki et al., 2009). Almost all of the neurons that we recorded in the subiculum showed distinct patterns of activity in the right-bound versus left-bound directions of travel; this is consistent with the known directional selectivity of neurons in area CA1 (McNaughton et al., 1983) and also agrees with previous recording studies of subicular neurons during directional running on an 8-arm radial maze (Barnes et al., 1990).

To compare the spatial firing patterns of neurons in the proximal, middle, and distal subiculum, we computed spatial activity fractions based on the firing-rate maps (Battaglia et al., 2004). A spatial activity fraction close to zero indicates that the neuron is sparsely active in space and only fires in a single location in the environment, whereas a fraction close to one indicates that the neuron fires uniformly over all spatial location. Plotting the spatial activity fractions of subicular neurons grouped by transverse location revealed a proximal-distal gradient that was consistent with the firing-rate gradient (Fig. 3*A*). The median (interquartile range) spatial activity fraction of principal neurons in the proximal, middle, and distal groups were, respectively, 0.31 (0.19–0.43), 0.54 (0.32–0.68), and 0.69 (0.59–0.74). The effect of proximal-distal location on spatial activity fraction was statistically significant (*p* < 10^{−4}, Kruskal–Wallis test).

Comparison between regions confirmed previous findings that neurons in the subiculum tend to fire over a greater proportion of the environment than neurons in area CA1 (Fig. 3*B*). The median (interquartile range) spatial activity fraction of putative principal neurons in distal CA1 was 0.26 (0.21–0.37), which was significantly lower than the overall median spatial activity fraction in the subiculum (*p* < 10^{−10}, Kruskal–Wallis test). Thus, the different firing-rate distributions in area CA1 and the subiculum are mirrored by different patterns of spatial selectivity. We found essentially identical results by estimating position-rate firing maps from data using simple Gaussian smoothing (the log-quadratic estimator in the methods section). Thus this observed proximal-distal gradient in firing properties is robust to details in the statistical procedures for estimating firing rate maps.

These differences could not be explained by biases in sampling of data from individual animals. Table 2 shows the number of units recorded in each rat and in each brain region, as well as the number of recording sessions. Since we recorded different numbers of cells from each rat in each brain region, we examined whether the variation found in firing properties across brain regions could be explained by biases due to individual variation across rats. Figure 3*C* shows four plots, one for each brain region, in which for each unit in each recording session, the overall unit's firing rate is plotted against its spatial activity fraction. Each unit/session is color coded according to the identity of the individual rat. Together, these plots do not reveal any systematic variation in firing properties across individual rats that could account for variation in firing properties across regions.

Similarly, we found that unit isolation quality could not explain anatomical differences in firing rate or sparsity. We took each recording session for each neuron and computed two firing rate properties, overall mean firing rate and spatial activity fraction, and two measures of unit isolation quality, the L-ratio and isolation distance (Schmitzer-Torbert et al., 2005). We examined data on a per-session basis because unit isolation quality could in principle change between recording sessions as the animal is moved between environments. We then examined the relationships between these values across all neurons and recording sessions. We found that only 1.5% and 4.8% of the total variation in firing rates across neurons and sessions could be explained by variations in L-ratio and isolation distance respectively. Similarly, we found that only 0.05% and 0.03% of the total variation in spatial activity fractions across neurons and sessions could be explained by variations in L-ratio and isolation distance respectively. Moreover, when we subsampled units from all four regions until the groups had approximately identical mean values for both cluster quality measures, we still found significant differences in overall firing rate (*p* < 10^{−16}, Kruskal–Wallis ANOVA) and mean spatial activity (*p* < 10^{−15}, Kruskal–Wallis ANOVA) across the four regions, indicating again that our results are not due to differences in the quality of recordings.

### Spatial information content of subicular firing-rate maps

What are the functional consequences of these spatial activity differences on neural information processing? We applied information theory to address this question. Using experimentally derived position-phase firing-rate maps, we simulated realistic spiking responses and computed the mutual information between spatial variables and spiking in short (20 ms) time windows. This model-based information-theoretic analysis allowed us to extrapolate from our experimental data to infer differences in the amount of spatial information that is conveyed by neuronal populations in the proximal, middle, and distal subiculum and in area CA1.

First, we compared the spatial information content of the firing-rate maps of single neurons. We found that, on average, single neurons in the subiculum convey more spatial information than do single neurons in distal CA1. The median (interquartile range) of mutual information for neurons in the subiculum was 0.057 (0.036–0.080) bits, versus 0.021 (0.0081–0.056) bits for neurons in distal CA1. This difference was statistically significant (*p* < 0.001, Kruskal–Wallis test). Comparisons within the subiculum revealed a statistically significant effect of transverse location as well. The median (interquartile range) mutual information for single neurons in the proximal, middle, and distal portions of the subiculum was, respectively, 0.030 (0.011–0.052), 0.062 (0.040–0.088), and 0.062 (0.043–0.083) bits (*p* = 10^{−4}, Kruskal–Wallis ANOVA). *Post hoc* multiple comparison with the Tukey–Kramer method showed that the median information for the proximal group was significantly lower than the median information for the middle and distal groups, at the 0.05 statistical significance level.

These differences in spatial information content can be understood by considering the relationship between mean firing rate and information (Fig. 4*A*,*B*). Neurons with high mean firing rates tend to fire over a large proportion of the environment, with numerous up-and-down spatial modulations of intensity; conversely, neurons with low mean firing rates tend to fire sparsely and remain silent in most of the environment (Fig. 3*A*,*B*). A key insight of our information-theoretic analysis is that distributed, spatially modulated firing throughout the environment conveys more information than sparse, spatially localized firing. For example, a spatially selective neuron that fires in only a single place field will be silent at most times, and these silences convey little information about the animal's current spatial location because they ambiguously code for any place that is not within the neuron's place field. Stated in information-theoretic terms, this mostly silent neuron has a low response entropy, and because the response entropy is an upper bound on mutual information, the mutual information must also be low. In agreement with this theoretical argument, we found a significant positive correlation between mean firing rate and the information conveyed by single neurons in the subiculum and in area CA1 (Spearman *r* = 0.66, *p* < 10^{−9}). Neurons with very low firing rates (<1 spikes/s) conveyed little information, whereas neurons with higher firing rates conveyed more information, even though they were less spatially selective. Hippocampal area CA1 and the proximal subiculum contain a substantial fraction of neurons that are silent or only fire in a small region of the environment; these relatively uninformative neurons explain the lower information content of these regions. In contrast, the middle and distal subiculum contain neurons with higher mean firing rates and larger activity fractions, and most of these neurons exhibit informative spatial firing-rate modulations.

Next, we compared the spatial information content of multineuronal combinations of firing-rate maps. Given an ensemble of neurons with diverse firing-rate maps, the amount of information that is conveyed in a combinatorial ensemble code can exceed the sum of information that each single neuron conveys independently. We found that, for a given ensemble size, randomly selected ensembles of subicular neurons, particularly the middle and distal subiculum, carry substantially more spatial information than randomly selected ensembles of CA1 neurons (Fig. 4*C*). For all groups, the relationship between ensemble size and mutual information was very close to linear. We therefore quantified the gain in information per additional neuron by linear regression of the median mutual information against ensemble size. The slopes (bootstrap standard error) of the best fit lines for the proximal, middle, and distal portions of the subiculum were, respectively 0.037 (0.00094), 0.059 (0.0010), and 0.064 (0.00080) bits/neuron. The slope of the best fit line for area CA1 and the CA1/subiculum transition zone was 0.036 (0.0013) bits/neuron. We performed a one-way ANOVA on these slopes, using bootstrap estimates of the variances of the regression coefficients. This ANOVA revealed a highly significant effect of transverse anatomical location (*p* < 10^{−9}). Pairwise multiple comparisons with the Tukey–Kramer method revealed that distal area CA1 and the proximal subiculum had significantly smaller slopes than the middle subiculum and distal subiculum, at the 0.05 significance level. Thus, neuronal ensembles in the proximal subiculum convey as much spatial information as neuronal ensembles in area CA1, and neuronal ensembles in the middle and distal subiculum convey more spatial information.

Our information-theoretic analysis points to the functional consequences of sparse versus non-sparse spatial representations in area CA1 and in the subiculum. Sparse, selective spatial representation has an information cost, because neurons (like CA1 place cells) that are active in only a small fraction of the environment convey little information when they are silent. In contrast, the non-sparse, high firing-rate spatial representation in the middle and distal subiculum results in greater average information per neuron, because neurons that exhibit spatial firing-rate modulations over a large proportion of the environment consistently convey information for most of the time.

### Remapping of subicular spatial representations between two geometrically identical environments

Most place cells in the hippocampus exhibit distinct firing-rate maps in easily distinguished environments (Muller and Kubie, 1987; Best et al., 2001). Given the projections from area CA1 to the subiculum, we would expect these changes in spatial representation to propagate to the subiculum. However, published studies of the subiculum suggested that the spatial representation in the subiculum is much less sensitive to changes in environment than the upstream spatial representation in area CA1 (Sharp, 2006; Lever et al., 2009). To investigate this puzzle, we examined the firing-rate maps of neurons that were recorded in both environment 1 and environment 2.

We found clear evidence that the spatial representation in subiculum remaps across two geometrically identical environments that differed only in their familiarity to the rats and their allocentric locations. Of the subicular neurons that we recorded, 59 putative principal neurons had stable, well isolated, identifiable spike waveforms in both environment 1 and environment 2. We found many examples of neurons whose firing-rate maps were completely different between environment 1 and environment 2. Examples of between-environment remapping by representative neurons in the proximal, middle, and distal subiculum are shown in Figure 5. The diversity of remapping patterns is remarkable. Some neurons (like the one shown in Fig. 5*A*) had clear place fields in one environment and rarely spiked in the other environment. Other neurons maintained homotopic place fields across both environments, but with significant differences in firing rate, which has been described as “rate remapping” in the hippocampus (Leutgeb et al., 2005). Such rate remapping could account for previous claims about the invariance of subicular receptive fields (Sharp, 2006). Most commonly, we observed a complex form of remapping, where a subset of place fields were maintained across both environments while other portions of the firing-rate map changed. Given a neuron's firing-rate map in one environment, we could not find any parameter that reliably predicted how similar that neuron's firing-rate map would be in the other environment.

To quantify these observations over our entire dataset, we first examined differences in the firing rates of neurons across environment 1 versus environment 2. Overall, neurons in the subiculum maintained a similar level of activity in both environments (Fig. 6*A*). There was a nonsignificant trend for neurons in the distal subiculum to have lower firing rates in environment 2 than in environment 1 (*p* = 0.15, two-sided binomial test); this difference may be attributable to that fact that environment 2 was less familiar and therefore rats tended to run more slowly in environment 2. The overall correlation between firing rate in environment 1 and firing rate in environment 2 was strong among subicular neurons (Spearman *r* = 0.87, *p* < 10^{−9}). In contrast, principal neurons in distal area CA1 had much lower firing rates, and a substantial proportion had firing rates close to zero in one or both of the environments (Fig. 6*B*). Among neurons in the distal CA1 group, the correlation between firing rate in environment 1 and firing rate in environment 2 was weak and not statistically significant (Spearman *r* = 0.19, *p* = 0.49). With a firing rate of 1 spikes/s as a threshold for categorizing a neuron as active/inactive, 7 of 15 (47%) principal neurons in distal CA1 and the CA1/subiculum transition zone were active in one environment but not the other. Only 4 of 59 (7%) principal neurons in the subiculum exhibited such a contrast in activity between the two environments. Fast-spiking putative inhibitory interneurons in both the subiculum and in CA1 fired consistently at the same high firing rate in either environment (data not shown).

We also quantified the within-neuron pairwise similarity of position-phase firing-rate maps in environment 1 and environment 2 (Fig. 6*C–F*). The first measure, cosine similarity of receptive fields (DeAngelis et al., 1999), is equivalent to treating a given neuron's firing-rate maps in the two environments as vectors and taking their dot product. The second measure, normalized overlap (Battaglia et al., 2004; Singer et al., 2010) measures the overlapping areas of the spatial firing profiles. Both the cosine similarity and normalized overlap measures are sensitive to overall place field rates and sizes, so we performed within-neuron comparisons only, comparing the similarity of firing-rate maps in environment 1 versus environment 2 to the similarity across two exposures to environment 1. If the spatial representation in the subiculum were invariant to changes in environment, then we would expect, by chance alone, approximately half of the neurons to have higher similarity scores for [environment 1–environment 2] than for [environment 1–environment 1].

In fact, we found that a significant majority of neurons in the subiculum exhibited more similarity in their firing-rate maps between two sessions in the same environment than between different environments (cosine similarity: 39 of 39 neurons, *p* < 10^{−10}; normalized overlap: 37 of 39 neurons, *p* < 10^{−8}; binomial test). Performing the same analyses with principal neurons in distal CA1 showed that they also remapped between environment 1 and environment 2, as expected (cosine similarity: 10 of 11 neurons, *p* < 0.01; normalized overlap: 10 of 11 neurons, *p* < 0.01; binomial test). These results clearly demonstrate that the spatial representation in the subiculum, like the spatial representation in the hippocampus, is not invariant to changes of environmental context. This is in contrast to previous studies which emphasized the qualitative, though not exact, invariance of subicular representations across environments (Sharp, 1997, 2006).

We again used our model-based information-theoretic approach to rigorously compare how much information neurons convey about the animal's current environmental context. Because the two tracks were geometrically identical, we could compare the spiking probabilities of neurons at equivalent linearized positions across both environments. In terms of information theory, we estimated the conditional mutual information between ensemble spiking responses and the identity of the environment, given the animal's current linearized position on the track. We found that neurons in the subiculum, on average, are as informative about spatial context as neurons in area CA1. Figure 6, *G* and *H*, shows the gain in conditional mutual information per additional neuron in the ensemble. The slopes of the best fit lines (bootstrap standard error) for the proximal, middle, and distal portions of the subiculum were, respectively 0.012 (0.00027), 0.016 (0.00036), and 0.022 (0.00034) bits/neuron. For comparison, the gain in conditional mutual information per additional neuron in distal area CA1 and was 0.015 (0.00052) bits/neuron. We performed a one-way ANOVA on these slopes, using bootstrap estimates of the sample variances. This ANOVA revealed a significant effect of transverse anatomical location (*p* < 10^{−9}). Pairwise multiple comparisons with the Tukey–Kramer method revealed that the gain in conditional mutual information was significantly lower in the proximal subiculum versus the middle subiculum, and significantly higher in the distal subiculum, at the 0.05 statistical significance level. Thus, there is a proximal-distal gradient in the amount of environment-specific information that is conveyed in the subiculum. Neurons in the proximal subiculum convey slightly less information than those in area CA1, while neurons in the distal subiculum convey more information. Overall, the population of neurons in the subiculum is as informative about the identity of the current environment as the upstream population in area CA1.

### Theta modulation in the subiculum

The spatial information analyses revealed that the firing-rate maps of subicular neurons contain ample information about both the animal's location in an environment and the identity of that environment. We then asked whether subicular neurons also had the potential to participate in a temporal code associated with phase precession. We first examined the relationship between the firing of neurons in the subiculum and the phase of the theta oscillations in the LFP (Fig. 7).

Previous investigators found that the firing of neurons in the subiculum was significantly modulated by LFP theta phase (Anderson and O'Mara, 2003). We confirmed this result in our data and also discovered that distal subicular neurons were more strongly theta modulated than their proximal or distal subicular counterparts. Overall, we found significant theta power during 74% of the total recording time across all rats. Of the 91 neurons that we recorded in the subiculum, 89 had a statistically significant unimodal theta phase preference at the 0.05 significance level (Rayleigh test for nonuniformity). Categorizing neurons in the subiculum by transverse location revealed significant differences in the degree theta modulation (Fig. 7*B*,*C*). The population average firing rate in the distal subiculum was not only greater than that in the proximal subiculum, but also more strongly modulated by theta phase. Differences in theta phase modulation were also apparent when analyzed on a per-neuron basis. We quantified the phase locking of each neuron to LFP theta oscillations by fitting the distribution of spike phases with a von Mises distribution and calculating the concentration parameter of fit (Siapas et al., 2005). The concentration parameter assumes positive values, and larger values indicate stronger phase locking. Excluding fast-spiking neurons, median (interquartile range) spike phase concentration values in the proximal, middle, and distal portions of the subiculum were, respectively, 0.25 (0.16–0.34), 0.17 (0.11–0.29), and 0.34 (0.28–0.50) (*p* = 0.00014, Kruskal–Wallis ANOVA). Pairwise multiple comparisons with the Tukey–Kramer method revealed that neurons and proximal and middle subiculum had significantly weaker theta phase modulation than neurons in the distal subiculum at the 0.05 level, but the proximal and middle groups did not significantly differ from each other. To control for possible theta phase offsets between recording sites, we repeated these analyses with each neuron referenced to the local phase of theta recorded on the same tetrode on which the neuron's spikes were recorded, instead of a common reference electrode in each hemisphere. This did not qualitatively alter any of the results.

### Phase precession in the subiculum

The position-phase firing-rate maps in Figures 2 and 5 indicate that neurons in the subiculum exhibit theta phase precession. Many of the unitary place fields in these firing-rate maps are obliquely oriented in position-phase space, with a negative correlation between forward displacement along the track and the preferred theta phase of firing. In addition, phase precession is apparent over multiple successive place fields along the track. Similar multifield phase precession has been described for neurons in the hippocampus (Maurer et al., 2006) and in the EC (Hafting et al., 2008; Mizuseki et al., 2009).

To quantify theta phase precession over the entire dataset, we segmented individual place fields in the firing-rate maps and measured their sizes and slopes in position-phase space, using a fully automated, unbiased procedure (Fig. 8*A*). This procedure was designed to be equally sensitive to positive and negative slopes. Using this algorithm, we segmented 143 unitary place fields that had a statistically significant correlation between theta phase and linearized position. The correlation between phase and position was negative—that is, there was theta phase precession in the expected direction—in 142 of 143 of these unitary place fields (Fig. 8*B*).

We also measured the spatial extent of the place field along the linearized position dimension (field length) and found a significant correlation with phase precession slope, as has been reported for the place fields of CA1 neurons (Huxter et al., 2003; Dragoi and Buzsáki, 2006). Larger place fields had shallower phase precession rates (Spearman correlation = 0.65, *p* < 10^{−9}). Similar statistically significant correlations (all *p* < 0.001) were obtained when we separately analyzed neurons in the proximal, middle, and distal groups, and we also obtained a similar correlation with the place fields of neurons that we recorded in distal area CA1 (Fig. 8*D*). Place fields were preferentially centered around the trough of the LFP theta oscillation (Fig. 8*C*,*E*). These findings indicate that although neurons in the subiculum and CA1 differ in their spatial activity fractions and mean firing rates, the phenomenology of theta phase precession within unitary place fields is similar.

Although this place-field segmentation analysis demonstrated qualitative similarities in phase precession between the subiculum and CA1, it only included a small subset of the place field that could be cleanly segmented in an automated fashion. Many neurons in the subiculum had firing-rate maps that contained characteristic “sloped” regions indicative of phase precession, but unitary fields in these regions could not be segmented because they were overlapping (for example, see the middle part of the firing-rate map in Fig. 8*A*). To give a more complete accounting of phase precession which did not require segmentation of receptive fields and include all of the spike data, we computed the spike phase spectrum (Mizuseki et al., 2009). This method quantifies the relative frequency shift between the neuron's oscillatory spiking and the ongoing theta oscillation in the local field potential; phase advance due to theta phase precession is manifested as a peak frequency that is faster than the frequency of the LFP theta, which is defined as unity in the spike phase spectrum. Because most of the neurons fired differently in the two directions of running along the track and across the two environments, we estimated separate spectra for each environment and direction of travel. We used a multitaper covariance method to locate the peak of each spike phase spectrum, and we used a shuffle test to determine which spectra had significant peaks. This peak-finding procedure, illustrated in Figure 9*A*, was designed to be equally sensitive to peak frequencies above or below the theta frequency.

We found that neurons in the subiculum had peak frequencies which were consistently shifted above unity, indicating that their spiking oscillated faster than the local LFP theta oscillation. This was true for neurons in the proximal, middle, and distal thirds of the subiculum, and also held true regardless of the neuron's spatial activity fraction (Fig. 9*B*). There was, however, a weak but significant negative correlation between the spatial activity fraction and the peak frequency of the spike phase spectrum (Spearman correlation = −0.14, *p* < 0.03).

When we performed the same spike phase spectrum analysis with neurons in distal area CA1 and the CA1/subiculum transition zone, we found that principal neurons exhibited a range of shifted peak frequencies similar to those observed in the subiculum, even though neurons in distal area CA1 had significantly smaller spatial activity fractions and lower mean firing rates (Fig. 9*C*). Both in the subiculum and in CA1, fast-spiking putative interneurons, with spatial activity fractions close to unity, had peak frequencies that were closer to unity. In summary, principal neurons in the subiculum and in area CA1 oscillate at similar frequencies—offset from the frequency of theta oscillations—despite pronounced differences in firing rate and spatial activity fraction between these two regions and at different proximal-distal levels within the subiculum. This result implies that the mechanism of phase precession is not sensitive to differences in overall activity level or distributed spatial coding, and it also suggests that neurons in the subiculum are frequency-coupled to upstream cell assemblies in the hippocampus (Geisler et al., 2010).

Given that neurons in the subiculum exhibit theta phase precession, we examined whether phase precession is essential for the spatial information content of subicular firing-rate maps. We repeated our information-theoretic analysis with synthetically generated firing-rate maps in which theta phase precession (that is, the interactions between phase and position) had been removed while preserving each neuron's overall theta modulation and spatial firing profile. We found that replacing phase precession with simple theta modulation had little effect on the amount of spatial information conveyed in patterns of spiking and silence in short (10–40 ms) time windows. However, due to the computational limitations, we only measured information in single time bins and did not examine temporal correlations in spiking across multiple time bins. Thus it remains possible that additional information may be conveyed by sequential spiking patterns associated with theta phase precession in the subiculum.

## Discussion

We gained several insights into how spatial information is propagated and transformed through the subiculum. First, we found that the spatial representation within the subiculum exhibits a proximal-distal gradient, with higher firing rates and more spatially distributed firing patterns in the distal part of the subiculum. Second, we found that this transverse gradient within the subiculum has important functional consequences: the distributed firing-rate maps of neurons in the distal subiculum contained more information, on average per neuron, about both spatial location and spatial context than the sparse firing-rate maps of neurons in the proximal subiculum or area CA1. Finally, we demonstrated that neurons at all proximal-distal locations within the subiculum exhibit theta phase precession.

### Proximal-distal differentiation within the subiculum

We found several significant proximal-distal differences within the subiculum. Neurons in the proximal subiculum (closer to area CA1) have lower mean firing rates and sparser spatial firing-rate maps than neurons in the distal subiculum (closer to the presubiculum). The proximal-distal difference in the sparseness of coding is larger than the difference that was reported in a earlier study (Sharp and Green, 1994). This discrepancy may be due to methodological differences, including the septal-temporal recording location within the subiculum, the behavioral paradigm and the measure that was used to quantify sparseness.

We also found that the transverse gradient in spatial representation corresponded to a gradient in the spatial information. Neurons in the proximal subiculum, which are more similar to neurons in area CA1, carry less spatial information than neurons in the distal subiculum, which have high mean firing rates and large spatial activity fractions. In addition, neurons in the distal subiculum tend to show stronger theta modulation than neurons in the proximal subiculum.

The proximal and distal portions of the subiculum send output projections to different territories (Naber and Witter, 1998). Thus, the observed proximal-distal gradients in firing rate and spatial activity may reflect different requirements for sparse versus distributed spatial coding in the target areas. Proximal-distal gradients may also relate to the topographical organization of entorhinal inputs to the subiculum. The distal subiculum connects reciprocally with the medial EC while the proximal subiculum is interconnected with the lateral EC (Witter, 2006). LFP theta oscillations and theta modulation of spiking are stronger in the medial EC than in lateral EC (Deshmukh et al., 2010). This is consistent with our observation that neurons in distal subiculum show stronger theta modulation of spiking than neurons in proximal subiculum. At the same time, neurons in medial EC exhibit greater spatial selectivity than those in lateral EC (Hargreaves et al., 2005). We suggest that the dorsal subicular neurons that receive medial EC inputs may integrate a greater number of inputs from EC or area CA1 than neurons in the proximal subiculum, so that the net resultant spatial firing profile per neuron is less selective in the distal subiculum despite the spatial selectivity of individual neurons in medial EC. Finally, recent work has found an analogous proximo-distal functional specialization within the CA1 field (Henriksen et al., 2010). These findings point toward the existence of a general proximo-distal organization governing hippocampal computations.

### Advantages of distributed coding in the subiculum

We replicated previous findings that neurons in the subiculum tend to have higher firing rates and poorer spatial selectivity than neurons in area CA1 (Barnes et al., 1990; Sharp and Green, 1994). However, using information theory, we demonstrated that subicular place fields with high mean firing rates and broad, multipeaked spatial firing profiles actually convey more spatial information than the more compact CA1 place fields. In effect, subicular neurons with multiple unitary fields act as “multiplexed” CA1 place cells, sacrificing spatial selectivity for higher information capacity.

Previous studies have noted that spatial representations become progressively less sparse at each stage of feedforward processing (Barnes et al., 1990; Jung and McNaughton, 1993; Leutgeb et al., 2004). We propose that one of the functions of the subiculum is to transform the still relatively sparse spatial representation in area CA1 (Thompson and Best, 1989; Karlsson and Frank, 2008) into a dense, distributed representation that is suitable for conveying information to downstream targets outside of the hippocampal formation. Sparse representations have several advantages: metabolic efficiency, reduced interference during learning, associative memory capacity and the ability to learn using simple local synaptic plasticity rules (Marr, 1971; Fiete et al., 2004; Olshausen and Field, 2004). However, sparse representations are also wasteful in terms of average information per neuron, because a large fraction of neurons in the population are silent for most of the time and therefore contribute little information. In contrast, in dense distributed representations any randomly selected neuron in the population is likely to be informative at any given time. From the perspective of a downstream decoding circuit that can sample only a small random subset of incoming axons, it may advantageous to receive inputs from less sparse subicular neurons because enough spikes will be received in a short time window to disambiguate the animal's current spatial location. Thus, the distributed spatial representation in the subiculum may be specialized for conveying hippocampally processed information to the rest of the brain.

Intriguingly, these results are qualitatively consistent with novel theories of long range brain communication based on compressed sensing (Isely et al., 2011). This theory shows how a downstream region (in this context the subiculum) can randomly subsample a sparse high dimensional representation (in this context CA1) through dense synaptic connections, yielding a compressed representation with minimal information loss. This compressed representation can then be sent to other brain regions with a small number of long-range axons.

### Subicular representation of spatial context

An earlier hypothesis was that the subiculum, unlike area CA1, provides an “invariant” spatial representation (Sharp, 2006). Our data suggest subicular representations are more informative about spatial context than earlier suspected. We found that individual neurons in the subiculum remap across two different environments, even when those environments are geometrically identical. Furthermore, subicular and CA1 neurons convey similarly information about the identity of the current environment.

However, we did observe a diversity of remapping patterns in the subiculum. Some neurons appeared to have completely different place fields in the two environments, whereas others seemed to have place fields in similar locations but with modest scaling of firing-rate. These different degrees of remapping at a single-neuron may support multiple levels of representation, as has been suggested by others (Leutgeb et al., 2005). Neuron that maintain the same firing-rate map may support a generalized representation of geometric similarities across environments, whereas neurons that exhibit completely different firing-rate maps could disambiguate between geometrically identical environments (Singer et al., 2010).

### Theta phase precession in the subiculum

We found that neurons at all proximal-distal locations within the subiculum exhibit robust theta phase precession, and that the range of oscillatory spiking frequencies is similar between the subiculum and nearby regions of area CA1. This finding implies that the mechanism of phase precession in the subiculum must be remarkably invariant to proximal-distal gradients in mean firing rate, spatial activity, and anatomical connectivity.

The manifestation of theta phase precession across both area CA1 and the subiculum, and the “compositional” appearance of subicular firing-rate maps with multiple unitary place fields might seem to suggest that neurons in the subiculum inherit their properties as a result of the convergence of multiple inputs from CA1 place cells (Barnes et al., 1990; Sharp and Green, 1994). This hypothesis, although appealing, is unlikely to be correct. The projection from CA1 to the subiculum is extremely dense (Amaral et al., 1991; Cenquizca and Swanson, 2007), so each neuron in the subiculum likely receives inputs from many CA1 place cells. Recent modeling work has shown that the summation of spiking outputs from phase-precessing CA1 place cells results in a net excitatory drive that is phase-locked, not precessing, relative to theta (Geisler et al., 2010). Therefore, massed convergence and summation of many CA1 inputs would produce weak or absent phase precession in the subiculum. More intricate circuit-level mechanisms are required to explain the persistence of theta phase precession despite the tremendous convergence of CA1 inputs onto subicular neurons.

What are the consequences of theta phase precession for spatial information coding in the subiculum? It is not clear whether temporal coding of distances by “temporal sequence compression,” which was developed for hippocampal place cells with single place fields (Skaggs et al., 1996; Dragoi and Buzsáki, 2006), is a plausible mechanism in the subiculum. Neurons in the middle and distal subiculum tend to multiplex their firing over multiple place fields, so that their spikes do not correspond to a single region in the environment. As a result, spike-timing differences between pairs of neurons are ambiguous with respect to spatial order and distances between the neurons' place fields. Instead, it may be more fruitful to think about theta phase precession as a signature of transiently synchronized cell assemblies that are temporally segregated from the larger population by their frequency of co-oscillation (Geisler et al., 2010). Coordinated phase precession of cell assemblies across area CA1 and the subiculum may facilitate the transmission of spatial information and spike timing-dependent plasticity within the hippocampal circuit.

## Footnotes

This work was supported by NIH Grant MH080283, the Burroughs-Wellcome Foundation, and the John Merck Foundation. We thank Sheri Harris, Maggie Carr, Annabelle Singer, Mattias Karlsson, and Caleb Kemere for advice and technical support. We also thank Mary Dallman and Peter Ohara for use of their cryostat microtome, and the Nikon Imaging Center at University of California, San Francisco for providing microscopy facilities.

- Correspondence should be addressed to Loren M. Frank, W.M. Keck Center for Integrative Neuroscience and Department of Physiology University of California, San Francisco, 513 Parnassus Avenue, S-859, Box 0444, San Francisco, CA 94143-0444. loren{at}phy.ucsf.edu