American Association for the Advancement of Science
Scalable representation of time in the hippocampus
DOI 10.1126/sciadv.abd7013 , Volume: 7 , Issue: 6 , Pages: 0-0
Article Type: research-article, Article History
•
•
• Altmetric

### Notes

Abstract

Hippocampal “time cells” encode specific moments of temporally organized experiences that may support hippocampal functions for episodic memory. However, little is known about the reorganization of the temporal representation of time cells during changes in temporal structures of episodes. We investigated CA1 neuronal activity during temporal bisection tasks, in which the sets of time intervals to be discriminated were designed to be extended or contracted across the blocks of trials. Assemblies of neurons encoded elapsed time during the interval, and the representation was scaled when the set of interval times was varied. Theta phase precession and theta sequences of time cells were also scalable, and the fine temporal relationships were preserved between pairs in theta cycles. Moreover, theta sequences reflected the rats’ decisions on the basis of their time estimation. These findings demonstrate that scalable features of time cells may support the capability of flexible temporal representation for memory formation.

Shimbo, Izawa, and Fujisawa: Scalable representation of time in the hippocampus

## INTRODUCTION

It has been proposed that animals have a cognitive map that serves as an internal representation of the relationships between entities in the world and supports flexible behavior (1, 2). The hippocampus has been thought to provide cognitive maps for spatial navigation (3) because hippocampal cell assemblies organize map-like representations of environments. Hippocampal place cells use both rate and temporal coding for spatial representation. Rate coding takes the form of receptive fields of locations [i.e., place fields (4)]. On the other hand, spike timings of place cells are strongly modulated by concurrent theta oscillations (4 to 10 Hz), exhibiting spike phase shifts as a function of distance relative to the center of the receptive field [i.e., phase precession (57)]. As a consequence of this unique relationship of spiking phase and relative position, temporal coding can be seen as cell assembly sequences in theta cycles that compressively represent sequences of past, current, and future positions [i.e., theta sequence (810)]. With these two forms of coding, spatial representations of the hippocampus involve both positional and relational information. The spatial representation of place cell assemblies for a familiar environment is often scalable, and they can rescale their place fields when the spatial environment is resized (6, 1113). This indicates that relative positional relationships may provide more fundamental information for spatial recognition in the hippocampus than absolute positions (1416).

Recent research has suggested that the hippocampus also provides a cognitive map for nonspatial entities, such as temporal information (2, 14, 1719). During temporally organized experiences, hippocampal “time cells” fire at specific moments and form self-organized sequence activity that provides map-like representations of temporal structures of episodes (1719). The representation of temporal information is also supported in the forms of rate and temporal coding frameworks (17, 20). However, not many studies have focused on investigating the systematic reorganizations of temporal representations in the hippocampus.

To address this, we investigated the scalability of temporal representations of time cells in the hippocampus. To this end, we conducted temporal bisection tasks, which required rats to discriminate between long- and short-interval times to get a reward. The key feature of the tasks was that the sets of interval times to discriminate between were designed to be extended or contracted across the blocks of trials. In this task, we found that the majority of the time cells in CA1 region showed scalable representations of temporal information across the trial blocks. To identify the mechanism underlying the scalable representation of time, we investigated fine temporal structures of sequence activity of time cells associated with theta oscillations in the hippocampus. Theta sequences of time cells were also scalable, and the fine temporal relationships between cell pairs within single theta cycles were preserved. Moreover, while analyzing the data during test trials, we found that activities of time cells reflected the rats’ decisions on the basis of their temporal perception. Our results indicated that temporal representations in the hippocampus might also be scalable with respect to rate and temporal coding frameworks. These findings suggest that a common circuit mechanism may underlie map-like representations of spatial and nonspatial information in the hippocampus because of the similarity between place cells and time cells.

## RESULTS

We recorded multiple extracellular single units from hippocampal CA1 in the rats during the task behavior using high-density silicon probes (Materials and Methods). Units were categorized as putative pyramidal cells or interneurons on the basis of the shapes of their waveforms and firing rates (fig. S1) (21, 22). The total number of units recorded in the tasks is summarized in table S1.

### Scalable representation of time in the temporal bisection tasks

Rats were trained on the temporal bisection task. During this task, the sets of time intervals to discriminate between were designed to be scaled up and back across blocks of trials (task 1; Fig. 1A and fig. S1). In each trial, the rats were forced to run on a treadmill for long or short intervals of time. After the forced run, the rats were required to select the left or right arm in the Y-maze, which were associated with long or short time intervals, respectively. Water drops were provided as a reward at the end of the correct arm. A single session consisted of three blocks of trials. During the first block (block 1), the sets of interval times to discriminate between were 10 (long) and 5 s (short). During the second block (block 2), the sets of interval times were changed to 20 (long) and 10 s (short), i.e., extended two times compared with block 1. During the third block (block 3), the sets of intervals were returned to those used in block 1.

Fig. 1
(A) Schematic of the task (task 1). The rats ran on the treadmill for long or short intervals and then selected the left or right arm associated with water reward for long or short intervals, respectively. Sets of intervals were 10 s (long) and 5 s (short) in blocks 1 and 3 and 20 s (long) and 10 s (short) in block 2. (B, a to c) Firing patterns of three representative time cells in blocks 1 (left), 2 (center), and 3 (right). Top: Raster plots of long- and short-interval trials. Bottom: PETHs of long (solid)– and short (dotted line)–interval trials. (C) Firing patterns of time cells in blocks 1 (left), 2 (center), and 3 (right) (n = 454 units from seven rats). Each row represents PETHs of single neurons during long-interval trials in each block. The color scale represents the firing rate of each neuron. The neurons were ordered according to the peak time in block 1. (D) Population vector analyses of the data (C). Left: Autocorrelation of the population vector matrix of block 1. Center and right: Cross-correlations of blocks 1 and 2 and blocks 1 and 3, respectively. The color scale represents Pearson’s correlation coefficient. (E) Distributions of the scaling factors of time cells between blocks 1 and 2 (blue) and blocks 1 and 3 (green).Scalable time representation of CA1 neurons in the temporal bisection task.

To assess the firing activities of CA1 single pyramidal cells during the interval period, we compared peri-event time histograms (PETHs) of the discharge activities of single units in each block. First, we assessed the time-dependent activities of pyramidal cells in blocks 1 and 3 (Materials and Methods). In this study, we defined time cells as pyramidal cells that preferentially and robustly represented the information of elapsed time during the interval period in blocks 1 and 3 using the following criteria: First, to select units that represented elapsed-time information during the interval period, we selected units with maximum firing rates of >2 Hz and information rates (23) of >0.1 bit per spike during the intervals in blocks 1 and 3. Second, to select neurons that preferentially coded temporal information (rather than distance or spatial information), we fitted the spiking activity of each neuron to generalized linear models (GLMs) using interval time, running distance, or spatial position as covariates (19) (note that the speed of the treadmill was assigned as either 12 or 24 cm/s in each trial randomly and independently from the total duration of the trial, for the purpose of controlling for the influence of running distance). Then, we selected units whose spiking activity was best fitted to the model with elapsed time compared to the other models with running distance or spatial position (Materials and Methods). Third, to select neurons that robustly coded temporal information across the trial blocks, we selected units whose peak firing times of the PETHs in blocks 1 and 3 were within a 50% difference. Using these criteria, 22.2% of pyramidal cells (454 of 2041 units) were categorized as time cells (Fig. 1, B and C, and table S1).

Next, we assessed how activity patterns of the time cells changed when the set of interval times was extended two times in block 2. Population vector analysis of the time cells indicated that activity patterns of the neuronal population were extended in block 2 compared to block 1 (Fig. 1D; Materials and Methods). To investigate the scaling properties of the time cells, we quantified the scaling factors of the PETHs for each time cell. In this analysis, we made the firing rate models for each time cell. The model assumed that the cell’s firing activity as a function of elapsed time was scaled across the block; thus, the PETH in block 2 or 3 should be extended from that in block 1 with respect of time (Materials and Methods). We fitted the parameters of the scaling factor of the model using maximum likelihood estimation. With this method, the scaling factors of the time cells in block 2 were 1.81 ± 0.71 (means ± SD), while those in block 3 were 1.17 ± 0.52 (k = 0.65 and P = 3.1 × 10−86 for block 2 versus block 3, Kolmogorov-Smirnov test; Fig. 1E). This result indicated that a major portion of the time cells had scalable representations of elapsed-time information.

Then, we tested whether the time representation of CA1 pyramidal cell were also scalable to a shrinking direction. To this end, another set of rats were trained in a temporal bisection task with a different block paradigm (task 2; fig. S2A). This time, the sets of interval times to bisect were scaled down and back across the blocks. In blocks 1 and 3, the sets of time intervals to bisect were 20 (long) or 10 s (short). In block 2, the sets of time intervals were changed to 10 (long) or 5 s (short), i.e., the intervals were half of those used in block 1. With the same criteria above, we found that 18.9% of the pyramidal cells (241 of 1276 units) were classified as time cells (fig. S2, B and C, and table S1). We also investigated the scaling properties of the time cells using the same procedure in task 1. Population vector analysis of the time cells indicated that activity patterns of the neuronal population were shrunken in block 2 compared to block 1 (fig. S2D). The scaling factors of the time cells in block 2 were 0.61 ± 0.25 (means ± SD), while those in block 3 were 1.04 ± 0.27 (k = 0.76 and P = 3.3 × 10−62 for block 2 versus block 3, Kolmogorov-Smirnov test; fig. S2E). These results indicated that the time representation of CA1 pyramidal cells was scalable to both expanding and shrinking directions.

### Scalability of theta phase precession and theta sequences of time cells

During spatial navigation, hippocampal cell assemblies generate sequential structures in single theta cycles, termed theta sequences, which compressively represent trajectory sequences of past, current, and future positions (810). In this study, we investigated whether and how the time cell assemblies generate theta sequences of temporal information during the interval period.

First, we checked for the presence of theta phase precession as a function of elapsed time during the interval period in the temporal bisection task (task 1). The simultaneously recorded local field potential (LFP) was dominated by a robust theta frequency (4 to 10 Hz) oscillation during the interval period (17). We analyzed the phase-time relationships of single neurons and found that 77.3% of the time cells (351 of 454 units) displayed phase precession based on our criteria (Fig. 2A; Materials and Methods). To quantify the slope of the phase precession of each time cell in each block, we computed circular-linear fits of data samples (Materials and Methods) (24). Figure 2B plotted the relationships of the slopes of phase precession of each cell between blocks 1 and 2 and those between blocks 1 and 3. Group analysis revealed that the slope of phase precession in each cell in block 2 divided by that in block 1 was 0.63 ± 0.03 (means ± SEM), which was significantly smaller than that in block 3 divided by that in block 1 (0.96 ± 0.04) (paired t test, t350 = 6.8 and P = 2.2 × 10−11 ; Fig. 2C). Overall mean phases of those neurons were not significantly different across the blocks (177.4° ± 2.7° in block 1, 184.9° ± 2.7° in block 2, and 190.4° ± 2.8° in block 3, means ± SEM; F1,350 = 1.84 and P = 0.16, Watson-Williams test). This indicated that the phase-time relationship between each cell also reflected the scaling of temporal receptive fields across the blocks.

Fig. 2
(A, a to f) Firing rates and theta phases of three representative time cells in blocks 1 (left column), 2 (center column), and 3 (right column) in task 1. Top: Raster plots of short and long-interval trials. Middle: PETHs of long (solid line)– and short (dotted line)–interval trials. Bottom: Theta phase plots as a function of elapsed time. Each dot represents an action potential from the long (light color)– and short (dark color)–interval trials. The black line in each panel represents circular-linear fits of the data (see Materials and Methods). (B) Relationships of the slopes of phase precession between blocks 1 and 2 (blue dots) and between blocks 1 and 3 (green dots) in each unit in task 1 (n = 351 units). Solid lines represent the fit of the linear regression of these data. (C) Slope of phase precession in each cell in block 2 divided by that in block 1 (blue). Slope of phase precession in each cell in block 3 divided by that in block 1 (green). Means ± SEM. ***P < 0.001, paired t test.Scalability of theta phase precession of time cells.

Previous studies have reported that theta phase precession of single neurons does not always guarantee the concurrent presence of theta sequences of cell assemblies (25, 26). To address whether time cell assemblies formed theta sequences (810) during the task (task 1), we investigated the temporal structure of spike sequences of time cell assemblies in single theta cycles. First, we applied pairwise analyses of theta sequences to spiking activities of cells during the interval periods (9, 17). If theta sequences are present during spatial navigation, the distance of the place fields of a neuronal pair should be reflected in their distance of spiking phases in single theta cycles; this can be estimated by cross-correlation analysis during single trials (9, 22). Here, we applied this pairwise analysis to assess whether the difference in PETH peak times of a neuronal pair was correlated with differences in spiking theta phases (Fig. 3; Materials and Methods). In the analysis, we selected the cell pairs in which neurons were both activated and phase modulated in the same theta cycles during the trials (Materials and Methods). First, we checked the distance of the peak times of PETHs of each pair of time cells (Fig. 3, A and C), which can also be computed from the peak shift of the cross-correlation of the PETHs of the cell pair. Second, we calculated the cross-correlogram (CCG) of spike times of each pair, tested the theta modulation, and quantified the peak shift from point zero (Fig. 3, A and B). Figure 3D shows the relationship between CCG peak shifts and PETH peak time differences of all of the cell pairs that had significant theta modulations (Materials and Methods). Group analysis of the pairs revealed a significant positive correlation between the CCG peak shifts and peak time differences of the PETHs of the pairs in blocks 1 and 2 (n = 694 pairs; correlation coefficient r = 0.42 and P = 1.2 × 10−31 in block 1 and r = 0.31 and P = 6.8 × 10−17 in block 2; Materials and Methods and Fig. 3D), suggesting that theta sequences were present during the time intervals in the temporal bisection task. We fitted a linear regression model to these data and found that the slope of the regression line in block 2 was about half of that in block 1 (red lines in Fig. 3D; block 1, 29.2; block 2, 15.0 ms/s). These data indicated that the compression ratio of the temporal sequence information in single theta cycles in block 2 was about two times higher than that in block 1. To assess the differences in the CCG of peak shifts between cell pairs in blocks 1 and 2, we also applied a linear regression analysis to the data of the relationship between the CCG of peak shifts from blocks 1 and 2 (Fig. 3E) and found that the slope of the regression line was 0.98 (red line in Fig. 3E). Statistical analysis revealed that there were no significant differences in the CCG of peak shifts of cell pairs between blocks 1 and 2 (Z = 1.13 and P = 0.26, signed-rank test), demonstrating that temporal relationships of cell pairs in single theta cycles were preserved across the blocks.

Fig. 3
(A, a and b) Left: Firing rates and theta phases of two representative pairs of time cells (red and blue) in blocks 1 and 2 (task 1). The arrows indicate the peak firing rates. Right: Cross-correlograms (CCGs) of the pairs in blocks 1 (orange) and 2 (purple). The arrows indicate the peaks of the CCGs. Note that the CCG peak time shifts in (Ab) were larger than those in (Aa) (right), reflecting larger distances in PETH peak times of the pairs in (Ab) than those in (Aa) (left). a.u., arbitrary units. (B) Cross-correlations of PETHs of the pairs of time cells. (C) CCGs of the pairs of time cells in blocks 1 and 2 (n = 694 pairs). The orders of the pairs are the same in (B) and (C). (D) Pairwise analysis of the temporal compression of spike sequences in single theta cycles in blocks 1 and 2. Each dot represents a single neuronal pair. Horizontal axes show peak distances of PETHs, and vertical axes show peak shifts of the CCGs. The slopes of the regression line (red) were 29.2 and 15.0 ms/s in blocks 1 and 2, respectively. (E) Relationship between the time shifts of the CCGs from blocks 1 and 2. The slope of the regression line (red) was 0.98.Pairwise analysis of the temporal compression of spike sequences in single theta cycles.

We also investigated the theta sequences of time cell assemblies using a Bayesian decoding framework (22, 25, 27). In spatial navigation, rats’ positions can be decoded from the spiking activity of their place cells (2729). Using this method, previous studies have demonstrated that the spiking activity of place cell assemblies represents trajectory sequences of past, current, and future positions in single theta cycles (25, 26). We applied this method to the spiking activity of time cell assemblies to investigate neuronal representations in single theta cycles. First, we divided the spike trains of simultaneously recorded time cells into 20-ms time bins in each trial. Second, we estimated the probability density of temporal information from the spike counts of the cells in each time bin using the Bayesian decoding method (Materials and Methods). Figure 4A represents the decoded elapsed-time information at each time point in a single example trial. Then, we aligned the x axis of the decoded temporal information with concurrent theta phases and the y axis to the time that the spikes were sampled. The probability densities of the decoded time at each theta phase demonstrated robust theta sequence formation of the time cells (Fig. 4B). On descending phases of theta (0° to 180°), the decoded probability densities of time reflected past time relative to the time when the spikes were sampled. On the other hand, on ascending phases of theta (180° to 360°), decoded probability densities of time reflected the future time relative to the time when the spikes were sampled. The theta sequence of block 2 was stretched in the vertical direction compared to that of block 1 (Fig. 4B), indicating that the theta sequence of block 2 represented longer time periods in single theta cycles than that of block 1. To quantify this, we assessed scaling factors of the theta sequences in the direction of time (y axis) between blocks 1 and 2 (Fig. 4C; Materials and Methods). Group analysis revealed that the scaling factor of the theta sequences between blocks 1 and 2 were 1.62, whereas that between blocks 1 and 3 was 1.13 (paired t test, t25 = 3.11 and P = 0.0031). This indicated that the temporal representation of theta sequence in single theta cycles was significantly more extended in block 2 than in block 1.

Fig. 4
(A) Decoded elapsed-time information at each time point in a single example trial in task 1. Top: Simultaneously recorded LFPs (wide band and 4 to 10 Hz) and units from 10 pyramidal neurons ordered by the peak times of their firing rates from bottom to top. Bottom: Probability density of temporal information was computed from simultaneously recorded pyramidal cell activity in the trial. Bayesian decoding was performed at every 1-ms time step with a 20-ms time window. (B) Probability densities of the decoded time at each theta phase in blocks 1, 2, and 3 in the example single session in task 1. Horizontal axes represent theta phases in which spiking data were sampled, and vertical axes represent time decoded from spiking activity relative to the sampled time. The result indicates that decoded relative time was extended vertically in block 2. (C) Scaling factors of theta sequences between blocks 1 and 2 and blocks 1 and 3 (n = 26 sessions; sessions that contained more than eight time cells were selected). Means ± SEM. **P < 0.01, paired t test.Scalability of theta sequences of time cells.

Last, we assessed whether the theta sequences of the time cells were merely the consequence of phase precession of individual neurons or had higher correlative structures than those predicted from phase precession (9, 10, 30). To address this, we applied the phase resampling method to the single theta sequence events of the time cells (fig. S4) (10). In each theta sequence event, we constructed surrogate datasets while preserving the distributions of phase precession of individual neurons and compared the correlative structures of the spike sequences between the original and surrogate events (fig. S4A). We found that 12.6% of theta sequence events in block 1 and 10.5% of events in block 2 had significantly higher correlations (P < 0.05, permutation test for each event; fig. S4A) than those of surrogate events, and these ratios were significantly higher than chance level as the group data (Clopper-Pearson CIs for 99%, CI = 12.0 to 13.1% in block 1 and CI = 10.1 to 10.9% in block 2; chance level is 5%; fig. S4B). These results showed that the organizations of theta sequences of time cells have higher correlative structures than those predicted from phase precession and that these highly correlative structures are preserved across the trial blocks, indicating that the theta sequences of the time cells were not merely the consequence of phase precession of individual neurons. Together, these results demonstrated that the theta sequences of time cells were also scalable across the trial blocks, and the fine temporal relationships between cell pairs were conserved within single theta cycles.

### Activities of time cell assembly reflect animals’ decisions

Last, we asked how the assembly dynamics of time cells reflect recognition of elapsed time in the rats. To this end, we examined the temporal bisection task with an additional condition (task 4; Fig. 5A). Along with the normal trials of short and long time intervals, we randomly inserted test trials using the intervals that were the geometric means of short and long ones. During the test trials, the rats were also required to select the left or right arm of the Y-maze after forced running during the interval period; however, a reward was omitted in both arms. Because the left and right arms were associated with long- or short-interval times in the normal trials, respectively, we hypothesized that the rats’ choices in the test trials would reflect their judgment about short or long time intervals depending on their estimation of the time. The task sessions also consisted of three blocks. In block 1, the sets of interval times were 10 (long), 5 (short), and 7.07 s (test). In block 2, the sets of interval times were 20 (long), 10 (short), and 14.14 s (test). In block 3, the sets of interval times were 10 (long) and 5 s (short).

Fig. 5
(A) Schematic of the task (task 4). The task consisted of normal trials with long and short intervals (rewarded) and test trials with geometric means of short and long intervals (not rewarded). The sets of the intervals were 10 s (long), 5 s (short), and 7.07 s (test) in block 1; 20 s (long), 10 s (short), and 14.14 s (test) in block 2; and 10 s (long) and 5 s (short) in block 3. (B) Probability densities of decoded time at each theta phase in selected-long and selected-short test trials in block 1. Horizontal axes represent theta phases from which spiking data were sampled, and vertical axes represent time decoded from spiking activity. Spiking activity was sampled for 1 s (between 6 s and 7 s, indicated by dotted lines) of the interval of test trials. (C) Decoded time in selected-long and selected-short test trials in every 1-s duration of the interval in block 1. (D) Probability densities of decoded time at each theta phase in selected-long and selected-short test trials in block 2. Spiking activity was sampled for 2 s (between 12 s and 14 s). (E) Decoded time in selected-long and selected-short test trials in every 2-s duration in block 2. Means ± SEM. *P < 0.05 and **P < 0.01. Two-way analysis of variance (ANOVA) with Bonferroni posttest for (C) and (D). n = 13 sessions for (B) to (E).Activities of time cell assembly reflect animal’s decisions.

We investigated whether activities of the time cells reflected the choices of the rats in the test trials. First, we estimated the PETHs of the time cells during the normal trials in each block and used this information for decoding templates using the same procedures of the previous tasks. Using a Bayesian decoding technique (Materials and Methods), we estimated the differences in the theta sequences of cell assemblies between the test trials in which the rats selected the left (“selected long”) or the right arm (“selected short”) of the Y-maze. The theta sequences in the selected-long and selected-short trials (Fig. 5B) revealed that temporal representation decoded from spiking activity in the selected-long trials was later than that in the selected-short trials. Then, we quantified differences in decoded time between these two trials. We estimated theta sequences of the selected-long and selected-short trials from the spiking activity during every 1-s duration in the interval period in block 1 in each session. Then, we computed the peak times of the averaged theta sequences of each 1-s duration to clarify in which period the differences arose. There were significant differences in decoded temporal information between the selected-long and selected-short trials [n = 13 sessions from three rats; two-way analysis of variance (ANOVA) followed by Bonferroni posttest; interaction of trial type and period, F6,72 = 5.39 and P = 0.0001; simple effects for interaction of trial type and period, simple effects of trial type at 0 to 1 s, F1,12 = 0 and P = 1; at 1 to 2 s, F1,12 = 1.87 and P = 0.20; at 2 to 3 s, F1,12 = 1.37 and P = 0.26; at 3 to 4 s, F1,12 = 9.17 and P = 0.011; at 4 to 5 s, F1,12 = 1.99 and P = 0.18; at 5 to 6 s, F1,12 = 4.42 and P = 0.057; and at 6 to 7 s, F1,12 = 9.33 and P = 0.01; Fig. 5C]. We also quantified differences in decoded time between these two trials in every 2-s duration in block 2 and found the same tendency (two-way ANOVA followed with Bonferroni posttest; interaction of trial type and period, F6,72 = 3.05 and P = 0.0103; simple effects for interaction of trial type and period, simple effects of trial type at 0 to 2 s, F1,12 = 0.43 and P = 0.53; at 2 to 4 s, F1,12 = 1.18 and P = 0.30; at 4 to 6 s, F1,12 = 2.28 and P = 0.16; at 6 to 8 s, F1,12 = 6.52 and P = 0.025; at 8 to 10 s, F1,12 = 6.86 and P = 0.022; at 10 to 12 s, F1,12 = 7.31 and P = 0.019; and at 12 to 14 s, F1,12 = 9.35 and P = 0.0099; Fig. 5, D and E). These results indicated that the activity of time cell assemblies actually reflected animals’ decisions during the test trials.

## DISCUSSION

In this study, we investigated the scalability of time cells in the hippocampus with respect to rate and temporal coding. We monitored CA1 neuronal activity during the temporal bisection tasks and found that temporal representations of pyramidal neurons were scalable when the sets of time intervals were varied (Fig. 1 and fig. S2). We tested the task dependency of temporal representation in the hippocampus and found that the proportion of pyramidal cells that were time cells was significantly reduced in the light discrimination task, during which the rats were not required to estimate time, as compared with that in the temporal bisection task; however, the time cells displayed scalable representations of elapsed time when the sets of time intervals were varied (fig. S3). We also investigated fine temporal structures of sequential activity of time cells associated with theta oscillations in the hippocampus. Theta sequences of time cells were also scalable, and the fine temporal relationships between cell pairs in single theta cycles were conserved (Figs. 2 to 4). Moreover, we investigated the activity of the time cells during test trials and found that decoded temporal information from the theta sequences of the time cells was correlated to the choices of rats (Fig. 5), indicating that the activity of time cells reflected rats’ decisions on the basis of their estimation of the time. These results demonstrated that temporal representation in the hippocampus is scalable with respect to rate and temporal coding (Fig. 6).

Fig. 6
Representations of time information were scalable in rate and temporal coding in the hippocampus. Regarding fine temporal structures of spiking activity in single theta cycles, the temporal relationships of units are preserved between short and long episodes. Regarding neuronal representations of elapsed-time information in single theta cycles, the temporal representations of spike sequences are scaled between short and long episodes.Scalability of rate and temporal coding of time cells.

### Reorganization of neuronal representations of space and time

Previous research on spatial navigation has reported reorganizations of spatial representations of CA1 pyramidal neurons following changes in the environment. When the boundaries of a familiar environment are expanded or contracted, place cell assemblies often rescale their place fields to fit the new configuration of the environment, enabling the scalable representation of the environment (6, 11, 12). On the other hand, when animals experience two different environments, the neuronal representations for the environments are completely remapped (31). When one environment is gradually morphed into a different environment, the neuronal representations are progressively transformed from one environment to the other (32). These observations suggest that hippocampal cell assemblies organize separate map representation for each environment, which can be scalable as long as the animals recognize that the environment is the same (33, 34).

Representations and reorganizations of temporal information of CA1 pyramidal neurons have also been the focus of several recent studies. Pastalkova and colleagues (17) found that reliably and continually changing CA1 cell assemblies appeared during wheel running. Their results indicated that CA1 pyramidal cells can encode not only spatial but also temporal information. MacDonald et al. (18) also reported that CA1 neurons represented temporal information during the interval period in the object-odor association task. To clarify whether CA1 time cells represent information of elapsed time or distance during running on the treadmill, Kraus et al. (19) compared firing patterns of neurons during periods of fixed time intervals with different running speeds and those during variable time intervals with fixed running speeds. They found that a group of neurons represented temporal information even if the running speed was varied and that another group of cells represented running distance information rather than temporal information. In the current study, we investigated the reorganizations of temporal representations of CA1 pyramidal cells when the interval periods were varied. We conducted the temporal bisection task that required rats to measure the interval time whose duration was associated with the rewarded side. The strength of using the temporal bisection task is that we can scale the interval times while preserving the association between the duration and the rewarded side. The design of the task enabled us to investigate the reorganization of temporal representations of the cells systematically. Here, we demonstrated that temporal representations of CA1 neurons were scalable when the temporal structure of familiar episodes was systematically enlarged or shortened (Fig. 1 and fig. S2). Theoretical studies have proposed several mathematical frameworks that enable scalable representations of both spatial and temporal information (15, 3538).

### Scalability of phase precession and theta sequences of time cells

To clarify the mechanisms underlying scalable representations of time, we investigated fine temporal structures of sequence activities of time cells associated with theta oscillations in the hippocampus. During spatial navigation, spike timings of place cells exhibit spike phase shifts as a function of distance relative to the center of the receptive field [i.e., phase precession (57)]. As a consequence of this unique relationship of spiking phase and relative position, hippocampal place cell assemblies generate sequential structures across single theta cycles (i.e., theta sequences) that compressively represent sequences of past, current, and future positions (810, 39). Theta phase precession is the form of temporal coding that appears at single-cell levels (40), and theta sequences are the form of temporal coding that appears at cell ensemble levels (9, 10).

Similar to the place fields, when the boundaries of a familiar environment are expanded or contracted, theta sequences also rescale to fit the new configuration of the environment (13, 41). The fine temporal relationships between the pairs of place cells are preserved even after the theta sequences are rescaled (13). Similar to previous research, the results from this study revealed that phase precession and theta sequences of time cells were present during the interval period (17, 20). We found that phase precession and theta sequences of time cells were also scaled when the sets of time intervals were changed, and the fine temporal relationships between the pairs of the time cells were preserved (Figs. 2 to 4 and 6). These results indicated that a common mechanism underlies scalable representations for both spatial and temporal information. The scalable feature of time cells may be advantageous for a flexible representation of time when future plans change.

### Time cells and time perception

Previous research reported that several brain regions, such as prefrontal cortex, striatum, midbrain, entorhinal cortex, and hippocampus, are involved in perception and representation of elapsed time (4248). In the hippocampus, information of elapsed time is represented by the sequential activity of time cells, which may be suitable for precise representations of time on a scale of tens of seconds (1719). Lesion studies have also supported the involvement of the hippocampus in time estimation during the interval period (49, 50). However, the involvement of sequential activity in time perception is still unknown. In this study, we investigated the roles of time cells during the recognition of time. First, our data demonstrated that the temporal representations of theta sequences of time cells were significantly different in the selected-long and selected-short trials in the test trials (Fig. 5). These results indicated that temporal representations of theta sequences with time cells might have actually reflected the rats’ recognition of elapsed-time information. Second, we demonstrated that the organization of time cells was dependent on task demands. The proportion of pyramidal cells that were time cells was significantly higher in the temporal bisection task than that in the light discrimination task, although both of the tasks contained an interval period before decision-making (fig. S3D). A previous study using the object-odor association task showed that the representations of the majority of time cells were randomly remapped when the interval duration was extended (18). The difference in the task demands could also explain the differences in the time cell reorganization in our study. These observations indicate that the formation of receptive fields for temporal information may be promoted when the task requires animals to measure elapsed time. We hypothesize that CA1 pyramidal cells can organize cognitive maps of either spatial or temporal information depending on the requirements of the environment (16, 37).

In conclusion, the results of this study demonstrated that hippocampal pyramidal cells can organize scalable representations of time in both rate and temporal coding frameworks, suggesting that a common circuit mechanism may underlie map-like representations of spatial and nonspatial information in the hippocampus (16, 37). Such scalable, map-like representations of spatial and temporal information may allow for the formation of flexible memory to adapt to complex environments.

## MATERIALS AND METHODS

### General behavioral methods

All experimental protocols were approved by the RIKEN Institutional Animal Care and Use Committee. Six-week-old male Long-Evans rats were purchased from Japan SLC Inc. (Shizuoka, Japan). Rats were housed in a temperature-controlled room (20° to 22°C) under a light-dark cycle (lights on from 8:00 to 20:00). When the rats weighed more than 250 g, they were subjected to a water deprivation schedule and began training on the temporal bisection task (8 to 10 weeks old at the start of training). Under the water deprivation schedule, the amount of water intake was restricted to 10 ml outside of the task, but the food was available ad libitum.

All behavioral experiments were conducted in a soundproof chamber (W120 cm by D120 cm by H200 cm; S1212, Star lite, Fukui, Japan) located in an experimental room. The apparatus was installed on the center of the chamber (60 cm above the room floor). The apparatus consisted of a treadmill (W15 cm by D40 cm; SVKR type, Misumi Group Inc., Tokyo, Japan) and Y-maze (20 cm length for the central arm and 24 cm length for the side arms), which were partitioned by a transparent acrylic automated door (W40 cm by H40 cm; fig. S1B). The treadmill was surrounded by transparent acrylic walls (50 cm in height) except for the door side (right side to the running direction). Two water ports (D2 cm by H1.5 cm) for reward were placed at 5 cm from the ends of the arms of the Y-maze, and one port (L-shape; 2 cm by 2 cm by 1.5 cm) was placed on the front walls of the treadmill. These water ports were connected to micropumps (Type 7615; Christian Bürkert GmbH & Co. KG, Ingelfingen, Germany) that delivered 0.2% saccharin water (M5N1991; Nacalai tesque, Kyoto, Japan) as a reward. Two light-emitting diode (LED) lights (OXL/CLH/80 Series; Oxley, Singapore) were positioned 10 cm above and 5 cm outside the end of the arm. For system control, four infrared sensors were set in the apparatus to detect the rats’ positions. One infrared sensor (PZ-G51P; Keyence, Osaka, Japan) was set on the treadmill, one was located at the entrance of Y-maze, and two were set at a 12-cm distance from the end of the arm. A charge-coupled device (CCD) camera (Q2F-00021; Microsoft, Redmond, WA) was placed on the ceiling to track the positions of the rats, which was used for analyzing data from the electrophysiology experiments.

Behavioral experiments were automatically controlled by the custom-written LabVIEW (National Instruments, Austin, TX) programs running on a Windows personal computer (PC) located outside of the chamber. The treadmill was driven with a stepper motor and the program controlled its speed. Images from the CCD camera and electrophysiological data were recorded by the dedicated recording software (Amplirec; Amplipex Ltd., Szeged, Hungary) that was running on another Windows PC outside of the chamber.

The rats (rat ids: s15, s18, s20, s23, s26, s32, and s39) were trained on a modified version of the temporal discrimination task (51). In this task, the rats were required to discriminate between long and short intervals on the treadmill. The experimental procedures of a single trial are as follows: When the rats came onto the treadmill from the Y-maze area, the door was closed, and a small water reward (15 μl) was provided at the water port on the wall of the treadmill (fig. S1A). After the sound stimulus (8-kHz sine wave, 1 s) was presented, the treadmill began to move, and the rats were forced to run for long or short intervals of time. The durations of long or short intervals were different across the blocks of trials (described below). The speed of the treadmill was randomly assigned to either 12, 18, or 24 cm/s in each trial (18 cm/s was used only for training periods before surgery.) After the forced run during the interval periods, the door opened and two LED lights were turned on. The rats were required to select either the left or right arm in the Y-maze, which were associated with long or short time intervals, respectively. If the rats selected the correct side, the LED lights were turned off, the sound stimulus (4-kHz sine wave, 1 s) was presented, and a water reward (20 μl) was provided at the water port of the correct arm. If the rats selected the incorrect side, the LED lights were turned off and the sound stimulus (11-kHz square wave, 1 s) was presented; however, the rats did not receive the water reward. After an incorrect response, the correction trials, which consisted of the same sets of intervals and speeds of the incorrect trials, were performed until the rats selected the correct side. One session consisted of 60 trials or a duration of 60 min.

Training sessions started after pretraining that included habituation to apparatus, running on the treadmill, and learning to drink water reward from water ports. A training session consisted of one block of trials, and the sets of times intervals to discriminate between were 10 (long) and 5 s (short). When the correct response rate was >80% for three successive sessions (not including the correction trials), the surgery was conducted. Following recovery from surgery (4 to 7 days), the water deprivation schedule and training session started again. After the animals demonstrated stable behavioral performances (>80% correct response rate for two consecutive sessions), a recording session started.

A recording session consisted of three blocks of trials with the standard-extended-standard block paradigm (Figs. 1 to 4). Similar to the training session, the sets of time intervals to discriminate between during the first block (block 1) was 10 (long) and 5 s (short). During the second block (block 2), the sets of interval times were changed to 20 (long) and 10 s (short) (i.e., extended twice from block 1). During the third block (block 3), the sets of intervals returned to those that were used during block 1. Each block consisted of 60 trials. For the sake of providing information about the block change to the animals, 1-min intervals were inserted before the start of the first trials of new blocks in the experiments. The speed of the treadmill was either set to 12 or 24 cm/s. Four sets of intervals and speeds in the trials (2 intervals × 2 speeds) were presented pseudorandomly in each block. Neural activity was analyzed when rats reached a correct response rate of >80% during both the long and short trials in each of the three blocks. The performance of the rats in the task is summarized in fig. S1C.

The rats (rat ids: s25, s28, and s33) were trained to discriminate between long and short time intervals on the treadmill using the standard-contracted-standard block paradigm (fig. S2). The procedure of a single trial was the same as task 1, except for the length of intervals. During the training sessions, the rats were trained to discriminate between 20- (long) and 10-s (short) intervals. During the recording sessions, the sets of interval times to bisect were 20 and 10 s in blocks 1 and 3. During block 2, the sets of time intervals were changed to 10 (long) and 5 s (short) (i.e., half of blocks 1 and 3).

The rats (rat ids: s30, s34, and s36) were trained on light discrimination task using the standard-extended-standard block paradigm (fig. S3). The procedure of a single trial was the same as task 1, except for the reward contingency. The rats were forced to run on a treadmill during long or short time intervals (5 or 10 s), but the intervals were not relevant to the rewarded arm. After the time intervals, the door opened, and the LED light of either the left or the right arm of the Y-maze turned on, and the rats were required to select the arm of the Y-maze according to the side of the lighted LED. The block paradigm was the same as that used in task 1. During blocks 1 and 3, the sets of interval times to bisect were 10 (long) and 5 s (short). During block 2, the sets of interval times were 20 (long) and 10 s (short).

The three rats (rat ids: s26, s32, and s39) performed the temporal bisection task with test trials (Fig. 5). The task procedure was the same as task 1, except for the random insertion of test trials. During the test trials, we used the geometric means of the short and long intervals as the time intervals. This is because previous studies reported that bisection points are located around the geometric means of intervals (51). The rats were also required to select the left or right arm of the Y-maze after the forced run, but rewards and sound stimuli were not provided in either of the arms. The block paradigm was similar to task 1. During block 1, the sets of interval times were 10 (long), 5 (short), and 7.07 s (test). During block 2, the sets of interval times were 20 (long), 10 (short), and 14.14 s (test). During block 3, the sets of interval times were 10 (long) and 5 s (short).

### Surgery and recording

After the training periods (15 to 46 weeks), we implanted silicon probes in 13 rats (rat ids: s15, s18, s20, s23, s25, s26, s28, s30, s32, s33, s34, s36, and s39) for the chronic recording of neuronal activities during the tasks. General surgical procedures for chronic recordings have been described previously (22). In this study, we used Buzsáki64sp silicon probes (NeuroNexus, Ann Arbor, MI), which consisted of six shanks (200-μm shank separation). Each shank had 10 recording sites (160 μm2 each site; ~1-megohm impedance), staggered to provide a two-dimensional arrangement (20-μm vertical separation; see fig. S1E inset). The rats were implanted with one or two silicon probes in the CA1 region of the hippocampus [rat s15, one probe located at anterior-posterior (AP) = −3.8 mm and medial-lateral (ML) = −2.8 mm; other rats, two probes located at AP = −3.8 m and ML = ±2.8]. The shanks were aligned parallel to the septotemporal axis of the hippocampus (i.e., 45° parasagittal). The silicon probe was attached to a micromanipulator and moved gradually to its desired depth position. During the recording sessions, the wide-band neurophysiological signals were acquired continuously at 20 kHz on a 256-channel Amplipex system (52) (KJE-1001, Amplipex Ltd., Hungary). The wide-band signal was downsampled to 1.25 kHz and used as the LFP signal. Spike sorting was performed semi-automatically using KlustaKwik2 (53) followed by the manual adjustment of the clusters. The tips of the probes moved spontaneously or manually by the experimenter between sessions. However, we could not exclude the possibility that some neurons that were recorded in different sessions were identical because spikes from sessions recorded on different days were clustered separately. Thus, the total number of independent neurons may have been overestimated (54).

### Data analysis

All analyses after spike sorting were performed using custom-written tools in MATLAB with Signal Processing, Optimization, and Statistic toolboxes (MathWorks, Natick, MA). The Circular Statistics Toolbox (55) was used for circular analyses.

### Selection of time cells

In this study, we defined time cells as the units that preferentially and robustly represented the information of elapsed time during the interval period in blocks 1 and 3 (Fig. 1C). Here, we used the following selection criteria. First, to select units that represented time elapsed information during the interval period, we selected the units whose maximum firing rates were >2 Hz, and information rates (23) were >0.1 bit per spike during the time intervals in blocks 1 and 3. Second, to select neurons that preferentially coded time information rather than distance or spatial information, we fitted the spiking activity of each neuron to GLM using elapsed time, running distance, or spatial position as the covariates (19). We selected units whose spiking activity was best fitted to the model that used elapsed time as the covariate as compared with the other models that used running distance and spatial position as the covariates. The details of the GLM analysis are described below. Third, to select neurons that robustly coded time information, we selected neurons whose peak firing times in blocks 1 and 3 were within a 50% difference. The numbers of selected cells that met these criteria are summarized in table S1.

### GLM analysis

To quantify the effects of elapsed time, running distance, and spatial position on neural activity during the interval period, we applied a GLM to the spiking activity of each neuron (18, 19). We modeled the spiking activity during the interval period as an inhomogeneous Poisson process with the firing rate as a function of covariates (elapsed time, distance, or position) that modulates spiking activity (18, 19). We compared three models. The “time” model was a fifth-order polynomial of elapsed time τ(t) described as

${\mathrm{\lambda }}_{\text{time}}\left(t\right)=\text{exp}\left({\sum }_{i=1}^{5}{\mathrm{\alpha }}_{i}\mathrm{\tau }{\left(t\right)}^{i}\right)$
where αi are parameters to fit. The “distance” model was a fifth-order polynomial of running distance on the treadmill (d(t)) described as
${\mathrm{\lambda }}_{\text{distance}}\left(t\right)=\text{exp}\left({\sum }_{i=1}^{5}{\mathrm{\beta }}_{i}d{\left(t\right)}^{i}\right)$
where βi are parameters to fit. The “position” model was a Gaussian-shaped place field described as
${\mathrm{\lambda }}_{\text{position}}\left(t\right)=\text{exp}\left({\mathrm{\gamma }}_{1}x\left(t\right)+{\mathrm{\gamma }}_{2}x{\left(t\right)}^{2}+{\mathrm{\gamma }}_{3}y\left(t\right)+{\mathrm{\gamma }}_{4}y{\left(t\right)}^{2}+{\mathrm{\gamma }}_{5}x\left(t\right)y\left(t\right)\right)$
where γi are parameters to fit. Then, the spiking activity of each neuron was fitted to each model. Here, the interval durations in block 2 were divided into 20-ms bins, and spike counts of the cells in each bin were estimated from real data. We used the glmfit function of MATLAB for GLM to fit the spiking data to each model and compared the deviances of the fit across the models. We selected the best model by comparing the deviances of the fit.

### Estimations of the scaling factors of time cells

To investigate the scaling properties of time cells, we quantified the scaling factors of PETHs for each time cell. In this analysis, we made the firing rate models for each time cell in blocks 2 and 3. The model assumed that the cell’s firing activity as a function of elapsed time was scaled across the block; thus, the PETH in block 2 or 3 should be extended from that in block 1 with respect of time.

The detailed procedure of the computation is as follows. First, we made firing rate models for each cell type in block 2 or 3, using the PETH of the cell in block 1. In this analysis, we defined the firing rate model f(t) in block 2 or 3 as

$f\left(t\right)={a}_{1}·{\text{PSTH}}_{\text{block}1}\left({a}_{2}·t\right)$
where a1 and a2 are the parameters for fitting, and t is elapsed time in each trial. PETHblock1(x) is defined as PETH of the cell in block 1 computed from real data when x is less than the duration of the longer trials (i.e., 10 s in tasks 1 and 3 and 20 s in task 2) and is defined as the minimum firing rate of the PETH when x is more than the duration of the longer trials. a2 represents the scaling factor of time cells. The setting range for a2 was from 0.1 to 3.5 for tasks 1 and 3 and from 0.1 to 2.0 for task 2.

Second, the interval durations in block 1 were divided into 20-ms bins, and spike counts of the cells in each bin were estimated from real data in block 1. Here, we hypothesized that the spike number of the cell in each bin would follow the Poisson distribution

$P\left(n\mid t\right)=\frac{{\left(\mathrm{\tau }f\left(t\right)\right)}^{n}}{n!}\text{exp}\left(-\mathrm{\tau }f\left(t\right)\right)$
where n is the spike number of each 20-ms time bin, t is elapsed time in each trial, and f(t) is the firing rate model described above. The likelihood function was defined as
$L\left(a\right)=\prod _{i=1}^{N}P\left(n\mid a\right)$
where N is the total number of time bins, and a is the parameter(s) in each model. We fitted the parameters of the model using maximum likelihood estimation.

### Population vector analysis

First, the normalized PETH for each cell was computed in each block. Then, the population vector matrix of all of the time cells (i.e., a single row vector represented a normalized PETH in a single neuron, and a column vector represented a population vector of all neurons at a single time bin) was constructed in each block. Last, the autocorrelation or cross-correlation of the population vector matrix in each block or between the blocks was computed (Figs. 1D and 2D). Here, the element at points i and j represented Pearson’s correlation coefficient between the population vectors computed at the i-th and j -th time bin in the trials (12).

### Theta phase precession

For the theta phase extraction, LFPs at the center of the CA1 pyramidal cell layers (where the amplitudes of the ripples were maximum) were filtered with a Butterworth filter with a pass-band range of 4 to 10 Hz. Instantaneous theta phases were estimated using the Hilbert transformation of the filtered signals. The theta phases estimated from the LFP at the center of the CA1 pyramidal cell layers were used as reference theta phases for all neurons.

The relationship between elapsed time and the theta phase of spiking activity in each unit during the interval period was analyzed. Theta phase precession was characterized by a significant linear-circular correlation (P < 0.01) between time and theta phases during the interval period (55).

To quantify the slope of phase precession in each cell, we applied a circular-linear fit to the data (24). Spike samples with elapsed time (ti) and theta phase (ϕi) for each neuron in each block (i = 1…n) were fitted to a linear model

$\mathrm{\varphi }\left(t\right)=2\mathrm{\pi }·a·t+{\mathrm{\varphi }}_{0}$
where the slope (a) and the phase offset (ϕ0) are the parameters for fitting. To obtain an optimized slope (a), we maximized the mean resultant length [R(a)] where
$R\left(a\right)=\sqrt{{\left[\frac{1}{n}\sum _{i=1}^{n}\text{cos}\left({\mathrm{\varphi }}_{i}-2\mathrm{\pi }·a·{t}_{i}\right)\right]}^{2}+{\left[\frac{1}{n}\sum _{i=1}^{n}\text{sin}\left({\mathrm{\varphi }}_{i}-2\mathrm{\pi }·a·{t}_{i}\right)\right]}^{2}}$

The phase offset ϕ0 was calculated as

${\mathrm{\varphi }}_{0}\left(a\right)=\text{arctan}*\frac{{\sum }_{i=1}^{n}\text{sin}\left({\mathrm{\varphi }}_{i}-2\mathrm{\pi }·a·{t}_{i}\right)}{{\sum }_{i=1}^{n}\text{cos}\left({\mathrm{\varphi }}_{i}-2\mathrm{\pi }·a·{t}_{i}\right)}$
where arctan* is the quadrant-specific inverse of the tangent (24).

### Pairwise analysis of theta sequences

Pairwise analysis of theta sequences was applied to the cell pairs in the interval period (Fig. 3) (22). In this analysis, we selected the cell pairs in which neurons were both activated and phase modulated in the same theta cycles during the trials. To this end, we used the following criteria: (i) Both neurons had phase precessions, and (ii) the CCG of the cell pairs had significant theta modulation. For the second criterion, we first computed the Fourier power of the CCG of the cell pair. Then, we made the surrogate dataset by shuffling the trials (i.e., shift predictor; 100 surrogate datasets for each pair) and computed the Fourier power of their CCG. If the Fourier power of the CCG of the real data was significantly larger than that of the surrogate dataset (P < 0.05) in the theta frequency band, then we considered that the CCG of the cell pairs had significant theta modulation.

### Reconstruction of time information from neuronal activities

A memoryless Bayesian decoding algorithm was used to estimate the information of elapsed time during the interval period from the spike trains (Fig. 4) (22, 25, 27, 28, 56). On the basis of the Bayes’ theory, the posterior probability of time from a trial onset (time) given spike trains from single neurons (spikes) was estimated as

$P\left(\text{time}\mid \text{spikes}\right)=\frac{P\left(\text{spikes}\mid \text{time}\right)·P\left(\text{time}\right)}{P\left(\text{spikes}\right)}$

The prior probability was estimated under the assumption of Poisson firing statistics and independent rates as

$\begin{array}{ccc}\hfill P\left(\text{spikes}\mid \text{time}\right)& & =\prod _{i=1}^{N}P\left({n}_{i}\mid \text{time}\right)\hfill \\ & & =\prod _{i=1}^{N}\frac{{\left(\mathrm{\tau }{f}_{i}\left(\text{time}\right)\right)}^{{n}_{i}}}{{n}_{1}!}\text{exp}\left(-\mathrm{\tau }{f}_{i}\left(\text{time}\right)\right)\hfill \end{array}$
where τ is the time window of sampling spike trains (20 ms, moving with 1-ms step, was used in this study), fi(time) is the PETH of i-th unit, ni is the number of spikes of i-th unit in the time window, and N is the total number of units. Combining these equations, the posterior probability of elapsed time was computed as
$P\left(\text{time}\mid \text{spikes}\right)=C·P\left(\text{time}\right)\left(\prod _{i=1}^{N}\frac{{\left(\mathrm{\tau }{f}_{i}\left(\text{time}\right)\right)}^{{n}_{i}}}{{n}_{i}!}\right)\text{exp}\left(-\mathrm{\tau }\sum _{i=1}^{N}{f}_{i}\left(\text{time}\right)\right)$
where C is a normalization factor that depends on τ and the numbers of spikes of each neuron (27).

Then, we assessed the relationships between decoded time information and theta phases (i.e., theta sequence). Decoded probabilities computed from Eq. 1 over trials in each block were aligned to the theta phases calculated from concurrently recorded LFPs (x axis) and the current time (y axis) and averaged to construct $P\left(\stackrel{^}{\text{time}}\mid \text{phase}\right)$, which is the probability density of relative time information decoded from spiking activity at a given theta phase; $\stackrel{^}{\text{time}}$ is the relative elapsed time and phase is theta phase (Fig. 4B) (25, 28). In this analysis, we analyzed theta sequence events in which more than three neurons were activated in single theta cycles.

In Fig. 4C, we assessed the compression rate of theta sequences between blocks 1 and 2 from each session. We designed the model of the theta sequence of block 2 by scaling the time domain of the theta sequence from block 1 using the following equation

$f\left(\stackrel{^}{\text{time}}\mid \text{phase}\right)={P}_{\text{block}1}\left(a·\stackrel{^}{\text{time}}\mid \text{phase}\right)$
where a is a parameter of scale factor and ${P}_{\text{block}1}\left(\stackrel{^}{\text{time}}\mid \text{phase}\right)$ is the theta sequence from block 1. We estimated the parameter (a) to minimize the sum of squared residuals between the real data and the model using the following equation
$\text{SSR}=\sum _{\text{time}}\sum _{\text{phase}}{\left({P}_{\text{block}2}\left(\stackrel{^}{\text{time}}\mid \text{phase}\right)-f\left(\stackrel{^}{\text{time}}\mid \text{phase}\right)\right)}^{2}$
where ${P}_{\text{block}2}\left(\stackrel{^}{\text{time}}\mid \text{phase}\right)$ is the theta sequence from block 2. We also estimated the compression rate of theta sequences between blocks 1 and 3 in the same way and compared the scale factors across the blocks. For these analyses, we selected the sessions in which more than six simultaneously recorded pyramidal neurons had phase precession (n = 26 sessions).

### Reconstruction of time information from neuronal activity in the test trials (task 4; Fig. 5)

In task 4, we performed the temporal bisection task with the insertion of test trials. Along with the normal trials of short- and long-interval times, we randomly inserted test trials. The geometric means of the short and long intervals were used as the time intervals during the test trials, and the rats were also required to select the left or right arm (associated with long or short intervals in the normal trials, respectively) of the Y-maze after forced running; however, a reward was omitted in both arms. We assessed differences in neuronal activity between the test trials in which the rats selected the arm that was correct in the long-interval trails (selected long) and the arm that was correct in the short-interval trials (selected short).

To this end, we decoded the elapsed-time information during the interval period from the spike trains separately in the test trials using the Bayesian decoding algorithm described above. We computed the posterior probability P (time|spikes) (Eq. 1) separately in the selected-long and selected-short test trials. Here, we used fi(time) as the PETH of i-th unit estimated in the normal trials and ni as the number of spikes of i-th unit in the test trials during the last 1 s before the end of the running periods (6 to 7 s for block 1 and 12 to 14 s for block 2). Then, decoded probabilities over trials in each block were aligned to the theta phases calculated from concurrently recorded LFPs (x axis) and averaged to construct P (time|phase), which is the probability density of elapsed-time information decoded from spiking activity during the last 1 s before the end of the running periods at a given theta phase (Fig. 5, B and D). To evaluate differences in the theta sequences between selected-long and selected-short test trials, we first estimated the theta sequences of the selected-long and selected-short trials from the spiking activity during every 1-s duration in the block 1 interval period in each session. Then, we computed the peak times of the averaged theta sequences of the two types of trials in each 1-s duration (n = 13 sessions; Fig. 5C). For block 2, theta sequences were estimated from the spiking activity during every 2-s duration in the block 2 interval time in each session. Then, we computed the peak times of the averaged theta sequences of the two types of trials in each 2-s duration (n = 13 sessions; Fig. 5E). For statistical comparison of the data between selected-long and selected-short trials, we applied two-way ANOVA followed by Bonferroni posttest using R software (57).

## Acknowledgments

We thank the RIKEN Advanced Manufacturing Support Team for supporting the fabrication of the custom-made experimental apparatus. Funding: This study was supported by the JSPS KAKENHI grants (15H05876, 16H01519, 18H02711, and 18H05525 for S.F.; 20K16480 for A.S.), the Mitsubishi Foundation (S.F.), the RIKEN JRA fellowship (A.S.), and the Keio University Doctorate Student Grant-in-Aid Program (A.S.). Author contributions: A.S., E.-I.I., and S.F. designed the experiments. A.S. performed experiments. A.S. and S.F. performed data analysis. A.S., E.-I.I., and S.F. wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

## REFERENCES AND NOTES

E. C. Tolman, Cognitive maps in rats and men. Psychol. Rev. 55, 189208 (1948).

T. E. J. Behrens, T. H. Muller, J. C. R. Whittington, S. Mark, A. B. Baram, K. L. Stachenfeld, Z. Kurth-Nelson, What is a cognitive map? Organizing knowledge for flexible behavior. Neuron 100, 490509 (2018).

J. O’Keefe, L. Nadel, The Hippocampus as a Cognitive Map (Oxford University Press, 1978).

J. O’Keefe, J. Dostrovsky, The hippocampus as a spatial map. Preliminary evidence from unit activity in the freely-moving rat. Brain Res. 34, 171175 (1971).

J. O’Keefe, M. L. Recce, Phase relationship between hippocampal place units and the EEG theta rhythm. Hippocampus 3, 317330 (1993).

J. Huxter, N. Burgess, J. O’Keefe, Independent rate and temporal coding in hippocampal pyramidal cells. Nature 425, 828832 (2003).

G. Buzsaki, Theta oscillations in the hippocampus. Neuron 33, 325340 (2002).

W. E. Skaggs, B. L. McNaughton, M. A. Wilson, C. A. Barnes, Theta phase precession in hippocampal neuronal populations and the compression of temporal sequences. Hippocampus 6, 149172 (1996).

G. Dragoi, G. Buzsaki, Temporal encoding of place sequences by hippocampal cell assemblies. Neuron 50, 145157 (2006).

10

D. J. Foster, M. A. Wilson, Hippocampal theta sequences. Hippocampus 17, 10931099 (2007).

11

J. O’Keefe, N. Burgess, Geometric determinants of the place fields of hippocampal neurons. Nature 381, 425428 (1996).

12

K. M. Gothard, W. E. Skaggs, B. L. McNaughton, Dynamics of mismatch correction in the hippocampal ensemble code for space: Interaction between path integration and environmental cues. J. Neurosci. 16, 80278040 (1996).

13

K. Diba, G. Buzsaki, Hippocampal network dynamics constrain the time lag between pyramidal cells across modified environments. J. Neurosci. 28, 1344813456 (2008).

14

H. Eichenbaum, Hippocampus: Cognitive processes and neural representations that underlie declarative memory. Neuron 44, 109120 (2004).

15

M. W. Howard, C. J. MacDonald, Z. Tiganj, K. H. Shankar, Q. Du, M. E. Hasselmo, H. Eichenbaum, A unified mathematical framework for coding time, space, and sequences in the hippocampal region. J. Neurosci. 34, 46924707 (2014).

16

G. Buzsaki, D. Tingley, Space and time: The hippocampus as a sequence generator. Trends Cogn. Sci. 22, 853869 (2018).

17

E. Pastalkova, V. Itskov, A. Amarasingham, G. Buzsaki, Internally generated cell assembly sequences in the rat hippocampus. Science 321, 13221327 (2008).

18

C. J. MacDonald, K. Q. Lepage, U. T. Eden, H. Eichenbaum, Hippocampal “time cells” bridge the gap in memory for discontiguous events. Neuron 71, 737749 (2011).

19

B. J. Kraus, R. J. Robinson II, J. A. White, H. Eichenbaum, M. E. Hasselmo, Hippocampal “time cells”: Time versus path integration. Neuron 78, 10901101 (2013).

20

Y. Wang, S. Romani, B. Lustig, A. Leonardo, E. Pastalkova, Theta sequences are essential for internally generated hippocampal firing fields. Nat. Neurosci. 18, 282288 (2015).

21

J. Csicsvari, H. Hirase, A. Czurko, G. Buzsaki, Reliability and state dependence of pyramidal cell-interneuron synapses in the hippocampus: An ensemble approach in the behaving rat. Neuron 21, 179189 (1998).

22

S. Terada, Y. Sakurai, H. Nakahara, S. Fujisawa, Temporal and rate coding for discrete event sequences in the hippocampus. Neuron 94, 12481262.e44 (2017).

23

W. E. Skaggs, B. L. McNaughton, K. M. Gothard, E. J. Markus, An information-theoretic approach to deciphering the hippocampal code, in Advances in Neural Information Processing Systems Vol. 5, S. J. Hanson, J. D. Cowan, C. J. Giles, Eds. (Morgan Kaufmann, 1993), pp. 1030–1037.

24

R. Schmidt, K. Diba, C. Leibold, D. Schmitz, G. Buzsaki, R. Kempter, Single-trial phase precession in the hippocampus. J. Neurosci. 29, 1323213241 (2009).

25

T. Feng, D. Silva, D. J. Foster, Dissociation between the experience-dependent development of hippocampal theta sequences and single-trial phase precession. J. Neurosci. 35, 48904902 (2015).

26

S. J. Middleton, T. J. McHugh, Silencing CA3 disrupts temporal coding in the CA1 ensemble. Nat. Neurosci. 19, 945951 (2016).

27

K. Zhang, I. Ginzburg, B. L. McNaughton, T. J. Sejnowski, Interpreting neuronal population activity by reconstruction: Unified framework with application to hippocampal place cells. J. Neurophysiol. 79, 10171044 (1998).

28

T. J. Davidson, F. Kloosterman, M. A. Wilson, Hippocampal replay of extended experience. Neuron 63, 497507 (2009).

29

B. E. Pfeiffer, D. J. Foster, Hippocampal place-cell sequences depict future paths to remembered goals. Nature 497, 7479 (2013).

30

K. D. Harris, J. Csicsvari, H. Hirase, G. Dragoi, G. Buzsaki, Organization of cell assemblies in the hippocampus. Nature 424, 552556 (2003).

31

T. J. Wills, C. Lever, F. Cacucci, N. Burgess, J. O’Keefe, Attractor dynamics in the hippocampal representation of the local environment. Science 308, 873876 (2005).

32

J. K. Leutgeb, S. Leutgeb, A. Treves, R. Meyer, C. A. Barnes, B. L. McNaughton, M. B. Moser, E. I. Moser, Progressive transformation of hippocampal neuronal representations in “morphed” environments. Neuron 48, 345358 (2005).

33

L. L. Colgin, E. I. Moser, M.-B. Moser, Understanding memory through hippocampal remapping. Trends Neurosci. 31, 469477 (2008).

34

G. Buzsaki, E. I. Moser, Memory, navigation and theta rhythm in the hippocampal-entorhinal system. Nat. Neurosci. 16, 130138 (2013).

35

H. T. Blair, A. C. Welday, K. Zhang, Scale-invariant memory representations emerge from moire interference between grid fields that produce theta oscillations: A computational model. J. Neurosci. 27, 32113229 (2007).

36

N. Burgess, C. Barry, J. O’Keefe, An oscillatory interference model of grid cell firing. Hippocampus 17, 801812 (2007).

37

M. E. Hasselmo, Grid cell mechanisms and function: Contributions of entorhinal persistent spiking and phase resetting. Hippocampus 18, 12131229 (2008).

38

K. H. Shankar, M. W. Howard, A scale-invariant internal representation of time. Neural Comput. 24, 134193 (2012).

39

A. S. Gupta, M. A. van der Meer, D. S. Touretzky, A. D. Redish, Segmentation of spatial experience by hippocampal theta sequences. Nat. Neurosci. 15, 10321039 (2012).

40

J. R. Huxter, T. J. Senior, K. Allen, J. Csicsvari, Theta phase-specific codes for two-dimensional position, trajectory and heading in the hippocampus. Nat. Neurosci. 11, 587594 (2008).

41

S. McKenzie, G. Buzsáki, Hippocampal mechanisms for the segmentation of space by goals and boundaries, in Micro-, Meso-and Macro-Dynamics of the Brain (Springer, ed. 21, 2016), chap. 1.

42

C. V. Buhusi, W. H. Meck, What makes us tick? Functional and neural mechanisms of interval timing. Nat. Rev. Neurosci. 6, 755765 (2005).

43

J. G. Heys, D. A. Dombeck, Evidence for a subcircuit in medial entorhinal cortex representing elapsed time during immobility. Nat. Neurosci. 21, 15741582 (2018).

44

B. J. Kraus, M. P. Brandon, R. J. Robinson II, M. A. Connerney, M. E. Hasselmo, H. Eichenbaum, During running in place, grid cells integrate elapsed time and distance run. Neuron 88, 578589 (2015).

45

G. B. Mello, S. Soares, J. J. Paton, A scalable population code for time in the striatum. Curr. Biol. 25, 11131122 (2015).

46

S. Soares, B. V. Atallah, J. J. Paton, Midbrain dopamine neurons control judgment of time. Science 354, 12731277 (2016).

47

A. Tsao, J. Sugar, L. Lu, C. Wang, J. J. Knierim, M. B. Moser, E. I. Moser, Integrating time from experience in the lateral entorhinal cortex. Nature 561, 5762 (2018).

48

J. J. Paton, D. V. Buonomano, The neural basis of timing: Distributed mechanisms for diverse functions. Neuron 98, 687705 (2018).

49

W. H. Meck, R. M. Church, D. S. Olton, Hippocampus, time, and memory. Behav. Neurosci. 98, 322 (1984).

50

S. K. Tam, C. Bonardi, Dorsal hippocampal involvement in appetitive trace conditioning and interval timing. Behav. Neurosci. 126, 258269 (2012).

51

R. M. Church, M. Z. Deluty, Bisection of temporal intervals. J. Exp. Psychol. Anim. Behav. Process. 3, 216228 (1977).

52

A. Berenyi, Z. Somogyvari, A. J. Nagy, L. Roux, J. D. Long, S. Fujisawa, E. Stark, A. Leonardo, T. D. Harris, G. Buzsaki, Large-scale, high-density (up to 512 channels) recording of local circuits in behaving animals. J. Neurophysiol. 111, 11321149 (2014).

53

S. N. Kadir, D. F. Goodman, K. D. Harris, High-dimensional cluster analysis with the masked EM algorithm. Neural Comput. 26, 23792394 (2014).

54

K. Mizuseki, A. Sirota, E. Pastalkova, G. Buzsaki, Theta oscillations provide temporal windows for local circuit computation in the entorhinal-hippocampal loop. Neuron 64, 267280 (2009).

55

P. Berens, CircStat: A MATLAB toolbox for circular statistics. J. Stat. Softw. 31, 121 (2009).

56

T. Danjo, T. Toyoizumi, S. Fujisawa, Spatial representations of self and other in the hippocampus. Science 359, 213218 (2018).

57

R. C. Team, R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2013).
https://creativecommons.org/licenses/by-nc/4.0/This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

Citing articles via