Boundary Vector Cells Encode a Future-Biased Spectrum of Positions in the Rat.

Introduction

Cellular and oscillatory activity in the hippocampal formation exhibits robust correlations with navigational states, including position (O'Keefe and Dostrovsky1971; Ekstrom et al.2003; Ulanovsky and Moss2007; Hafting et al.2005), head direction (Taube et al.1990), boundary proximity/direction (Lever et al.2009; Solstad et al.2008; Alexander et al.2020), and locomotion, including running speed (Sargolini et al.2006; Caplan et al.2003; Wells et al.2013). These correlations are typically derived by comparing the neural variable in question to the subject's instantaneous state, such as the firing rate of a place cell to the instantaneous position of the rat locomoting in an open field. Intriguingly, however, such neural correlations are often enhanced when spiking activity, for example, is aligned not with the animal's instantaneous state, but rather with its time‐shifted or time‐integrated state (Muller and Kubie1989; Blair and Sharp1995; Dannenberg et al.2019; Bright et al.2020; Chaudhuri‐Vayalambrone et al.2023).

For example, the spatial tuning of CA1 neurons is more precise if aligned with the animal's position around 120 ms after the spike occurred, instead of its position at the time of the spike (Muller and Kubie1989; Chaudhuri‐Vayalambrone et al.2023). Similarly, spatial tuning in the subiculum (Sharp1999) and in grid cells of the medial entorhinal cortex (Chaudhuri‐Vayalambrone et al.2023) is enhanced by aligning spikes to future positions. These observations are consistent with the idea that the entorhinal‐hippocampal system is generating predictions about future states (Stachenfeld et al.2017; Momennejad2020; Geerts et al.2023). However, not all coding in the entorhinal‐hippocampal system appears to be prospective. For instance, entorhinal speed cells, though obviously linked to locomotion, provide a retrospective record, encoding the historical running speed of the animal (Dannenberg et al.2019). Furthermore, cells in the lateral entorhinal cortex also track the recent past, coding elapsed time since salient event boundaries (Bright et al.2020; Tsao et al.2018). Accordingly, it is not straightforward to suggest that entorhinal‐hippocampal codes are simply predictive.

Beyond the prospective‐retrospective dimension, a further important dimension of functional interest concerns multiple timescales. For instance, considering prospective coding, entorhinal grid cells encode the rat's future position with temporal offsets proportional to their spatial scale, thereby providing a multiscale representation of the animal's trajectory (Chaudhuri‐Vayalambrone et al.2023). Considering retrospective coding: Speed cells encode the historical running speed of the animal at widely differing timescales ranging from hundreds of milliseconds to hundreds of seconds, depending on the neuron (Dannenberg et al.2019); and lateral entorhinal cells encode the amount of time that has passed since critical event boundaries at multiple timescales (Bright et al.2020; Tsao et al.2018).

Theoretical modeling shows that multiscale encoding provides a rich basis upon which to form memories and predict future states (Howard et al.2014; Shankar et al.2016). A key prediction of that modeling is the presence of a spectrum of encoded timescales. Such a spectrum enables encoding of a memory timeline ofwhat happened whenand forming an estimate of the timeline of the future. The multiscale encodings observed in entorhinal neuronal activity (e.g., Tsao et al.2018; Dannenberg et al.2019; Bright et al.2020; Chaudhuri‐Vayalambrone et al.2023) exemplify functional types of spectral encoding. Importantly, in this theory, the type of information encoded determines the types of inferences obtainable. For example, multiscale encoding of time since event boundaries supports retrospective reconstruction of a timeline.

Boundary vector cells (BVCs) are spatially tuned neurons found in the subiculum that exhibit increased firing in the presence of environmental boundaries (e.g., wall or drop edge) or sufficiently large objects at specific allocentric distances and directions (Lever et al.2009; Stewart et al.2014; Poulter et al.2021). The existence of such cells was first predicted by computational models of inputs to place cells (O'Keefe and Burgess1996; Hartley et al.2000; Lever, Burgess, et al.2002) seeking to explain how place cell firing fields are shaped by environmental geometry (O'Keefe and Burgess1996; Lever, Wills, et al.2002). BVCs show some plasticity: for example, their firing fields can be influenced by the (e.g., rectangular) geometry of the environment (Muessig et al.2024); a subtype of BVCs, called vector trace cells, also exhibitmemoryfor the locations of boundaries and objects (Poulter et al.2021). See (Wang and Bicanski2026) for a computational model of how vector trace cells may emerge. Generally, less is known about subicular BVCs than about other cell types in the entorhinal‐hippocampal system, see (Poulter et al.2018; Bicanski and Burgess2018,2020; Wang and Bicanski2026; Lever et al.2026) for reviews and discussion.

Whether subicular BVCs encode boundary vectors retrospectively or prospectively, and whether this code operates over a range of timescales, are open questions. To address these questions, we first applied rate‐map methods, like those used to reveal multiscale temporal coding in grid cells (Chaudhuri‐Vayalambrone et al.2023), to estimate the temporal offset between BVC spiking and position for each cell that maximizes spatial tuning in the rate‐map. In a complementary approach, we modeled the spiking of each BVC as a function of a time‐shifted and/or time‐integrated representation of the rat's position. Model variants with a free parameter for time‐shift test whether BVC spiking is more predictable when the position was shifted forward or backwards in time and, if so, provide insight into the offsets that maximize that predictability. Model variants with a free parameter enabling time‐integration of position test whether BVC spiking is most consistent with a narrow (e.g., instantaneous) estimate of the rat's position or whether it tracks a position estimate obtained by averaging over a temporal window with some width. If the latter, these models also provide insight into the nature and scale of that temporal horizon. Our results, across analysis approaches, demonstrate that BVC spiking is best explained as time‐integrated estimates offutureposition but that, across cells, there is a spectrum of temporal offsets and temporal scales encoded.

Methods

New analyses were performed on recordings obtained by and first described by Lever et al. (2009). Subjects, surgery, and data collection are paraphrased here for ease of reference.

Subjects and Surgery

Data were obtained from six male Lister Hooded rats, weighing 315–390 g at the time of surgery. Rats were maintained on a 12:12 h light:dark schedule (with lights off at 15:00). Food deprivation was maintained such that subjects weighed 85%–90% of free feeding weight. Recordings were made from the dorsal subiculum. Details of the implantation and confirmation of electrode placement were described previously (Cacucci et al.2004; Lever et al.2009). Briefly, a bundle of four HM‐L coated platinum‐iridium wire tetrodes mounted to a microdrive was chronically implanted dorsal to the subiculum under deep anesthesia. Tetrodes were lowered to isolate cells after surgery.

Electrophysiology, Location and Head Orientation

Tetrodes were lowered toward the subiculum region and then left to stabilize before the start of the recording. The head stage was connected to a pre‐amplifier by 3‐m lightweight wire. The outputs of the pre‐amplifier passed through a switching matrix, and then to the filters and amplifiers of the recording system. Each channel was continuously monitored at a sampling rate of 50 kHz and action potentials were stored as 50 points per channel whenever the signal from any of the 4 channels of a tetrode exceeded a given threshold. EEG signals were amplified 10–20 K, band‐pass filtered at 0.34–125 Hz and sampled at 250 Hz (See details in Lever et al.2009).

Head position and orientation were tracked using an overhead video camera and tracking hardware/software by tracking the position of two arrays of head‐mounted small, infrared LEDs, one array brighter and more widely projecting than the other. Cluster cutting to isolate single units was performed manually using custom made software (TINT, Axona, UK). (See details in Lever et al.2009). For all analyses, the complete record of position was used. Though there is precedent for exclusion of data from epochs of low or no velocity, in the context of time shifting and time integrating behavior, this practice would have introduced undue complexity to the methods and reduced overall transparency.

Experimental Design and Training

Recordings were collected as rats foraged for cereal in multiple distinct open field arenas. Rats were transferred to the arenas in a consistent fashion to maintain orientation. At the end of each trial, the rat was removed from the recording environment and placed back on the holding platform until the next trial. The length of the trial was proportional to the size of the arena and was in the range 10–16 min. The inter‐trial interval was 20 min. The arena specifics are described by Lever et al. (2009). In all recordings, an external white cue card (102 cm high, 77 cm wide) provided directional constancy.

Rate Map Creation

Firing rate maps were constructed from 2.6 × 2.6 cm binned data, smoothed using a Gaussian smoothing kernel with a standard deviation of 5.2 cm. Spike count was divided by occupancy time to obtain a firing rate per bin. To display firing rate maps in figures, the rate maps were auto‐scaled false color maps, each color representing a 15% band of peak firing rate, from dark blue (0%–15%) to red (85%–100%). Peak rate (after smoothing) is shown next to each rate map.

Rate Map Tuning Based Approaches Estimating Time‐Shift

Following previously used approaches for estimating the temporal lag of spatial coding, we examined results from spatial information (s.i.) (Skaggs et al.1993) and zero‐lag spatial autocorrelation (ZLAC) (Chaudhuri‐Vayalambrone et al.2023) approaches with variable temporal offsets between the tracking and spiking. S.i. is computed as,

$$s . i . = \underset{n}{\sum} p_{n} \frac{r_{n}}{\bar{r}} log_{2} \left(\frac{r_{n}}{\bar{r}}\right)$$

where$p_{n}$is the probability of occupying bin$n$,$r_{n}$is the smoothed mean firing rate at bin$n$, and$\bar{r}$is the mean overall firing rate. This form of s.i. has the units of bits per spike. Note, however, that because the optimization was performed for each recording separately, preserving the recording duration and total spikes, it would not have changed the qualitative results had it been performed with the bits per second variant. ZLAC in this context is equivalent to the mean squared value of each rate map, computed as,

$$\text{ZLAC} = \frac{1}{N} \underset{n}{\sum} r_{n}^{2}$$

where$r_{n}$is the smoothed mean firing rate at bin$n$(i.e., the nth bin of the rate map). The maximum s.i. and ZLAC scores are obtained when spikes are least uniformly distributed over spatial bins.

In both approaches, each recording for each cell was analyzed separately. For each, we shifted the timing of the spiking relative to the tracking one video frame (20 ms) at a time for up to 2 s into the past and into the future. For each step, we created a rate map and then computed the s.i. and ZLAC scores. Over all steps, a curve of scores was constructed. These curves were then smoothed with a boxcar kernel with a width of 6 (120 ms). Representative rate maps and curves are shown in Figure1. Significance of s.i. and ZLAC at each lag was established by performing 500 pseudo‐random rotations between the tracking and spiking to establish the null distributions of s.i. and ZLAC for each recording. The rotations were constrained to shift the relative alignment of the data types by at least 5% of the recording length (30 s in the case of the 10 min recordings). The ‘rotation’ here is the wrapping of the misaligned data from one end of the recording around to the other. The permuted ZLAC and s.i. scores were used to compute a z‐score type normalization of the empirical ZLAC and s.i. values, respectively. Normalized values were considered to reflect significant spatial tuning if they were larger than 3.4808, reflecting$p < 0.00025$(the Bonferroni corrected$p < 0.05$over 200 comparisons). From each of the normalized curves with significant values, we identified the local maxima closest to zero as the ‘peak shift’ for that recording. The peak shifts across different recordings of the same cell were then averaged to establish a per‐cell‐estimate so as to perform statistics across cells instead of recordings.

Time shifting spikes relative to rat position improves spatial tuning. (A) Normalized (z‐scored) spatial information (s.i.; solid line) and normalized zero‐lag auto‐correlation(ZLAC; dashed line) for five representative BVCs are shown, plotted as a function of the temporal shift between the spiking and spatial position. Black dots mark the lag that resulted in the peak value for each cell. (B) Spatial firing ratemaps for a subset of the points shown for each of the aligned curves in A. A fixed color map is used across all ratemaps scaled to the peak firing rate observed at any shift, indicated at the right of each row. (C, D) same as A‐B but for representative non‐BVCs recorded in the subiculum.

Maximum Likelihood Estimation of Spatio‐Temporal Tuning

A concern with the approach described above is that it maximizes how non‐uniform the rate map is without consideration for whether that resulting rate map accurately captures the firing statistics of the neuron. To address this, we performed a complementary analysis approach that modeled the spiking of each cell as a function of a time‐shifted and/or time‐integrated representation of the rat's position. In this approach, temporal tuning parameters were selected so that the observed firing in a recording was maximally likely given knowledge of the rat's trajectory. Model variants with a free parameter for time‐shift test whether BVC spiking is more predictable when the position was shifted forward or backwards in time and, if so, provide insight into whether the neuron has prospective or retrospective coding. Model variants with a free parameter enabling time‐integration of position test whether BVC spiking is most consistent with a narrow (e.g., instantaneous) estimate of the rat's position or whether it tracks a position estimate obtained by averaging over a temporal window with some width. If the latter, these models also provide insight into the nature and scale of that temporal horizon.

In this approach, the quality of a model using a given set of parameters was assessed as follows:

$$\left(\hat{p}\right)_{\text{spike}} = 1 - e^{- r \left(\left(\overset{\sim}{x}\right)_{t}\right) \cdot \mathit{dt}}$$

where$r \left(\overset{\sim}{x}\right)$is the rate map value at surrogate position$\overset{\sim}{x}$and$\mathit{dt}$is the bin‐width used in discretizing the predictions, set to 0.02 s (one frame of the tracking).

$$\begin{matrix}\mathcal{L} \left(y , \left(\hat{p}\right)_{\text{spike}}\right) = - \underset{t}{\sum} \left[y_{t} log \left(\left(\hat{p}\right)_{\text{spike} , t}\right) + \left(1 - y_{t}\right) log \left(1 - \left(\hat{p}\right)_{\text{spike} , t}\right)\right]\end{matrix}$$

where the sum is over all time bins$t$in the recording,$y$is the vector of observed binary spike data (where$y_{t} \in \left\{0 , 1\right\}$), and$\left(\hat{p}\right)_{\text{spike}}$is the corresponding vector of the model's predicted spike probabilities. This calculation yields a single scalar value that quantifies the goodness‐of‐fit for a given set of parameters,$\tau$and$\kappa$. Lower$\mathcal{L} \left(y , \left(\hat{p}\right)_{\text{spike}}\right)$values indicate better predictions of the observed spike train.

This approach enables testing of whether a given parameter significantly improves how well the rate map obtained with given values of$\tau$and$\kappa$predicts the observed spiking. Broadly, this is done by testing whether$\mathcal{L} \left(y , \left(\hat{p}\right)_{\text{spike}}\right)$is significantly better when a parameter is treated as a free‐parameter (i.e., one that can be optimized via a gradient descent algorithm) instead of being fixed to a default value. Details of statistical testing are described below.

Time‐Shift

The time‐shift parameter,$\tau$, temporarily offsets the spiking with respect to the behavioral tracking. A value of$0$corresponds to how standard rate maps, with no temporal shift, are constructed. Positive$\tau$values align spiking with future positions (i.e., predict spiking from where the rat will be). Negative values align spiking to past positions (predict spiking from where that rat was). We tested whether optimizing$\tau$generated rate maps that reliably predicted spiking better than those generated with no time‐shift, and if so, whether that offset was reliably offset toward the future or the past as described in Statistical Testing.

Time‐Integration

Time‐integration estimates each position based on a time‐integrated version of position over a recent window with width proportional to$\kappa$. The form of the time‐integration depends on the shape of the kernel used.

We compared three kernel shapes: boxcar, Gaussian, and exponential. Each was compared to a baseline model with no integration.

The generalized form for the smoothed estimate at time$t$is:

$$\overset{\sim}{x} \left(t;\tau,\kappa\right) = \frac{\underset{s \in \mathcal{W} \left(t,\tau\right)}{\sum} K_{\kappa} \left(s - \left(t + \tau\right)\right) x \left(s\right)}{\underset{s \in \mathcal{W} \left(t,\tau\right)}{\sum} K_{\kappa} \left(s - \left(t + \tau\right)\right)}$$

where$K_{\kappa} \left(x\right)$is a given kernel function with width controlled by the parameter$\kappa$. To prevent information leakage between training and test folds during cross‐validation, the integration window$\mathcal{W} \left(t,\tau\right)$was bounded to include only samples within 15 s of the shifted time point$t + \tau$. For all models, if a kernel's effective window extended beyond this bound, it was truncated to$\mathcal{W} \left(t,\tau\right)$and renormalized by the denominator in Equation (6).

Boxcar (Uniform)

This model implements the simpler hypothesis that a neuron is equally sensitive to all behavioral states within a defined window. This corresponds to a uniform kernel applied retrospectively from the time point$t + \tau$:

$$\left(\overset{\sim}{x}\right)_{\text{Boxcar}} \left(t;\tau,\kappa\right) = \frac{1}{\kappa} \sum_{k = 0}^{\kappa - 1} x \left(t + \tau - k\right)$$

where the window width,$\kappa$, is the number of time bins included in the average. During optimization,$\kappa$was bounded to 15,000 ms.

Gaussian

This model tests the hypothesis that neural activity is maximally sensitive to behavior at latency$\tau$, with decreasing sensitivity for times further from this point. This is implemented with a Gaussian kernel:

$$K_{\kappa} \left(u\right) = exp \left(- \frac{u^{2}}{2 \sigma^{2}}\right) , \text{where} \sigma \equiv \kappa$$

This symmetric kernel integrates states both before and after the central time point$t + \tau$. Such a formulation is consistent with theoretical work on scale‐invariant representations of time that can be adapted to encode space O'Keefe and Burgess (1996). During optimization, the standard deviation$\sigma$was bounded to$\left[10 , 100000\right]$ms.

Exponential

To capture potentially asymmetric temporal sensitivity (i.e., a retrospective or prospective bias), we used a one‐sided exponentially decaying kernel. The model's direction (past vs. future) and decay rate were governed by a single fitted parameter,$\omega$, which was bounded to$\left[- 12 , 12\right]$. This parameterization avoids the numerical instability that can arise when fitting the direction and decay rate separately. Specifically, we mapped$\omega$to a decay constant$\kappa$and a direction$d$:

$$\kappa = \frac{1}{\mid \omega \mid + \epsilon} , d = sign \left(\omega\right) \in \left\{- 1,+ 1\right\}$$

where a small constant$\epsilon$ensures a minimal decay rate. The estimated behavioral state is then given by:

$$\left(\overset{\sim}{x}\right)_{Exp} \left(t;\tau,\omega\right) = \frac{\sum_{k = 0}^{K_{max}} \lambda^{k} x \left(t + \tau + d k\right)}{\sum_{k = 0}^{K_{max}} \lambda^{k}} , \text{where} \lambda = exp \left(- \frac{\Delta t}{\kappa}\right)$$

here$k$indexes bins,$K_{max} = \left⌊\frac{15000 ms}{\Delta t}\right⌋$to set the number of bins to sum over to enforce the 15,000 ms bound,$d = + 1$applies the kernel to future samples (prospective model) and$d = - 1$applies it to past samples (retrospective model).

Parameter Optimization

To optimize the free parameters for each model, we used a particle swarm gradient descent approach. Specifically, we used the Matlab functionparticleswarm.mwithSwarmSizeset to 50 andHybridFcnset tofmincon, which implements a local gradient‐based search after the initial swarm exploration. To ensure robust fitting, a random restart procedure was used wherein the particle swarm was reinitialized iteratively until five consecutive restarts failed to produce an improved fit. The parameter set that yielded the lowest final negative log‐likelihood ($\mathcal{L} \left(y , \left(\hat{p}\right)_{\text{spike}}\right)$) on the training data was kept for subsequent analysis.

Cross‐Validation

To test whether optimized parameters were generalizable to data not used in their estimation, and thereby support the critical evaluation of the goodness of the model fits, we used a K‐fold cross‐validation approach with$K = 10$. That is, recordings were split into 10 consecutive segments along the time variable. On each of 10 separate runs, a different segment was omitted from the fitting process and the remaining$K - 1$data segments were used for parameter estimation. The omitted segment was used as the validation dataset, used to assess the quality of the resulting fit.

Time‐shifting and time‐integrating procedures can blur the boundaries between the fitting and validation segments and introduce segments at the start or end of the recording without aligned behavior and spiking. To address this, buffers were cut from the start and end of the recording and around the validation dataset. Critically, a buffer of fixed size was used across all parameterizations. This was so that the number of omitted data points was constant across parameterizations and the likelihood score could not be affected by varying numbers of samples across parameter values.

Parameter Aggregation

Individual parameters were separately estimated for each cell dozens of times as a result of each of the recordings of that cell being subjected to a 10‐fold cross validation testing protocol. After establishing that the parameter estimate stability was good to excellent (see section below onParameter reliability with ICC), a single estimate of each parameter was obtained by taking the median value over all the estimates to enable analysis of parameter distributions across cells.

Statistical Testing

We employed different statistical approaches across analyses as follows: We used the non‐parametric Wilcoxon signed‐rank test to test for significant improvements in the rate map based analyses. For the maximum‐likelihood (ML) models, we used the Bayesian Information Criterion (BIC) for within‐dataset model selection. The BIC score for a model improves (reduces) in two circumstances: (1) the model is more accurate at predicting the data, and (2) the model uses fewer free parameters. Thus, it helps select for the models that are most parsimonious yet account for the data well. Population‐level effects were assessed with a Linear Mixed‐Effects (LME) model. Parameter estimate reliability was assessed with intra‐class correlations (ICC). Each is described here.

Rate Map–Based Analysis

The Wilcoxon signed‐rank test (reported as a$W_{\mathit{sr}}$statistic) was used to compare distributions to zero or to perform paired tests. The Wilcoxon rank‐sum test (reported as a$W_{\mathit{rs}}$statistic) was used for testing between independent samples.

Maximum‐Likelihood Model Comparison

With the goal of selecting a method of building rate maps that were maximally predictive of spiking that was not overly complicated, we compared different models on the basis of their Bayesian Information Criterion (BIC) scores. The BIC score for a model improves (reduces) in two circumstances: (1) the model is more accurate at predicting the data, and (2) the model uses fewer free parameters. It is computed as follows:

$$BIC = k ln \left(n\right) - 2 ln \left(\left(\hat{L}\right)_{\text{test}}\right)$$

where$k$is the number of free parameters,$n$is the number of time bins, and$\left(\hat{L}\right)_{\text{test}}$is the likelihood of the spiking observed in the test dataset under the optimized parameters. The number of free parameters was either zero in the case of the standard rate map, one models with either a time shift ($\tau$) or a time‐integration window ($\kappa$), and models with both terms added two.

To assess whether one parameter or the other reliably improved the BIC, a two‐factor repeated‐measures analysis of variance was conducted to examine the effects of adding a shift and varying the time‐integration approach. Time shift (2 levels: none vs. optimized) and kernel (4 levels: none, uniform, Gaussian, exponential) were included as fixed factors, and their interaction was tested to determine whether the effect of temporal shifting depended on kernel type. Because multiple estimates were obtained from the same neuron, neuron identity was included as a random intercept in a linear mixed‐effects framework to account for the non‐independence of repeated estimates within cells. This specification models cross‐validation and trial‐level observations while performing statistical inference at the level of independent neurons.

Population‐Level Inference With LME

To ask whether particular temporal components (time shift and/or the choice of integration kernel) reliably improved predictive performance across a cell population (BVCs or non‐BVCs), we fit a Linear Mixed‐Effects (LME) model to the cross‐validated test‐set negative log‐likelihood (NLL), reducing sensitivity to overfitting:

$$\mathit{nll}_{\text{test}} \sim T_{\text{shift}} \times T_{\text{kern}} + \left(\left(\bar{r}\right)_{spk} | \textit{cellID}\right)$$

where$T_{\text{shift}} \in \left\{0 , 1\right\}$indicates whether a free time‐shift parameter was included, and$T_{\text{kern}} \in \left\{\text{none},\text{boxcar},\text{exponential},\text{normal}\right\}$indexes the temporal‐integration kernel,$\left(\bar{r}\right)_{spk}$is the average spike rate for each recording, and cellID indicates which cell each recording is from. In practice, this asks whether our ability to predict the spiking of a cell for a withheld portion of data ($\mathit{nll}_{\text{test}}$) varied significantly when we optimized$T_{\text{shift}}$and/or used particular temporal‐integration kernel. Performing this LME fits beta weights for each term (i.e., how much$\mathit{nll}_{\text{test}}$changes as a function of each). Whether the beta weights are significantly different from zero indicates whether those terms had a significant impact on$\mathit{nll}_{\text{test}}$. Negative beta weights indicate that a term improved the fit. The$\left(\left(\bar{r}\right)_{spk}| \textit{cellID}\right)$term defines a spike rate and cellID as random effects. As a result, statistical inference is performed at the level of cells, rather than trials or fits, while still leveraging within‐cell information to improve estimation precision.

The specific equation used to perform the LME analysis is itself another model and can be more or less complex. In the form we use here, the model included a random slope for the mean firing rate of the cell,$\left(\bar{r}\right)_{spk}$, to factor out variance in$\mathit{nll}_{\text{test}}$related to firing rate before assessing the impact of each of the other terms. This is more complex than the simplest possible model,

$$\mathit{nll}_{\text{test}} \sim T_{\text{shift}} \times T_{\text{kern}} + \left(1|\textit{cellID}\right)$$

We used it nonetheless because the BIC was lower (BIC$= 2.69 \times 10^{5}$vs.$2.94 \times 10^{5}$), indicating that the cost of the extra term was outweighed by the improvement predicting the$\mathit{nll}_{\text{test}}$. The qualitative pattern of results reported here did not depend on including mean spike rate as a random effect.

Parameter Reliability With ICC

To establish whether parameter fits were reliable (i.e., consistent) across cross‐validation folds and recordings of the same cell, we used intra‐class correlation (ICC) (Shrout and Fleiss1979; Koo and Li2016). We used the form to test for one‐way random effects, absolute agreement, multiple raters/measurements (i.e., ICC(1,k) approach). To test for stability at the broadest level, all parameter estimates for a given cell (10 k‐fold fits for each of the separate recordings of the same cell) were treated as independent raters. To test for stability across recordings of the same cell, only one of the K‐fold estimates were used from each of the recordings (the first) for each cell. Per convention, values less than 0.5 indicate poor reliability, between 0.5 and 0.75 indicate moderate reliability, between 0.75 and 0.9 indicate good reliability, and greater than 0.9 are indicative of excellent reliability (Koo and Li2016). We found the parameter estimates had good to excellent reliability (i.e., broadly consistent) across the k‐fold cross‐validation runs and between separate recordings (see Table1).

Table: Intra‐class correlations (ICC) for the time lag and time integration parameter estimates across model types to assess fit reliability.

Results

BVC Spatial Tuning Is Greater for Future Positions

Our first approach to assessing if BVCs encode boundary vectors retrospectively or prospectively was to ask whether, for each, the spatial tuning observable in the rate map was more structured when the spiking was temporally offset in one direction or the other relative to the behavior. This was repeated at varying magnitudes of offset. Figure1shows representative examples of how the rate maps changed as a function of temporal offset (note that only a subset of the offsets tested are illustrated here). To quantify changes in spatial tuning across offsets, we used spatial information (s.i.) and the zero‐lag auto‐correlation (ZLAC), metrics that increase when spiking is less uniform. Specifically, we used the offset that yielded the maximum score to assess whether each cell encoded boundary vectors prospectively or retrospectively.

The mean offsets over recordings for each cell that generated peak s.i. values were reliably positive across the 35 BVCs (33 of 35 were positive, median [CI]$=$120$ms$[110,130], Wilcoxon sign‐rank$W_{\mathit{sr}} = 594 , Z = 5.07 , p < 0.001$). Likewise for ZLAC, the offsets that generated the maximum values were also reliably positive (35 of 35 were positive, median [CI]$=$150$ms$[120,170],$W_{\mathit{sr}} = 630 , Z = 5.16 , p < 0.001$). Thus, both s.i. and ZLAC had peaks at an offset of about one theta cycle into the future. These results are shown in Figure2. The fact that rate maps had the largest s.i. and ZLAC at positive offsets indicates that BVC spiking is most informative about future positions.

Temporal offsets that align spiking with future behavior yield rate maps with peak normalized spatial information (s.i.) and zero‐lag autocorrelation (ZLAC). (A) Mean normalized s.i. as a function of temporal lags plotted separately for all recordings from each BVC with at least one significant normalized s.i. (black) and recordings without significant normalized s.i. (gray). Positive aligns spikes with future positions. (B) Mean standardized s.i. for recordings with significant (black) and non‐significant (gray) spatial tuning. Confidence intervals indicate standard error. (C) Histogram of the temporal offsets that resulted in the peak normalized s.i. scores. (D–F) Same as A–C but for ZLAC instead of s.i.

The analyses above assume that offsets between BVC spiking and position are temporal (i.e., a fixed time shift) but the offset could, alternatively, be spatial (e.g., displaced by the distance between the head and the body). The key difference is that spatial shifts, but not temporal shifts, depend on running speed. The observed temporal offsets (120 to 150 ms) would, given the median running speed of 15.70 cm/s (95% CI = [14.87 16.85]), correspond to a spatial offset of 2.04 cm [1.93 2.19]. Testing what spatial shift, ranging from −20 to 20 cm in 0.5 cm increments, maximized s.i., we found the peak s.i. did not differ significantly from zero (12 of 35 were positive,$W_{\mathit{sr}} = 156 , Z = 1.43 , p > 0.05$). The median optimal shift was 0.0 cm (95% CI = [0.0 0.3]). ZLAC was largest when spiking was aligned with future positions (32 of 35 were positive,$W_{\mathit{sr}} = 555 , Z = 4.93 , p < 0.001$) albeit at larger offsets (median [CI] = 2.5 cm [2.0 3.5]). The results are shown in Figure3.

Spatial offsets that yielded peak normalized s.i. and ZLAC aligned BVC spiking with future positions. (A) normalized s.i. as a function of spatial lags plotted separately for all recordings from each BVC with at least one significant normalized s.i. (black) and recordings without significant normalized s.i. (gray). Positive aligns spikes with future positions. (B) Mean standardized s.i. for recordings with significant (black) and non‐significant (gray) spatial tuning. Confidence intervals indicate standard error. (C) Histogram of the spatial offsets that resulted in the peak normalized s.i. scores. (D–F) Same as A–C but for ZLAC instead of s.i.

Comparing the spatial tuning improvement magnitudes between the temporal versus spatial offsets, we found that the temporal offsets produced significantly larger improvements in s.i. (median difference [CI] = 0.92 bits [0.76 1.09],$W_{\mathit{sr}} = 560 , Z = 4.01 , p < 0.001$) and ZLAC (median difference [CI] = 0.55 [0.23 0.90],$W_{\mathit{sr}} = 442 , Z = 2.08 , p < 0.05$). These results indicate that the spatial offsets are not better than temporal offsets in terms of improving the apparent spatial tuning of the rate map.

To gain perspective on the temporal tuning of these BVCs, we examined the tunings of non‐BVCs that were co‐localized with the BVCs analyzed above. Of the 78 unique non‐BVCs co‐located with BVCs, 74 ($95 \%$) had rate maps with significantly greater s.i. than chance for at least one temporal offset. The temporal offsets that generated the greatest s.i. scores for the non‐BVCs were, like the BVCs, reliably positive (60 of 74,$82 \%$, were positive, median [95% c.i.]$=$90 ms [60, 116],$W_{\mathit{sr}} = 2464 , Z = 6.13 , p < 0.001$). Using ZLAC to assess spatial tuning yielded a similar result. Of the 72 cells that had significant ZLAC scores at least one temporal offset ($92 \%$), the temporal offsets that generated the greatest ZLAC scores were reliably positive (63 of 72,$92 \%$were positive, median [95% c.i.]$=$130 ms [110140],$W_{\mathit{sr}} = 2442 , Z = 6.67 , p < 0.001$). Comparing the median optimal shifts from the BVCs and non‐BVCs, we found no significant differences for either s.i. ($W_{\mathit{rs}} = 3780 , Z = - 1.88 , p = n . s .$) or ZLAC ($W_{\mathit{rs}} = 3657 , Z = - 1.53 , p = n . s .$). See Figure4for a comparison. Thus, the BVCs share a similar magnitude bias toward future coding with other non‐BVC neurons of the subiculum, indicating coherent temporal coding across the area.

Spatial tuning in BVCs and non‐BVC ratemaps was maximized with similar temporal shifts. (A) Normalized histograms of shifts that yielded the maximum ratemap spatial information (s.i.) are shown for non‐BVCs (blue) and BVCs (red). (B) Same as A, but for shifts that yielded the maximum zero‐lag auto‐correlation (ZLAC). (C) Box‐and‐whisker plot comparisons of non‐BVCs BVCs in terms of the shifts that yielded the peak s.i. score. No significant difference was found. (D) Same as C, but for shifts that yielded the peak ZLAC score. No significant difference was found.

BVC Spiking Is Better Predicted by Future Position

The results above, showing that spatial tuning in rate maps increases when spiking is aligned with future positions, do not address whether those rate maps are more accurate descriptions of the spiking. Addressing this, we tested whether rate maps constructed with a temporal shift between the behavior and spiking were better able to predict neuronal spiking than the standard, non‐shifted rate maps.

For each recording, we estimated the optimal temporal shift and constructed a corresponding rate map using 90% of the data and then tested how well that rate map predicted spiking in the remaining data. Prediction accuracy was quantified using Bayesian Information Criterion (BIC), a metric which rewards prediction quality while punishing model complexity and, thereby, assesses whether the extra complexity is warranted given the change in prediction quality. See methods for full details. BIC was significantly improved (lowered) by building rate maps with an optimized temporal offset (F(3, 15,250) = 5.21;$p$= 0.02). See Table2for the full set of results. That is, when we used 90% of a recording to estimate the optimal temporal offset and then built a rate map from that data with that offset, our ability to predict the spiking in the remaining 10% was significantly better than if the rate map was constructed without a temporal offset even after adjusting for the complexity increase associated with adding the temporal offset.

Table: Linear mixed‐effects modeling results of how BVC spike prediction depended on temporal lag and integration approach.

The temporal shifts that were selected by the optimization process were reliably positive (Median$= 84.5 ms$,$W_{\mathit{sr}} = 566 , N = 35 , Z = 4.11 , p < 0.001$; Figure5A). However, though the optimized offsets for most BVCs (24 of 35; 69%) were positive, for a subset of the BVCs (10 of 35; 29%) the optimized offset was effectively zero (i.e., within one frame of synchronized to the coincident tracking). One cell (3%) had a negative offset. This finding that the activity of 30% of the BVCs was most predictable without an offset applied could reflect different sub‐populations. With the idea that the subset with positive offsets may reflect a distinct subset of BVCs, we recomputed the median offset for that set after removing the BVCs with negative or zero offset and found it to be 146.0 ms [113.7161.5], approximating one theta cycle. Optimized offset estimates for each cell were reliable across cross‐validations and recordings (Table1). These results build on those in the prior section by showing that not only does aligning BVC spiking with future position improves spatial tuning, the resulting rate map also provides a more accurate model from which to predict BVC spiking.

Temporal offsets that create ratemaps with the maximum likelihood of generating the observed spiking are reliably positive. For visibility, the x‐axis was truncated to±2s. Offsets of outliers are indicated in the caption. (A) The histogram of time offsets over BVCs when no temporal integration is allowed. (B) Temporal offsets when the boxcar integration kernel was used. Points at −7.4 and −5.4 s are not visible. (C) Temporal offsets when a Gaussian integration kernel was used. Points at −7.3 and 3.8 s are not visible. (D) Temporal offsets when an exponential integration kernel was used. A point at −5.8 s is not visible.

Time‐Integration Further Improves BVC Spiking Predictions

Published evidence showing that entorhinal cells encode behavioral variables with a spectrum of timescales, with some neurons tracking rapid changes while others integrate over longer timescales (Chaudhuri‐Vayalambrone et al.2023), which motivated us to ask whether BVCs similarly encode position with a spectrum of timescales. Central to this line of inquiry is the idea that BVC spiking may not simply encode a specific instantaneous position (albeit at some lag). Rather, BVC spiking may encode states with a broad temporal horizon. To gain perspective on this, for each cell, we optimized a parameter controlling the temporal integration window width. As above, this was done using 90% of each recording and then tested on the remaining 10%.

Critically, the nature of the information being integrated and the method of integration affect how spatial information gets encoded by BVC spiking. Time‐integration kernels with different shapes (e.g., boxcar, Gaussian, or exponential) will result in different predicted spiking dynamics. The spiking of a neuron computing a simple moving average of position estimates would be well approximated by a boxcar integration kernel. Or, the spiking of a neuron computing the average over a population of noisy estimates would be well approximated by a Gaussian integration kernel. Finally, the spiking of a neuron computing temporally‐discounted estimates of future states, used in some reinforcement learning frameworks, would be well approximated by an exponentially decaying integration kernel. To handle this potential dependency, we estimated the temporal integration window separately using each of these kernels (boxcar, Gaussian, and exponential). We then compared the predictability of the BVC spiking as a population under each of the kernels. The full analysis used a two‐by‐four two‐way design, crossing temporal‐offset type [none, optimized] and temporal‐integration type [none, boxcar, Gaussian, exponential]. This permitted us to ask whether there were main effects of including an offset, a time‐integration window, and whether there was an interaction between the two on the ability to predict the observed spiking as quantified by BIC.

Figure6illustrates how manipulating the integration window and type of integration kernel affects the nature of the encoded position estimate. Two example BVCs are shown, one for which the optimized integration window was relatively short (Figure6A) and one for which the integration window was long (Figure6B). Trajectory plots and rate maps are shown for wide ranging integration windows for each cell. Temporal integration smooths out small fluctuations in the tracked position at short timescales. At long timescales, the trajectory tracks where the rat lingers.

Optimizing a time integration parameter improves the predictability of BVC spiking illustrated by two representative examples (A, B). For each A & B, (a) Top: Example trajectory plots when smoothed with a boxcar time‐integration kernel with seven different widths (as indicated above each plot). The black line shows the time‐integrated position under each parameter value. Red dots mark the estimated position under that integration parameter at the time of each spike. Bottom: The occupancy normalized rate map for the trajectory plot above each panel. Warm colors indicate higher firing rates. All panels are normalized to the same maximum rate, which is indicated at the far right of the row. (b) The negative log likelihood (NLL) for each of the seven sub‐panels shown in (a). Lower is better. The height of the lowest bar across (b), (e) and (h) is marked with a red line. (c) The trajectory (top) and rate map (bottom) for the optimized width computed by the maximum‐likelihood parameterization. (d–f) Same as (a–c) but for varying values of sigma for a Gaussian time‐integration kernel. (g–i) Same as (a–c) but varying values of tau for an exponential time‐integration kernel.

In all cells of this analysis design, the parameter estimates were highly reliable over cross‐validations and recordings for individual BVCs (Table1). Comparing BIC values in the two‐way design, we found significant main effects for both temporal offset type (as reported above, see Figure5B–D) and temporal integration type (F(3, 15,250) = 183.27,$p < 0.001$) and found a significant interaction effect (F(3, 15,250) = 64.43,$p < 0.001$). Examination of the estimated effects of different levels of temporal offset and integration type on the BIC score, shown in Table2, reveal that both parameters made rate maps more predictive of spiking.

Improvements in BIC score obtained from adding a temporal integration were substantial across all three variants (integration improved BIC by$\sim$812–940). To test whether one kernel outperformed the others for predicting BVC spiking, we performed post hoc comparisons of BIC under each. The Gaussian kernel outperformed the other two, generating significantly lower BIC scores than the boxcar and exponential kernels (boxcar: est. = −94.6 [−185.0–5.3],t(15,432) = 2.08,$p < 0.05$; exponential: est. = −128.1 [−217.5–38.8],t(15,432) = 2.8,$p < 0.005$). There was no significant difference between the boxcar and exponential kernels (est. = 33.5 [−55.8, 122.8],t(15,432) = 0.73,$p = n . s .$). For the Gaussian kernel to better explain the variance in BVC spiking than the boxcar or exponential kernels suggests it is a better approximation of what the BVCs encode. The Gaussian kernel approximates a model wherein the BVC spiking encodes the average of noisy estimates of the rat's location.

The significant interaction effect was driven primarily by further improvements in BIC when the temporal offset is combined with the exponential integration approach. Neither of the other combinations resulted in further significant effects on BIC. The specificity of the interaction to the exponential kernel may be a consequence of this kernel placing exponentially greater weighting on the origin (i.e., the point set by the offset) over other time points. In contrast, neither the boxcar nor Gaussian kernels weight the origin so substantially.

The distributions of time integration kernel sizes are shown in Figure7. Across all three kernel types, a cluster of the kernels had small magnitudes, indicating that the most predictive ratemaps were generated by limited time integration. However, in each of the cases there was also considerable spread, including values that reflect integration over large temporal windows. The median [95%] sigma for the Gaussian kernels was 2511 ms [795 3307] and a coefficient of variation (CV) of 0.71. The median [95%] width for the boxcar kernels was 2905 ms [1006, 5768] and had a CV of 0.92. The median [95%] tau for the exponential kernels was 246 ms [173, 657] and had a CV of 1.43. The inherent differences across kernels regarding the meaning of these timescale parameters precludes direct comparison across kernels (e.g., the width of the boxcar defines 100% of the window whereas, for the Gaussian kernel, 95.45% of the window is defined by taking two‐sigma from before and after the center point). Yet, what is nonetheless clear is that there is substantial variability (CV$> 0.5$) in the temporal scales that position is encoded with across BVCs. Figure8illustrates what these parameters mean in the context of each of the kernels, making apparent the variability in temporal integration width that existed across the optimized kernels.

Optimized time integration parameters across kernel types. (A) Histogram of widths used in constructing boxcar kernels. (B) Sigma values used for constructing Gaussian kernels. (C) Tau equivalents of the omega parameter used to construct the exponentially decaying integration kernels. Positive values constructed kernels that decayed into the future, while negative values constructed kernels that decayed into the past.

Reconstructed time‐integration kernels across model variants. (A) Horizontal bars show the optimized boxcar kernel integration window width for each cell when no temporal offset was permitted. (B) False color heat maps show Gaussian kernel width for each cell when no temporal offset was permitted. Warmer colors indicate portions of kernel weighted more heavily. (C) Like B but for exponential kernels. (D–F) Same as A–C but for models with temporal offsets. In all large panels, the ordering of cells is the same. The inset panels plot the same data reordered based on the temporal offset and integration parameters for each specific model.

Time Shift and Integration Windows of BVCs and Non‐BVCs Are Similar

Finally, we asked if the temporal tuning properties were distinct between BVCs and colocalized non‐BVCs. Overall, we found the same qualitative pattern of results for non‐BVCs and found no significant differences in the median optimized parameter values. Precise statistics can be found in Table3. There was a significant main effect for temporal offset (F(1, 25,428) = 156.5;$p < 0.001$) and temporal integration type (F(3, 25,428) = 3517.3,$p < 0.001$) and a significant interaction effect (F(3, 25,428) = 27.4,$p < 0.001$). Estimated effects on BIC across all levels (see Table3) show that non‐BVC spiking was more predictable when an optimized temporal offset and integration window were used to model the spatial tuning.

Table: Summary of maximum likelihood results across model types for fits of non‐BVCs.

The optimized temporal offsets of the non‐BVCs, like the BVCs, were reliably positive (median = 122.7 ms [67.1166.7];W_{\mathit{sr}} = 1978 , N = 69 , Z = 4.6 ,$$p < 0.001). Of the 69 non‐BVCs colocalized with BVCs, the optimized temporal offsets of 51 (73.9%) were greater than 1 frame of behavioral tracking into the future ($>$20 ms), 9 (13.0%) were within 1 frame of zero offset, and 9 (13.0%) had offsets of more than 1 frame into the past. Comparing the distributions of optimized offsets between BVCs and colocalized non‐BVCs, we found no significant difference (102.9vs. 122.7 ms;W_{\mathit{rs}} = 1658 , Z = - 0.77 ,$$p = n . s .).

Comparing the different temporal integration kernels, as with the BVCs, we found that the Gaussian kernel yielded position estimates that made the spiking most predictable. BIC values for the Gaussian kernel were significantly lower than those from the boxcar and exponential kernels (Gaussian vs. boxcar: est. = −73.2 [−92.9–53.6],t(25,742) = 7.32,$p < 0.001$; Gaussian vs. exponential: est. = −137.4 [−157.0–117.8],t(25,742) = 13.73,$p < 0.001$). The boxcar kernel, in turn, produced BIC values that were significantly better than the exponential kernel (est. = −64.2 [−83.8–44.5],t(25,742) = 6.4,$p < 0.001$). The finding that the Gaussian temporal integration kernel made the spiking most predictable is consistent with non‐BVCs encoding the mean of a noisy population.

Across the non‐BVCs, there was considerable variability in the scale of the optimized time‐integration windows. For the Gaussian kernel, the median [95% CI] sigma value was 2571 ms [2363, 3674] and a coefficient of variation (CV) of 0.66 (data not shown). For the boxcar kernel, the median [95% CI] width was 5127 ms [3552, 6884] and a CV of 0.77. For the exponential kernel, the median [95% CI] tau (collapsing over forwards and backwards kernels) was 513 ms [298, 873] and a CV of 1.52. As was observed among the BVCs, there is high variability (CV$> 0.5$) in the integration windows across cells irrespective of which kernel type was used to estimate the temporal tuning width. This reflects the diversity of temporal encoding present in the subiculum.

Discussion

Spatial coding in the subiculum is less well understood than that in other entorhinal‐hippocampal regions such as CA1 and medial entorhinal cortex. See Place et al. (2026); Lee and Lee (2026), and Lever et al. (2026) for subiculum‐centred reviews of hippocampal spatial coding. The present study aimed to better understand the temporal coding properties of boundary vector cells (BVCs) recorded in the rat dorsal subiculum. To this end, we analyzed recordings of BVCs previously described by Lever et al. (2009) using two approaches. In the first, we asked what temporal offset between the rat's position and BVC spike timing yielded rate maps with the maximum apparent spatial tuning. Across BVCs, we found that positive offsets, aligning spikes to future positions, resulted in rate maps with the greatest spatial tuning. In the second approach, we modeled the spiking of each BVC as a function of the rat's position. The models could introduce a temporal offset between the behavioral tracking and spike timing and/or a time‐integrated representation of the position. The BVC spiking was best predicted by models with both the temporal offset and time‐integration of position. The offsets that resulted in the best predictions were reliably positive, aligning spikes with future tracking. The time‐integration windows that generated the best predictions varied widely across cells, ranging from 0.48 s up to 12.0 s. The variability from one cell to the next indicates that there is a multiscale encoding of future states across these BVCs. To establish whether these temporal coding properties are unique to BVCs among subiculum cells, we repeated the same analyses on non‐BVC cells recorded from the same tetrodes. The resulting offsets and parameter estimates were not significantly different between the BVCs and non‐BVCs, indicating that multiscale encoding of future states is likely a general property of the subiculum.

The two analysis approaches used here address related but distinct aspects of neural coding in the subiculum. The first analysis followed an approach that has been used previously to test for multi‐scale temporal coding in the entorhinal cortex (Chaudhuri‐Vayalambrone et al.2023). The analysis effectively asks, ‘At what temporal offset between the behavioral tracking and the spike timing is the correspondence between the neural activity and position least random?’ This is a computationally efficient analysis that effectively previewed the modeling results of the second approach. Yet, this approach does not address whether the offsets that maximize apparent spatial tuning also describe the spiking of the cells better. Describing the cell spiking is the primary objective of the second approach. The second approach used amodel selectionapproach to first ask what free‐parameters significantly improve our ability to predict the neural activity. Subsequent to establishing that a given parameter (e.g., the temporal offset) significantly improved our ability to model the spiking, we could examine the values of that parameter that were most likely to explain the observed spiking (e.g., offset aligned with future positions).

Another strength of the model selection approach is that it allowed us to test between alternate hypotheses for how neurons may encode position. By comparing models that predicted spiking by integrating the position with different kernels, we could test whether the observed spiking was most likely under one form of integration over another. We compared three kernel types: boxcar, Gaussian, and exponential. Each operationalizes a different implicit model for how the neurons integrate the behavioral state. A boxcar integration kernel weights all time points within a window equally. This captures an implicit model wherein the spiking of a neuron is the result of averaging over state estimates spanning a temporal range. The Gaussian integration kernel operationalizes an implicit model, related to the central limit theorem, wherein the spiking of a neuron results from sampling a noisy distribution. Finally, the exponential kernel operationalizes implicit models wherein the activity results from a process that discounts time points further away. Adding time‐integration to the models significantly improved our ability to predict the BVC spiking irrespective of kernel type. However, the Gaussian kernel outperformed the other two. This was the case for both the BVCs and the colocalized non‐BVCs. This is consistent with a mechanism wherein multiscale tuning in the subiculum emerges by each of the cells sampling a noisy distribution.

Given known heterogeneities among BVCs (e.g., trace vs. non‐trace; Poulter et al. (2021)), it is possible that there is variability in information being integrated by the different BVCs analyzed here. Though we lacked the power to subdivide the BVCs functionally, we did find that the spiking of a subset of the BVCs (10 of 35; 29%) was best predicted by models with no temporal offset despite that the spiking of most BVCs (24 of 35; 69%) was best predicted by future positions.

The multiscale temporal coding of behavior documented here is broadly consistent with what has been observed in the reciprocally interconnected entorhinal cortex. Entorhinal grid cells, speed cells, and cells of the lateral entorhinal cortex have all been found to exhibit multiscale encoding of state (Chaudhuri‐Vayalambrone et al.2023; Dannenberg et al.2019; Bright et al.2020; Tsao et al.2018). Yet, there are also important differences. The multiscale encoding of time in the lateral entorhinal cortex and of speed in the medial entorhinal cortex provideretrospectiverecords. Here, however, we find that BVCs areprospective. Entorhinal grid cells are also prospective, yet their hexagonal code affords different computations. Multiscale encodings among grid cells can, per established theory (Howard et al.2014; Shankar et al.2016), support inferential reasoning about trajectories between distal targets. Multiscale encoding of distance from boundaries could, by a similar reasoning, enable reasoning about the timing of events at a boundary.

What computations might an animal perform at environmental boundaries? One clue is that boundary encounters calibrate medial entorhinal grid cells, acting as an error‐correction signal (Hardcastle et al.2015). Such correction demands a calibrated estimate of distance to the boundary. A multiscale encoding of distance to a boundary would provide the necessary information for calculating this correction factor. Across a population of cells, there would be graded information about when a boundary was reached. In the moments leading up to and following an interaction with a boundary, there would be sufficient information across the population of BVCs to infer the relative trajectory taken. That estimate could then update the animal's current positional belief or perceived environmental scale.

The future bias of BVCs aligns them with the prospective coding of other spatially tuned neurons of the hippocampal formation (Muller and Kubie1989; Chaudhuri‐Vayalambrone et al.2023). Place cells, grid cells, and BVCs all tend to fire$100 - 200 ms$before the animal reaches a salient location. This prospective bias is in accord with the hypothesis that the hippocampus encodes a successor representation (Stachenfeld et al.2017; Momennejad2020; Geerts et al.2023). At its core, the successor‐representation view holds that hippocampal activity encodeswhat comes next—not a static map of space or context—conditioned on the current state and habitual transitions. Spike timing dependent plasticity mechanisms of the hippocampus are sufficient to establish successor representations (George et al.2023). The consistent direction and scale of prospective coding among spatially tuned neurons distributed across the hippocampal formation, including CA1, medial entorhinal cortex, and subiculum, point to the development of the successor representation being a circuit‐wide computation. Interesting directions for future work would be to test if this prospective coding develops as the internal model of an environment develops over familiarization, and would investigate the potential inter‐dependence of different circuit elements (e.g., CA3, CA1, subiculum) for prospective coding.

Multiscale encoding among spatially tuned neurons, including the BVCs examined here, in the context of a successor‐representation view of hippocampal function, may align with temporal discounting (Bos and McClure2013; McClure et al.2004; Masuda et al.2020). Neurons in CA1 encode time and other key parameters necessary for temporal discounting (Salz et al.2016; Masuda et al.2020), though it remains to be established that neurons that overtly encode spatial information are key contributors to decision making in the context of delay discounting. Temporal discounting is the idea that future states that will occur after a long delay are encoded more weakly than those that will occur soon but, importantly, are nonetheless encoded. This is, in effect, a way to interpret the multiscale integration windows we observed. Though many BVCs encoded boundaries that will be reached shortly, there were nonetheless BVCs that responded considerably earlier. As such, across the ensemble, there is information encoded about a spectrum of future horizons. This ensemble‐level of encoding is distinct from cellular‐level state encoding, which we failed to find support for when we tested whether exponential integration kernels could predict BVC spiking.

Departing from the successor‐representation view, however, we also observed substantial retrospective coding. In 15% to 30% of our fits, across approaches, we found that the temporal lags that made the empirical spiking most likely were offset into the past. Assuming that this does not reflect measurement error, this is not well accounted for by a standard successor‐representation perspective. Minor modifications could reconcile this mismatch. The Tolman‐Eichenbaum Machine (Whittington et al.2020), for example, contains elements of the Successor‐representation, including having an imposed bias toward predictive coding. Yet, it also includes associative learning that will form bidirectional associations between states. These associations could, in principle, reactivate prior states given a sensory cue, thereby introducing some retrospective coding.

The present work provided consistent evidence of future‐biased multiscale encoding of position based on the recordings from Lever et al. (2009), but there are several limitations important to acknowledge. One is that the present work did not overtly model the boundary‐oriented nature of the coding itself. Rather, the rate map‐based approach is agnostic to the form of spatial coding. Though this adds a level of generality, the generality runs the risk of missing key elements of the temporal tuning that may only appear in analyses optimized for boundary coding. Part of the reason for not taking that approach here is that accurate modeling of interactions with boundaries would require precise records of the timing of the moments of perceiving a boundary, which is not straightforward to obtain. It is also worth recognizing that there is no ground‐truth for the designations of BVC and non‐BVC. Consequently, cells may have been mislabeled in either direction. The risk of which, as it relates to the current work, is that temporal shifts observed for non‐BVCs may nonetheless reflect BVC coding and that some portion of the heterogeneity in temporal coding observed over BVCs was introduced by inclusion of non‐BVCs. Given the careful comparison of spatial coding over multiple sessions per cell performed by Lever et al. (2009), however, we do not expect this risk to be high, particularly with regard to the possible mis‐classification of non‐BVCs as BVCs. The current findings are unlikely to be merely due to somehow incorrect placement of LEDs upon the rat's head. Inspection of the timeshift curves for the BVCs (Figures1and2) shows the same phenomenon that (Muller and Kubie1989) observed for CA1 place cells. The rate of degradation of spatial tuning away from the optimal timeshift is not symmetric but rather is slower for shifts that are positive to the optimum (i.e., where spikes align with future positions). “We agree with (Muller and Kubie1989) view that this asymmetry is incompatible with merely incorrect LED‐placement.”

In conclusion, the current work showed through the use of multiple complementary analyses that BVCs are biased toward encoding future positions and that encoding is most consistent with a time‐integrated version of position. Across cells, we observed variability in the width of the window over which that integration was performed. From these findings, we conclude that BVCs encode a future‐biased spectrum of positions in the rat.

Funding

This work was supported by National Institutes of Health, R01AG076198; Biotechnology and Biological Sciences Research Council, BB/T014768/1.

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

We gratefully acknowledge support from the National Institutes of Health's National Institute on Aging, grant R01AG076198 to Z.T. and E.L.N. C.L. gratefully acknowledges grant funding from the BBSRC (BB/T014768/1). This research was also supported in part by Lilly Endowment Inc., through its support for the Indiana University Pervasive Technology Institute and by Shared University Research grants from IBM Inc. to Indiana University. ChatGPT, models 3.5 and 4 by OpenAI, were used to assist in the development and debugging of analysis code and for producing the LaTeX markup for equations and tables appearing in this manuscript.

Contributor Information

Ehren Lee Newman, Email: ehnewman@iu.edu.

Colin Lever, Email: colin.lever@durham.ac.uk.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Associated Data

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

1. Alexander, A. S. , Carstensen L. C., Hinman J. R., Raudies F., Chapman G. W., and Hasselmo M. E.. 2020. “Egocentric Boundary Vector Tuning of the Retrosplenial Cortex.” Science Advances 6, no. 8: eaaz2322. doi.org/10.1126/sciadv.aaz2322
2. Bicanski, A. , and Burgess N.. 2018. “A Neural‐Level Model of Spatial Memory and Imagery.” eLife 7: e33752. doi.org/10.7554/eLife.33752
3. Bicanski, A. , and Burgess N.. 2020. “Neuronal Vector Coding in Spatial Cognition.” Nature Reviews Neuroscience 34, no. 9: 171–470. doi.org/10.1038/s41583-020-0336-9
4. Blair, H. , and Sharp P.. 1995. “Anticipatory Head Direction Signals in Anterior Thalamus: Evidence for a Thalamocortical Circuit That Integrates Angular Head Motion to Compute Head Direction.” Journal of Neuroscience 15, no. 9: 6260–6270. doi.org/10.1523/JNEUROSCI.15-09-06260.1995
5. Bos, W. , and McClure S. M.. 2013. “Towards a General Model of Temporal Discounting.” Journal of the Experimental Analysis of Behavior 99, no. 1: 58–73. doi.org/10.1002/jeab.6
6. Bright, I. M. , Meister M. L. R., Cruzado N. A., Tiganj Z., Buffalo E. A., and Howard M. W.. 2020. “A Temporal Record of the Past With a Spectrum of Time Constants in the Monkey Entorhinal Cortex.” Proceedings of the National Academy of Sciences 117, no. 33: 20274–20283. doi.org/10.1073/pnas.1917197117
7. Cacucci, F. , Lever C., Wills T. J., Burgess N., and O'Keefe J.. 2004. “Theta‐Modulated Place‐By‐Direction Cells in the Hippocampal Formation in the Rat.” Journal of Neuroscience 24, no. 38: 8265–8277. doi.org/10.1523/JNEUROSCI.2635-04.2004
8. Caplan, J. B. , Madsen J. R., Schulze‐Bonhage A., Aschenbrenner‐Scheibe R., Newman E. L., and Kahana M. J.. 2003. “Human Theta Oscillations Related to Sensorimotor Integration and Spatial Learning.” Journal of Neuroscience 23, no. 11: 4726–4736. doi.org/10.1523/JNEUROSCI.23-11-04726.2003
9. Chaudhuri‐Vayalambrone, P. , Rule M. E., Bauza M., et al. 2023. “Simultaneous Representation of Multiple Time Horizons by Entorhinal Grid Cells and CA1 Place Cells.” Cell Reports 42, no. 7: 112716. doi.org/10.1016/j.celrep.2023.112716
10. Dannenberg, H. , Kelley C., Hoyland A., Monaghan C. K., and Hasselmo M. E.. 2019. “The Firing Rate Speed Code of Entorhinal Speed Cells Differs Across Behaviorally Relevant Time Scales and Does Not Depend on Medial Septum Inputs.” Journal of Neuroscience 39: 3434–3453. doi.org/10.1523/JNEUROSCI.1450-18.2019
11. Ekstrom, A. D. , Kahana M. J., Caplan J. B., et al. 2003. “Cellular Networks Underlying Human Spatial Navigation.” Nature 425, no. 6954: 184–188. doi.org/10.1038/nature01964
12. Geerts, J. P. , Gershman S. J., Burgess N., and Stachenfeld K. L.. 2023. “A Probabilistic Successor Representation for Context‐Dependent Learning.” Psychological Review 131, no. 2: 578–597. doi.org/10.1037/rev0000414
13. George, T. M. , Cothi W., Stachenfeld K. L., and Barry C.. 2023. “Rapid Learning of Predictive Maps With STDP and Theta Phase Precession.” eLife 12: e80663. doi.org/10.7554/eLife.80663
14. Hafting, T. , Fyhn M., Molden S., Moser M.‐B., and Moser E. I.. 2005. “Microstructure of a Spatial Map in the Entorhinal Cortex.” Nature 436, no. 7052: 801–806. doi.org/10.1038/nature03721
15. Hardcastle, K. , Ganguli S., and Giocomo L.. 2015. “Environmental Boundaries as an Error Correction Mechanism for Grid Cells.” Neuron 86, no. 3: 827–839. doi.org/10.1016/j.neuron.2015.03.039
16. Hartley, T. , Burgess N., Lever C., Cacucci F., and O'Keefe J.. 2000. “Modeling Place Fields in Terms of the Cortical Inputs to the Hippocampus.” Hippocampus 10, no. 4: 369–379. doi.org/10.1002/1098-1063(2000)10:4<369::AID-HIPO3>3.0.CO;2-0
17. Howard, M. W. , MacDonald C. J., Tiganj Z., et al. 2014. “A Unified Mathematical Framework for Coding Time, Space, and Sequences in the Hippocampal Region.” Journal of Neuroscience 34, no. 13: 4692–4707. doi.org/10.1523/JNEUROSCI.5808-12.2014
18. Koo, T. K. , and Li M. Y.. 2016. “A Guideline of Selecting and Reporting Intraclass Correlation Coefficients for Reliability Research.” Journal of Chiropractic Medicine 15, no. 2: 155–163. doi.org/10.1016/j.jcm.2016.02.012
19. Lee, S. , and Lee I.. 2026. “Encoding Patterned Spatial Behaviours in the Subiculum.” Hippocampus 36. doi.org/10.1002/hipo.70102
20. Lever, C. , Burgess N., Cacucci F., Hartley T., and O'Keefe J.. 2002. “What Can the Hippocampal Representation of Environmental Geometry Tell Us About Hebbian Learning?” Biological Cybernetics 87, no. 5: 356–372. doi.org/10.1007/s00422-002-0360-z
21. Lever, C. , Burton S., Jeewajee A., O'Keefe J., and Burgess N.. 2009. “Boundary Vector Cells in the Subiculum of the Hippocampal Formation.” Journal of Neuroscience 29, no. 31: 9771–9777. doi.org/10.1523/JNEUROSCI.1319-09.2009
22. Lever, C. , Poulter S., and Hollins A.. 2026. “The Subiculum, and Spatial and Episodic Memory.” Hippocampus.
23. Lever, C. , Wills T., Cacucci F., Burgess N., and O'Keefe J.. 2002. “Long‐Term Plasticity in Hippocampal Place‐Cell Representation of Environmental Geometry.” Nature 416, no. 6876: 90–94. doi.org/10.1038/416090a
24. Masuda, A. , Sano C., Zhang Q., et al. 2020. “The Hippocampus Encodes Delay and Value Information During Delay‐Discounting Decision Making.” eLife 9: e52466. doi.org/10.7554/eLife.52466
25. McClure, S. M. , Laibson D. I., Loewenstein G., and Cohen J. D.. 2004. “Separate Neural Systems Value Immediate and Delayed Monetary Rewards.” Science 306, no. 5695: 503–507. doi.org/10.1126/science.1100907
26. Momennejad, I. 2020. “Learning Structures: Predictive Representations, Replay, and Generalization.” Current Opinion in Behavioral Sciences 32: 155–166. doi.org/10.1016/j.cobeha.2020.02.017
27. Muessig, L. , Ribeiro Rodrigues F., Bjerknes T. L., et al. 2024. “Environment Geometry Alters Subiculum Boundary Vector Cell Receptive Fields in Adulthood and Early Development.” Nature Communications 15, no. 1: 982. doi.org/10.1038/s41467-024-45098-1
28. Muller, R. , and Kubie J.. 1989. “The Firing of Hippocampal Place Cells Predicts the Future Position of Freely Moving Rats.” Journal of Neuroscience 9, no. 12: 4101–4110. doi.org/10.1523/JNEUROSCI.09-12-04101.1989
29. O'Keefe, J. , and Burgess N.. 1996. “Geometric Determinants of the Place Fields of Hippocampal Neurons.” Nature 381, no. 6581: 425–428. doi.org/10.1038/381425a0
30. O'Keefe, J. , and Dostrovsky J.. 1971. “The Hippocampus as a Spatial Map. Preliminary Evidence From Unit Activity in the Freely‐Moving Rat.” Brain Research 34, no. 1: 171–175. doi.org/10.1016/0006-8993(71)90358-1
31. Place, R. E. X. , Sun Y., Olson J., and Nitz D.. 2026. “Subiculum Encodes Environmental Structure.” Hippocampus 2, no. 36: e70086. doi.org/10.1002/hipo.70086
32. Poulter, S. , Hartley T., and Lever C.. 2018. “The Neurobiology of Mammalian Navigation.” Current Biology 28, no. 17: R1023–R1042. doi.org/10.1016/j.cub.2018.05.050
33. Poulter, S. , Lee S. A., Dachtler J., Wills T. J., and Lever C.. 2021. “Vector Trace Cells in the Subiculum of the Hippocampal Formation.” Nature Neuroscience 24, no. 2: 266–275. doi.org/10.1038/s41593-020-00761-w
34. Salz, D. M. , Tiganj Z., Khasnabish S., et al. 2016. “Time Cells in Hippocampal Area CA3.” Journal of Neuroscience 36, no. 28: 7476–7484. doi.org/10.1523/JNEUROSCI.0087-16.2016
35. Sargolini, F. , Fyhn M., Hafting T., et al. 2006. “Conjunctive Representation of Position, Direction, and Velocity in Entorhinal Cortex.” Science (New York, N.Y.) 312, no. 5774: 758–762. doi.org/10.1126/science.1125572
36. Shankar, K. H. , Singh I., and Howard M. W.. 2016. “Neural Mechanism to Simulate a Scale‐Invariant Future.” Neural Computation 28, no. 12: 2594–2627. doi.org/10.1162/NECO_a_00891
37. Sharp, P. E. 1999. “Comparison of the Timing of Hippocampal and Subicular Spatial Signals: Implications for Path Integration.” Hippocampus 9, no. 2: 158–172. doi.org/10.1002/(SICI)1098-1063(1999)9:2<158::AID-HIPO7>3.0.CO;2-O
38. Shrout, P. E. , and Fleiss J. L.. 1979. “Intraclass Correlations: Uses in Assessing Rater Reliability.” Psychological Bulletin 86, no. 2: 420–428. doi.org/10.1037//0033-2909.86.2.420
39. Skaggs, W. E. , McNaughton B. L., Gothard K. M., and Markus E. J.. 1993. “An Information‐Theoretic Approach to Deciphering the Hippocampal Code.” Neural Information Processing Systems 5: 1030–1037.
40. Solstad, T. , Boccara C. N., Kropff E., Moser M.‐B., and Moser E. I.. 2008. “Representation of Geometric Borders in the Entorhinal Cortex.” Science (New York, N.Y.) 322, no. 5909: 1865–1868. doi.org/10.1126/science.1166466
41. Stachenfeld, K. L. , Botvinick M. M., and Gershman S. J.. 2017. “The Hippocampus as a Predictive Map.” Nature Neuroscience 20, no. 11: 1643–1653. doi.org/10.1038/nn.4650
42. Stewart, S. , Jeewajee A., Wills T. J., Burgess N., and Lever C.. 2014. “Boundary Coding in the Rat Subiculum.” Philosophical Transactions of the Royal Society, B: Biological Sciences 369, no. 1635: 20120514. doi.org/10.1098/rstb.2012.0514
43. Taube, J. , Muller R., and Ranck J.. 1990. “Head‐Direction Cells Recorded From the Postsubiculum in Freely Moving Rats. II. Effects of Environmental Manipulations.” Journal of Neuroscience 10, no. 2: 436–447. doi.org/10.1523/JNEUROSCI.10-02-00436.1990
44. Tsao, A. , Sugar J., Lu L., et al. 2018. “Integrating Time From Experience in the Lateral Entorhinal Cortex.” Nature 561, no. 7721: 57–62. doi.org/10.1038/s41586-018-0459-6
45. Ulanovsky, N. , and Moss C. F.. 2007. “Hippocampal Cellular and Network Activity in Freely Moving Echolocating Bats.” Nature Neuroscience 10, no. 2: 224–233. doi.org/10.1038/nn1829
46. Wang, F. , and Bicanski A.. 2026. “Dynamic Updating of Cognitive Maps via Traces of Experience in the Subiculum.” Hippocampus 2, no. 36: e70078. doi.org/10.1002/hipo.70078
47. Wells, C. E. , Amos D. P., Jeewajee A., et al. 2013. “Novelty and Anxiolytic Drugs Dissociate Two Components of Hippocampal Theta in Behaving Rats.” Journal of Neuroscience 33, no. 20: 8650–8667. doi.org/10.1523/JNEUROSCI.5040-12.2013
48. Whittington, J. C. , Muller T. H., Mark S., et al. 2020. “The Tolman‐Eichenbaum Machine: Unifying Space and Relational Memory Through Generalization in the Hippocampal Formation.” Cell 183, no. 5: 1249–1263. doi.org/10.1016/j.cell.2020.10.024

---

Republished from the open web under CC-BY. Authors: Newman EL, Mashanova-Galikova I, Tiganj Z, Lever C. Read the original.