## Abstract

Single hippocampal cells encode the spatial position of an animal by increasing their firing rates within “place fields,” and by shifting the phase of their spikes to earlier phases of the ongoing theta oscillations (theta phase precession). Whether other forms of spatial phase changes exist in the hippocampus is unknown. Here, we used high-density electrophysiological recordings in mice of either sex running back and forth on a 150-cm linear track. We found that the instantaneous phase of spikes shifts to progressively later theta phases as the animal traverses the place field. We term this shift theta “phase rolling.” Phase rolling is opposite in direction to precession, faster than precession, and occurs between distinct theta cycles. Place fields that exhibit phase rolling are larger than nonrolling fields, and in-field spikes occur in distinct theta phases in rolling compared with nonrolling fields. As a phase change associated with position, theta phase rolling may be used to encode space.

**SIGNIFICANCE STATEMENT** Theta phase precession is a well-known coding scheme in which neurons represent the position of the animal by the timing of their spikes with respect to the phase of ongoing theta oscillations. Here, we show that hippocampal neurons also undergo “theta phase rolling,” a phase change faster and opposite in direction to precession. As the animal advances in space, spikes occur at progressively later phases of consecutive theta cycles. Future studies may reveal whether phase rolling constitutes a novel coding mechanism of space.

## Introduction

How neurons encode information is a major question in systems neuroscience (Bialek et al., 1991; deCharms and Zador, 2000; Stanley, 2013). Single neurons carry information by their firing rates (Adrian and Zotterman, 1926; Barlow, 1972; Azarfar et al., 2018) or by the precise timing of spikes in relation to a reference. The reference can be another spike of the same neuron (Reich et al., 2000; Imaizumi et al., 2010) or an external event (Optican and Richmond, 1987; Chase and Young, 2007). Neuronal assemblies also carry information by spatiotemporal patterns of activation (Gray et al., 1989; Graf et al., 2011; Kohn et al., 2016). Unfortunately, only a small fraction of all neurons in each network can be monitored at a given time. However, the activity of a population of neurons can be estimated by local field potentials (LFPs), a measure of the summed electric fields created by the synaptic activity of many neurons (Buzsáki et al., 2012; Herreras, 2016). LFP signals often consist of periodic oscillations at discrete frequency bands (Buzsáki and Draguhn, 2004; Singer, 2018) which can be used to define an internal reference to spikes. In such “phase coding,” information is carried by the timing of the spikes of a single neuron in relation to ongoing LFP oscillations (Montemurro et al., 2008; Turesson et al., 2012; Siadatnejad et al., 2013). Thus, a phase code is a type of a temporal assembly code, in which the activity of one neuron is referred to the synchronized activity of many hidden neurons, represented by the LFP oscillations.

Rate and phase coding by spiking neurons have been studied extensively in the context of spatial representation in the hippocampus (Hartley et al., 2014). Hippocampal neurons called “place cells” encode space by firing at elevated rates when the animal is in a specific region, the place field (O'Keefe and Dostrovsky, 1971; Moser et al., 2008). Many cells also encode information by the timing of spikes in relation to the instantaneous phase of hippocampal theta oscillations (Vanderwolf, 1969; Lever et al., 2014). As the animal traverses the place field, spikes occur at progressively earlier theta phases (O'Keefe and Recce, 1993; Skaggs et al., 1996). This phase code, termed “theta phase precession,” carries spatial information in addition to that carried by the rate code, improving the accuracy of the inferred position of the animal (Jensen and Lisman, 2000; Reifenstein et al., 2012). Originally described in CA1 pyramidal cells (PYR), precession has also been described in CA1 interneurons (INT; Maurer et al., 2006; Ego-Stengel and Wilson, 2007), and is found in other hippocampal regions (O'Keefe and Recce, 1993; Tsodyks et al., 1996) as well as in nonhippocampal regions (Van Der Meer and Redish, 2011; Tingley et al., 2018).

To date, theta phase precession is the only thoroughly studied spatial theta phase code. Yet a priori, multiple other phase codes are possible. Even precession itself is not a purely linear phase shift and is often asymmetric. In the second half of the place field, the slope of precession decreases (Mehta et al., 2002) and the variability of the phase increases (Cei et al., 2014; Souza and Tort, 2017). In some cases, the slope of precession reverses to generate “procession” in the second part of a place field (Yamaguchi et al., 2002; Wang et al., 2020). Interneurons occasionally exhibit a positive phase shift, in which theta phase advances forward (Ego-Stengel and Wilson, 2007). Positive phase shifts were also reported in small fractions of cells recorded from the lateral septum and the hippocampus (Monaco et al., 2019).

To determine whether there are additional phase-position patterns in CA1, we used high-density extracellular recordings in freely-moving mice. We found that CA1 neurons undergo a previously undescribed phase shift that we term “theta phase rolling.” Akin to precession, phase rolling occurs with respect to the ongoing theta oscillations but is faster and opposite in direction to precession.

## Materials and Methods

### Experimental design

#### Experimental animals

Three freely moving C57BL/NJ mice and two hybrid C57BL/NJ x FVB/NJ mice (mC41 and mA234) were used in this study (Table 1). Four of the mice were males and one was female, aged 10–30 weeks at the time of implantation. Animals were healthy, were not involved in previous procedures, and weighed 24.2–33.7 g at the time of implantation. Mice were single housed to prevent damage to the implanted apparatus. All animal handling procedures were in accordance with Directive 2010/63/EU of the European Parliament, complied with Israeli Animal Welfare Law (1994), and approved by Tel Aviv University Institutional Animal Care and Use Committee (IACUC; #01-16-051).

#### Probes and surgery

Every animal was implanted with a multi-shank silicon probe attached to a movable microdrive and coupled with optical fibers following previously described procedures (Stark et al., 2012; Noked et al., 2021). The probes used were Stark64 (Diagnostic Biochips), Buzaski32 (NeuroNexus), and Dual-sided64 (Diagnostic Biochips). The Stark64 probe consists of six shanks, spaced horizontally 200 µm apart, with each shank consisting of 10–11 recording sites, spaced vertically 15 µm apart. The Buzaski32 probe consists of four shanks, spaced horizontally 200 µm apart, with each shank consisting of eight recording sites, spaced vertically 20 µm apart. The Dual-sided64 probe consists of two dual-sided shanks, spaced horizontally 250 µm apart, with each shank consisting of 16 channels on each side (front and back), spaced vertically 20 µm apart.

Probes were implanted in the neocortex above the right hippocampus (PA/LM, 1.6/1.1 mm; 45° angle to the midline) under isoflurane (1%) anesthesia. After recovery from anesthesia, animals were placed on a water-restriction schedule that guaranteed at least 40 ml/kg of water (corresponding to 1 ml per 25-g mouse) on every recording day. Recordings were conducted 5 d/week, and animals received free water on the sixth day. After every one to five recording sessions, the probe was translated vertically downwards by up to 70 µm. Analyses included only recordings from the CA1 pyramidal cell layer, recognized by the appearance of multiple high-amplitude units and iso-potential spontaneous ripple events.

#### Behavior and recording sessions

Neuronal activity was recorded in 4.5 [1.9 10.1] h sessions (median [interquartile range, IQR]). At the beginning of every session, neural activity was recorded while the animal was in the home cage. The animal was then placed on a 150-cm linear track that extended between two 10 × 10-cm square platforms. Each platform included a water delivery port. Mice were under water restriction and were trained to repeatedly traverse the track for a water reward of 3–10 µl. Mice ran 167 [131.5 200.24] one-direction trials over about 1 h (Table 1). Trials in which the mean running speed was below 10 cm/s were excluded from analyses. Animals were equipped with a 3-axis accelerometer (ADXL-335, Analog Devices) for monitoring head movements. Head position and orientation were tracked in real time using two head-mounted LEDs, a machine vision camera (ace 1300-1200uc, Basler), and a dedicated system (“Spotter,” Gaspar et al., 2019).

#### Spike detection and sorting

Neural activity was filtered, amplified, multiplexed, and digitized on the headstage (0.1–7500 Hz, x192; 16 bits, 20 kHz; RHD2132 or RHD2164, Intan Technologies), and then recorded by an RHD2000 evaluation board (Intan Technologies). Offline, spikes were detected and sorted into single units automatically using KlustaKwik3 (Kadir et al., 2014; Rossant et al., 2016) for shanks with up to 11 sites/shank, or KiloSort2 (Pachitariu et al., 2016) for 16 channel shanks. Automatic spike sorting was followed by manual adjustment of the clusters. Only well-isolated units were used for further analyses [amplitude >40 µV; L-ratio <0.05 (Schmitzer-Torbert et al., 2005); inter-spike interval index <0.2 (Fee et al., 1996)]. Units were classified into putative PYR or PV-like INT using a Gaussian mixture model (Stark et al., 2013). A total of 5661 PYRs and 909 INTs were recorded from CA1 of the five freely moving mice.

#### Linear track firing rates

To create trial-averaged firing rate maps (Fig. 1*B*), we defined a “trial” as a single run across the linear track, either left-to-right or right-to-left. Time segments that occurred while the speed of the animal was below 10 cm/s were omitted. The firing rate in each trial was binned into 2.5-cm spatial bins, and count and temporal occupancy maps were calculated and smoothed with a Gaussian kernel that had an SD of 5 cm. A firing rate map was constructed for each trial by dividing the count map by the temporal occupancy map. Finally, a trial-averaged firing rate map was constructed for each run direction separately by calculating the mean firing rate in every spatial bin, weighted by the value in the corresponding bin of the occupancy map.

#### Activity and stability criteria

Many of the recorded units did not fire while the animal was on the linear track. Other units exhibited unstable firing during track sessions. All the analyses in the current study were performed only on active and stable units. Units that fired at least 5 spikes in at least one 2.5-cm spatial bin, pooled over all trials, were considered active. A total of 3722/5661 PYR (61%) and 753/909 (79%) INT were active. To quantify stability, for each unit we calculated the rank correlation coefficient (cc) between the firing rate maps of all same-direction trial pairs. For every trial pair, statistical significance was estimated using a permutation test. The overall p-value was estimated by the geometric mean over all pairs in both running directions. Units with consistent cc were considered stable. A total of 1928/3722 (51%) active PYRs and 541/753 (72%) active INTs were stable.

#### Identification of place fields

Typically, place fields are defined as continuous regions in which the firing rate of a PYR is above 10–20% of the peak on-track firing rate (Schmidt et al., 2009; Mizuseki et al., 2011; Diamantaki et al., 2018). However, such definitions assume discrete firing fields and do not attach a probability value to each field. Furthermore, INT change their firing rate according to position, exhibiting “on” and “off” regions on the track (Maurer et al., 2006; Ego-Stengel and Wilson, 2007; Wilent and Nitz, 2007; Hangya et al., 2010). Yet peak firing rate methods do not allow the identification of INT spatial firing regions. To detect place fields in a rigorous statistical manner in both PYR and INT, we developed a method based on Poisson distributions for field detection. In the Poisson method, place fields are defined as regions in which the firing rate deviates from chance at a given level of probability. For that, an estimate of the on-track spontaneous mean firing rate λ of each unit is required. However, the mean firing rate cannot be estimated based on off-track spiking, since many units rarely fire off-track. On the other hand, the mean on-track firing rate is strongly biased upwards by spatial spiking. To balance these two opposing biases, we used an iterative approach for estimating λ:

Trial-averaged firing maps were calculated separately for right-to-left and left-to-right trials.

Initially, λ was defined as the mean firing rate over all right-to-left and left-to-right bins, weighted by the total occupancy time in each bin.

The mean firing rate in each bin was compared to λ. Outlier bins were defined as bins in which the Poisson probability to obtain the observed firing rates or higher (if the firing rate was larger than λ) or lower (if the firing rate was smaller than λ), was higher than expected by chance. The alpha level was Bonferroni corrected for the number of effectively independent bins, determined by the full-width at half-height of the autocorrelation of the firing-rate map.

λ was recalculated after excluding the outlier bins.

Steps 3–4 were repeated for five iterations. Initial testing showed that five iterations are sufficient for convergence.

After establishing λ for a given unit, we determined whether the unit exhibited one or more place fields. For that, we again estimated the Poisson probability to obtain the observed (or higher) mean firing rates in each spatial bin of the rate map, given λ. Place fields were identified as regions of 15–100 cm with above-chance firing rates. Finally, place fields that contained a total of <30 spikes were omitted from further analyses, since preliminary testing showed that circular-linear analyses of phase changes (see below) are unreliable when the number of spikes is smaller. Overall, a total of 1564 PYRs and 261 INTs place fields were identified in 1928 PYR and 541 INT units (Table 1).

As a control for results obtained using the Poisson field detection method, we also detected “classical” fields using a procedure employed in multiple previous studies (Schmidt et al., 2009; Mizuseki et al., 2011; Diamantaki et al., 2018; Sharif et al., 2021). Classical fields were defined as (1) continuous regions spanning at least 15 cm; (2) in which the firing rate is above 10% of the peak rate in the field; and (3) exhibit a peak firing rate of at least two spikes per second. The classical method yielded 1636 PYR fields (*p* = 0.93, *G* test). Out of 1564 PYR fields identified by the Poisson method, 1009 (65%) overlapped fields identified by the classical method by at least 50% of the Poisson field. As expected, the two methods yielded considerably different INT fields: the classical method identified only 22 fields in 541 INTs, whereas 261 fields were identified by the Poisson method. To determine whether the detection of phase rolling (see below) depends on the field identification method, we repeated analyses for spikes within classical fields. 555/1636 (34%) of the classical fields exhibited phase rolling, compared with 540/1564 (34%) of the Poisson fields (*p* = 0.8, *G* test). Of the 540 Poisson fields that exhibited rolling, 432 (80%) also exhibited fields according to the classical method. Of these, 298/432 (69%) fields exhibited rolling also according to the classical method (chance level, 11.7%; *p* < 1.11 × 10^{−16}, exact binomial test). Furthermore, 258/1636 (16%) of the classical PYR fields exhibited between-cycle rolling, compared with 241/1564 (15%) of the Poisson fields (*p* = 0.81, *G* test). Of the 241 Poisson fields that exhibited between-cycle rolling, 184 (76%) also exhibited fields according to the classical method. Of these, 101/184 (55%) fields exhibited rolling by both methods (chance level, 2.4%; *p* < 1.11 × 10^{−16}, exact binomial test).

#### Instantaneous theta phase

To estimate the instantaneous phase of theta oscillations, one CA1 pyramidal cell layer channel was chosen for each session. The theta reference was the channel that exhibited the best signal-to-noise ratio (SNR), defined as the theta (5–11 Hz) to δ (2–4 Hz) power ratio (“SNR” method). The theta phase at each time point was calculated using the waveform method (Belluscio et al., 2012; Maurer et al., 2014). The LFP from the CA1 PYR channel was filtered with a zero-phase two-pole Butterworth bandpass (1–60 Hz) filter. The local maxima and minima of the filtered signal were determined, and extrema closer than 71 ms (14 Hz) were pruned. Each peak was assigned a phase of 0 rad, and each trough was assigned a phase of π rad. Phases at all other time points were interpolated linearly between peaks and troughs. To minimize discontinuities, the interpolated phases were smoothed using a Gaussian low-pass filter with a corner frequency of 13.6 [13 14.7] Hz (corresponding to half the mean session-specific theta cycle duration). Every spike was assigned the theta phase that occurred at the time of the spike, referred to as “spike phase.”

As a control for the theta channel selection method, we used a theta reference channel from the center of the CA1 pyramidal cell layer (“CA1pyr” method). The center of the layer was identified for each shank by the channel with iso-potential spontaneous ripple events (Ylinen et al., 1995). For each session, the shank with the maximal power between 130 and 190 Hz was selected. For on-track instantaneous phases, the mean ± SD difference between the theta reference based on the SNR and the CA1pyr methods was 0.047π ± 0.33π rad. Out of 1564 PYR fields, 540 (34%) exhibited phase rolling (see below) using the SNR method, compared with 526 (34%) using the CA1pyr method (*p* = 0.71, *G* test); 378/540 (70%) fields exhibited rolling by both methods (chance level, 11.6%; *p* < 1.11 × 10^{−16}, exact binomial test). Furthermore, 241 (15%) PYR fields exhibited between-cycle rolling using the SNR method, compared with 230 PYR fields using the CA1pyr method (15%; *p* = 0.9, *G* test); 109/241 (45%) fields exhibited between-cycle rolling by both methods (chance level, 2.3%; *p* < 1.11 × 10^{−16}, exact binomial test). There was no consistent difference in the fraction of PYR fields undergoing precession between the two methods (power: 799/1564, layer: 821/1564, *p* = 0.66, *G* test). In multiple previous reports (Hafting et al., 2008; Royer et al., 2012; Sharif et al., 2021), theta phase was estimated by applying a narrow bandpass filter (5–11 Hz) followed by Hilbert transform. To determine whether filtering specifics affect the detection of phase rolling, we repeated the rolling analysis following narrow-band filtering. Out of 1564 PYR fields, 540 (34%) fields exhibited rolling following 1- to 60-Hz filtering, compared with 595 (38%) following 5- to 11-Hz filtering (*p* = 0.16, *G* test). 395/540 (73%) fields exhibited rolling following filtering at both bands (chance level, 13.1%; *p* < 1.11 × 10^{−16}, exact binomial test); 241 (15%) PYR fields exhibited between-cycle rolling following 1- to 60-Hz filtering, compared with 229 (15%) following 5- to 11-Hz filtering (*p* = 0.93, *G* test). A total of 96/241 (40%) fields exhibited between-cycle rolling following filtering at both bands (chance level, 2.3%; *p* < 1.11 × 10^{−16}, exact binomial test). Thus, the precise theta reference channel and filtering range do not consistently affect rolling prevalence.

#### Theta phase lock, precession, and rolling

To quantify on-track spike theta phase locking within each place field, we calculated the mean phase of all the spikes that occurred in the field. Phase locking fit was quantified using the length of the circular resultant vector (

Theta phase precession and phase rolling were quantified for each place field by a circular-linear analysis (Schmidt et al., 2009; Kempter et al., 2012). To quantify theta phase changes, we fitted a line to the instantaneous theta phase of each spike, *a* was the slope that yielded the largest R, corresponding to the slope for which the variance of the residuals was smallest. The value of the maximal resultant length was used to indicate the fit of the circular-linear model. Statistical significance of the optimized model was estimated by randomly permuting the {

Preliminary observations suggested that the same unit may simultaneously exhibit theta phase locking, precession, and phase rolling. However, the circular-linear analysis yields a single optimum, potentially causing one phenomenon to shadow another. Circular-linear analysis is especially problematic in the presence of phase locking, as often observed for INT (Csicsvari et al., 1999; Frank et al., 2001; Klausberger et al., 2003; Somogyi and Klausberger, 2005). If the range of slopes used for the circular-linear analysis includes zero (a horizontal line), the zero-slope may yield the best fit for a phase-locked unit, even if that unit co-exhibits phase changes. To enable detecting multiple phase phenomena in a single field, we tested each place field for two circular-linear models. Negative slopes (corresponding to potential phase precession) were optimized in the range between tan(−0.1) and tan(−0.005) cycles/cm, and positive slopes (corresponding to potential phase rolling) were optimized between tan(0.04) and tan(0.25) cycles/cm. A field was determined as undergoing theta phase precession or rolling if a significant model was found within the negative or positive range of slopes, respectively (*p* < 0.05, permutation test). The number of phase precession/rolling cycles in every field was estimated by *a*Δ*x*, where *a* is the precession/rolling slope, and Δ*x* is the spatial distance between the first and last spike within the field.

A cycle randomization test was developed to assess the hypothesis of “within-cycle” phase rolling. For a tested place field, the phase and position of every spike within each theta cycle were randomized, yielding a set of surrogate phase-position pairs. For every spike, the surrogate pair was generated by randomly drawing a new time point within the cycle from a uniform distribution. The spike was assigned a new phase and position, corresponding to the phase and position values that actually occurred at that time point (Fig. 3*B*, gray lines). Thus, cycle randomization preserved the phase-position relations within every cycle. The randomization was conducted for every spike independently 1000 times. For each of the 1000 surrogate datasets, the fit of the phase-position pairs to the original circular-linear model was computed, generating a null distribution. The probability to obtain the original or higher fit under the null hypothesis was estimated empirically from the null distribution. If the probability was lower than the α level, the within-cycle (null) hypothesis was rejected, and phase rolling was classified as “between-cycle.” As an independent check, the intra-cycle spike phase randomization test was applied to estimate theta phase precession. Of all CA1 PYR place fields, 772/1564 (49%) exhibited between-cycle phase precession, a fraction not consistently different from all fields with precession, 799/1564 (51%, *p* = 0.58, *G* test).

### Statistical analyses

In all statistical tests used in this study, a significance threshold of α = 0.05 was used. All descriptive statistics (*n*, median, IQR) can be found in Results and figure legends. To estimate whether a given fraction was smaller or larger than expected by chance, an exact binomial test was used. Differences between the proportions of observations of two categorical variables were tested with a likelihood ratio (*G*) test of independence. Differences between two group medians were tested with Mann–Whitney *U* test. Differences between medians of three groups were tested with Kruskal–Wallis test. Uniformity of the distribution of circular parameters was tested with Rayleigh's likelihood-ratio test. Differences between the circular medians of two groups were tested with Wheeler's nonparametric two-sample test. Rank (Spearman's) and partial rank correlation coefficients were tested using permutation tests. In cases of multiple comparisons, Bonferroni's correction was employed.

## Results

### CA1 hippocampal cells undergo theta phase rolling

To quantify the relations between spike timing and the phase of the ongoing theta oscillations, we performed high-density extracellular recordings as five mice ran back and forth on a 150-cm-long linear track. Overall, the mice ran 167 [131.5 200.24] trials per session (median [IQR] of 97 sessions; Fig. 1*A*). We recorded a total of 6570 well-isolated units from hippocampal region CA1 of the animals. Of these, 5661 were putative pyramidal cells (PYR) and 909 were interneurons (INT; Table 1). Here, we focus on 1928/5661 (34%) PYRs and 541/909 (60%) INTs which were active on the linear track and exhibited stable spatial firing rate patterns across trials, henceforth referred to as “active units.” Most active units exhibited increased firing rates in a specific region of the track (Fig. 1*B*). The peak on-track firing rate was 8.3 [3.94 16.06] spikes/s for PYR and 32 [7.13 53.26] spikes/s for INT (*p* = 3.53 × 10^{−211}, Mann–Whitney *U* test). Previous work showed that INT change their firing rate according to position, exhibiting “on” and “off” regions on the track (Maurer et al., 2006; Ego-Stengel and Wilson, 2007; Wilent and Nitz, 2007; Hangya et al., 2010). We therefore identified “place fields” in both PYR and INT as regions spanning 15–100 cm in which firing rate increased compared with the on-track spontaneous firing rate (*p* < 0.05, Bonferroni-corrected Poisson test). Place fields were identified in 1019/1928 (53%) of the active PYR and in 207/541 of the active INT (38%; Fig. 1*C*), yielding a total of 1564 PYR and 261 INT fields. The fraction of units with one or more fields was larger in PYR than in INT (*p* = 0.00029, *G* test of independence). Considering all active units, PYR had 0.69 ± 0.017 place fields per unit, while INT had 0.48 ± 0.029 place fields per unit (mean ± SEM; *p* = 3.57 × 10^{−10}, *U* test; Fig. 1*Da*). Place fields spanned 42 [28 58] cm in PYR, and were smaller in INT, spanning 30 [20 40] cm (*p* = 1.2 × 10^{−19}, *U* test; Fig. 1*Db*). The in-field firing rate gain was higher among PYR (5.7 [1.8 21]) than among INT (1.3 [1.1 1.6]; *p* = 4.5 × 10^{−56}, *U* test). Within place fields, PYR tended to fire slightly after the theta trough (mean ± SD: 1.26π ± 1.93π rad; *p* = 1.9 × 10^{−17}, Rayleigh test), while INT tended to fire earlier (*p* = 2.22 × 10^{–16}, Wheeler's test), slightly before the theta trough (0.85π ± 1.36π rad; *p* = 6.1 × 10^{−19}, Rayleigh test; Fig. 1*E*). To conclude, in mouse CA1, spatial regions with increased firing rates (“place fields”) were observed in both PYR and INT. However, PYR place fields were larger and exhibited higher in-field firing rate gain.

Many place fields exhibit theta phase precession (O'Keefe and Recce, 1993; Skaggs et al., 1996; Harvey et al., 2009), in which spikes occur earlier in consecutive theta cycles as the animal advances within the place field (Fig. 2*A*). To quantify precession, we fitted a circular-linear model to the instantaneous theta phase of each spike and the spatial position of the animal at spike time (Schmidt et al., 2009; Kempter et al., 2012). For every place field (1564 PYR and 261 INT fields), the slope *a* of a circular-linear model was varied systematically and the resultant length (R) of the residuals was calculated (Fig. 2*D*). The slope of the model was determined by optimizing the fit to the data. Of the PYR with place fields, 799/1564 (51%) fields exhibited theta phase precession (*p* < 0.05, permutation test), a fraction higher than expected by chance (*p* < 1.11 × 10^{−16}, exact binomial test). In contrast, only 16/261 (6%) of the INT fields exhibited precession, a fraction not higher than expected by chance (*p* = 0.24, binomial test) and smaller than among PYR (*p* < 2.93 × 10^{−10}, *G* test; Fig. 2*E*). Thus, theta phase precession was observed in about half of the PYR place fields but not in INT. These observations expand previous work (Maurer et al., 2006; Ego-Stengel and Wilson, 2007), showing that INT undergo precession specifically in regions corresponding to place fields of upstream PYR.

In addition to phase precession, we found that many fields exhibited a positive slope of phase changes: spikes occurred at a later theta phase on one cycle, compared with the previous cycle (Fig. 2*B*,*C*). Thus, spike phase shifted forward as the animal traversed the place field. We term this phenomenon, which is opposite in direction to theta phase precession, theta “phase rolling.” To quantify theta phase rolling, we employed the same circular-linear model fitting procedure used for quantifying phase precession, extending the range of allowed slopes to positive values (Fig. 2*D*). Studying these plots, three categories of slopes could be observed: (1) negative slopes, corresponding to phase precession; (2) horizontal (zero) slopes, corresponding to locking of spikes to a specific theta phase; and (3) positive slopes, corresponding to phase rolling. Some fields exhibited several fit peaks, implying the co-occurrence of spike phase lock, precession and rolling. To examine rolling even in fields with strong phase lock or precession, we performed the circular-linear model fitting procedure for positive slopes separately. Theta phase rolling occurred in 540/1564 (35%) of the PYR fields and in 80/261 (31%) of the INT fields (*p* < 0.05, permutation test; Fig. 2*E*). Both fractions were above expected by chance (*p* < 1.11 × 10^{−16}, binomial test). Over the five mice tested, phase rolling prevalence was limited to the 26–43% range in PYR fields and to the 25–53% range in INT fields (Table 2). Unlike precession, the prevalence of phase rolling did not differ consistently between PYR and INT fields (*p* = 0.38, *G* test). Some fields co-exhibited phase lock, precession and rolling (Table 3). The occurrence of phase rolling and precession was not consistently associated (*p* = 0.82, *G* test). However, the occurrence of rolling was associated with the occurrence of phase lock (*p* = 0.0097, *G* test). In sum, theta phase rolling is a newly described phenomenon observed in about a third of the place fields in CA1.

### Phase rolling occurs within and between theta cycles

The phenomenon of theta phase rolling may be of two types. First, in a cell with inter-spike intervals shorter than the theta cycle duration, consecutive spikes that occur during the same theta cycle will always occur at larger theta phases. This results in trivial theta phase rolling. Furthermore, if theta phase is locked to a specific position across multiple trials (Agarwal et al., 2014), phase rolling will appear whenever spikes are pooled over trials, even if every cycle includes only one spike. Together, we term such phase rolling as occurring “within-cycle.” Second, phase rolling can also occur between theta cycles, if individual spikes (or spike bursts) in consecutive theta cycles occur at consecutively larger theta phases: in the same trial, or in different trials. We term this “between-cycle” phase rolling. Any observed phase rolling may originate from within-cycle rolling only, between-cycle rolling only, or both. Since within-cycle phase rolling can arise trivially while between-cycle cannot, we aimed to focus on between-cycle phase rolling.

To determine whether place fields exhibit nontrivial, “between-cycle” rolling, we developed a cycle randomization test. In the test, the phase and position of every spike are randomly drawn from the phases and positions that actually occurred in the relevant theta cycle (Fig. 3*A*,*B*, gray lines), generating surrogate phase-position pairs. The fit of the surrogate pairs to the original circular-linear rolling model is then computed. Because of the randomization, within-cycle phase-position relations are preserved, whereas relations between spikes occurring at distinct theta cycles are disrupted. If the probability to obtain the fit of the original spikes is lower than expected under the “within-cycle” (null) hypothesis, phase rolling is classified as “between-cycle.” Of all CA1 place fields, 241/1564 (15%) of PYR fields and 19/261 (7%) of INT fields exhibited between-cycle phase rolling (Fig. 3*C*). In PYR, the fraction of between-cycle phase rolling was above expected by chance (*p* < 1.11 × 10^{−16}, exact binomial test), but smaller than the fraction of all phase rolling fields (*p* < 3.33 × 10^{−16}, *G* test). In INT, the fraction of between-cycle phase rolling was not above expected by chance (*p* = 0.067, binomial test). The effect size η of between-cycle rolling was estimated by the fit of each rolling field (R), divided by the median fit over all cycle randomizations. η of PYR fields undergoing between-cycle rolling (1.5 [1.2 2]) was larger than that of within-cycle rolling fields (1 [0.94 1.1]; *p* = 1.4 × 10^{−63}, Bonferroni-corrected Kruskal–Wallis test; Fig. 3*D*). To conclude, rolling in INT occurs exclusively within cycles. In contrast, in over a third of the phase rolling PYR fields, theta phase rolling occurs between cycles.

Two independent tests were used to determine whether fields identified by cycle randomization undergo “trivial” rolling. First, we estimated the lock of instantaneous phase to position across trials of each field. If rolling occurs mainly within-cycles, fields undergoing rolling would be expected to exhibit stronger lock to position. However, the mean fit of instantaneous lock to position did not differ consistently between rolling and nonrolling fields (nonrolling fields median [IQR]: 0.18 [0.15 0.21]; within-cycle rolling fields: 0.18 [0.15 0.21]; between-cycle rolling fields: 0.18 [0.15 0.2]; *p* = 0.07, Kruskal–Wallis test; Fig. 3*E*). Second, to assess whether phase-position relations suffice to generate rolling, we constructed artificial spike trains that fire at fixed intervals of 60 Hz, interposed on phase and position data of every real rolling field (Fig. 3*F*). We then conducted the circular-linear analysis and cycle randomization test, exactly as for the real spike trains. Phase rolling was detected in 316/540 (59%) of the artificial trains, above expected by chance (*p* < 1.11 × 10^{−16}, binomial test). However, none of the artificial trains exhibited between-cycle rolling (Fig. 3*Fb*). Furthermore, effect sizes in fields undergoing between-cycle rolling (1.5 [1.3 2.2]) were larger than in fields with 60-Hz artificial spike trains (1 [1 1], *p* = 5 × 10^{−84}, Bonferroni-corrected Kruskal–Wallis test; Fig. 3*D*). Thus, phase lock to position does not predict between-cycle rolling, and high firing rate alone does not generate between-cycle rolling.

An obvious limitation of the cycle randomization test is disruption of the ISI structure of spikes within individual theta cycles. To evaluate the prevalence of between-cycle rolling in a manner that maximally conserves the ISI structure, we applied a pattern jitter test to all spikes in every place field (Harrison and Geman, 2009). Spikes with ISI ≤ 10 ms were grouped together, and the groups were jittered within a window of up to 126 ms, representing the mean duration of a single theta cycle (Fig. 4*A*). Spike groups were then interval-jittered 1000 times to generate a null distribution, and the probability to obtain the original (or higher) fit under the null hypothesis was estimated empirically from the null distribution. If the probability was lower than the α level, the within-cycle (null) hypothesis was rejected, and phase rolling was classified as “between-cycle.” Out of 1564 PYR fields, 233 (15%) exhibited between-cycle phase rolling according to the pattern jitter method, a fraction higher than expected by chance (*p* < 1.11 × 10^{−16}, exact binomial test). The fraction of fields exhibiting between-cycle rolling according to pattern jitter was not consistently different from that yielded by cycle randomization (241 fields, 15%; *p* = 0.73, *G* test; Fig. 4*B*). 191/233 (79%) fields exhibited between-cycle rolling according to pattern jitter and to cycle randomization (chance level, 2.3%; *p* < 1.11 × 10^{−16}, exact binomial test). To test the sensitivity of the pattern jitter results, we repeated the analysis with three other ISI durations (6, 8, 12 ms) and two other windows (half-cycle: 63 ms; double cycle: 252 ms; Fig. 4*C*). In all cases, the fraction of fields exhibiting between-cycle rolling was above chance (*p* < 1.11 × 10^{−16}, exact binomial test). Thus, an alternative randomization method which preserves the ISI structure yields similar results to the cycle randomization method.

### Theta phase rolling is a rapid phase change

To characterize between-cycle theta phase rolling, we compared several properties of between-cycle phase rolling and phase precession (see example units in Figs. 2*A*, 5*A*, recorded simultaneously on the same shank). The fit of spikes to the circular-linear model was 0.24 [0.18 0.31] in PYR fields undergoing precession, compared with 0.14 [0.11 0.18] in PYR fields undergoing phase rolling (*p* = 7.42 × 10^{−51}, *U* test; Fig. 5*B*). Compared with precession, between-cycle phase rolling exhibited a much faster slope with respect to spatial position (Fig. 5*C*). In PYR fields, phase rolling slope was 0.15 [0.068 0.19] cycles/cm, while precession slope size was 0.021 [0.015 0.03] cycles/cm (*p* = 2.16 × 10^{−112}, *U* test). The variance of slope size was higher for phase-rolling fields than for fields exhibiting phase precession (*p* = 0.002, resampling test). Corresponding to the difference in slope size, spiking underwent multiple phase rolling cycles in every PYR place field (6.4 [2.7 11] cycles/field; Fig. 5*D*). The occurrence of multiple rolling cycles is in sharp contrast to precession, in which spikes underwent less than a single cycle of precession within every field (0.75 [0.57 1] cycles/field; *p* = 3.86 × 10^{−102}, *U* test). Thus, compared with theta phase precession, between-cycle phase rolling is a faster phase change, spanning multiple cycles per field.

Rolling slope variance could originate from behavioral and anatomic differences within or between sessions. To examine whether precession and rolling slopes are stable over same-session trials, we divided the trials in every session into “early” and “late.” The precession slopes in the two halves of the session were correlated (747 fields, cc: 0.48, *p* = 0.002, permutation test), and did not differ consistently (*p* = 0.79, Wilcoxon's paired signed-rank test; Fig. 5*E*). Similarly, rolling slopes were correlated between the two halves (241 fields, cc: 0.3, *p* = 0.002), and did not differ consistently (*p* = 0.22). To examine possible anatomic sources of variability, we accumulated all possible pairs of between-cycle rolling fields (28 920 pairs) over all sessions. For every pair, we calculated the pairwise slope difference. Field pairs were split into three mutually exclusive groups: (1) recorded during distinct sessions (28,313 pairs); (2) recorded during the same session but on distinct shanks (433 pairs); and (3) recorded during the same session and on the same shank (174 pairs). There was no disparity in the pairwise slope difference between the three groups (median [IQR] in different sessions: 0.07 [0.028 0.12]; different shanks: 0.062 [0.027 0.12]; same shank: 0.06 [0.025 0.11]; *p* = 0.49, Kruskal–Wallis test; Fig. 5*F*). To conclude, while rolling slope sizes are more variable than precession slope sizes, rolling slope appears to be consistent between trials. Furthermore, locally co-recorded units do not appear to show particularly similar slopes.

### Phase rolling occurs in larger place fields and is centered around a distinct firing phase

To understand what differentiates PYR fields that undergo between-cycle phase rolling from fields that do not, we first considered the possible contribution of microcircuits along the deep-superficial axis of CA1 (Mizuseki et al., 2011; Geiller et al., 2017; Sharif et al., 2021). Fields were grouped to bins of 5 µm along the depth of the pyramidal layer, and the fraction of between-rolling fields out of all fields was calculated within each bin. The fraction of fields exhibiting between-cycle rolling was not consistently correlated with depth (cc = −0.1, *p* = 0.43, permutation test). Furthermore, there was no consistent correlation between rolling slope size and depth (cc: 0.086, *p* = 0.3, permutation test). Thus, phase rolling prevalence and slope do not vary consistently along the deep-superficial axis of CA1.

Next, we assessed the correlation between a host of cellular-network features and the occurrence of phase rolling. A total of 21 features were considered: behavioral, in-field LFP, spiking, network, on-track spiking, in-field spiking, and in-field spike-LFP features (Table 4; Fig. 6*A*). A multiple regression analysis of these features explained part of the variance associated with the occurrence of between-cycle phase rolling (*R*^{2} = 0.14, *p* = 0.00067, permutation test). Next, for each feature, the rank cc and the partial cc (pcc) were computed and Bonferroni corrected for multiple comparisons. Five features were consistently correlated with the occurrence of between-cycle phase rolling. First, on-track firing rates were not consistently correlated with phase rolling occurrence (cc: 0.042 ± 0.0024, *p* = 0.096), but did exhibit consistent partial correlation with rolling occurrence (pcc: 0.095 ± 0.0024, *p* = 0.00067). Second, phase rolling occurrence was correlated with the fit of spikes to the rolling slope (pcc: 0.31 ± 0.002, *p* = 0.00067, permutation test). The fit of spikes to the rolling slope was higher in fields with between-cycle phase rolling, 0.14 [0.11 0.18] compared with 0.12 [0.089 0.16] in other fields (*p* = 8.71 × 10^{−11}, *U* test; Fig. 6*B*). This observation is consistent with a monotonous relation between phase rolling fit and the probability of rolling detection. Third, phase rolling occurrence was negatively correlated with spike phase lock fit (pcc: −0.1 ± 0.0026, *p* = 0.00067). The fit (R) of spikes to a fixed theta phase was 0.11 [0.073 0.15] in fields with phase rolling and 0.13 [0.083 0.2] in fields without phase rolling (*p* = 8.49 × 10^{−7}; Fig. 6*C*). Indeed, if spike theta phase shifts within the field, the fit to a horizontal line (phase lock) decreases. Fourth, phase rolling occurrence was associated with distinct preferred firing phases (sine of preferred phase: pcc: −0.087 ± 0.0028, *p* = 0.00067). Specifically, fields with between-cycle phase rolling preferred firing at the middle of the rising theta cycle (mean ± SD: 1.48π ± 1.83π rad, R = 0.19, *p* = 0.00019, Rayleigh test), while nonrolling fields exhibited locking to an earlier phase (1.22π ± 1.91π rad, R = 0.16, *p* = 1.5 × 10^{−15}), slightly after the theta trough (*p* = 0.01, Wheeler test; Fig. 6*D*). Finally, phase rolling occurrence was correlated with place field size (pcc: 0.099 ± 0.0024, *p* = 0.00067). Fields with phase rolling were 50 [34 62] cm, larger than fields without rolling (40 [28 58] cm; *p* = 3.34 × 10^{−5}, *U* test; Fig. 6*E*). This suggests a relation between phase rolling and the way cells encode spatial position.

## Discussion

We found that putative pyramidal cells (PYR) and interneurons (INT) in mouse hippocampal region CA1 exhibit a previously undescribed type of phase changes with respect to space, theta phase rolling. In phase rolling, spike phase shifts forward between theta cycles as a function of space: as the animal traverses the place field, spikes on distinct theta cycles occur at gradually later phases. Theta phase rolling is a fast phase change that tends to occur in larger fields.

### Within-cycle versus between-cycle phase rolling

Theta phase rolling may appear because of the advancement of spikes with relation to theta phase either within or between cycles (Fig. 3). Within-cycle theta phase rolling can arise trivially, directly from the definition of phase, either if multiple spikes occur within a given theta cycle, or if theta phase itself is similar at a given position across multiple trials (Agarwal et al., 2014). In addition to within-cycle rolling, phase rolling may occur between spikes in distinct theta cycles, i.e., at distinct spatial positions. To determine whether place fields undergo between-cycle rolling, a cycle randomization test was developed, which identified fields in which rolling could not be explained by the within-cycle hypothesis. The results of the cycle-randomization test were corroborated by a pattern-jitter test, and by analyzing artificial spike trains. Using both cycle randomization and pattern jitter, we found that a total of 15% of the PYR fields exhibited between-cycle rolling. However, this may be considered a lower bound on the prevalence of phase rolling, since other fields may have exhibited both types of rolling.

### Possible functions of between-cycle phase rolling

Our results uncover a newly described phase change in which spikes shift to later theta phases as the animal moves forward in space. Previous studies have shown that a positive phase shift can follow the negative slope of phase precession in the second part of the place field (Yamaguchi et al., 2002; Wang et al., 2020). A positive phase shift was also described in other reports (Ego-Stengel and Wilson, 2007; Monaco et al., 2019). However, in the previous studies, positive phase changes exhibited slope sizes similar to those reported for precession. In contrast, the between-cycle phase rolling reported here is prevalent among PYRs (Fig. 3*C*) and, compared with phase precession, exhibits higher fidelity with respect to spatial position (Fig. 5).

Akin to place cell firing rates and to precession, theta phase rolling may contain information about the position of the animal. The instantaneous firing rate of a place cell indicates whether the animal is within the place field, and to some extent, where the animal is within the field (O'Keefe and Dostrovsky, 1971; Wilson and McNaughton, 1993; Brown et al., 1998; Moser et al., 2008). However, for single-peaked place fields, the instantaneous firing rate can only indicate that the animal is in one of two ranges. In place cells undergoing precession, the theta phase of spikes also contains information about position (O'Keefe and Recce, 1993; Skaggs et al., 1996), such that spike theta phase can indicate that the animal is within a single spatial range. The combination of theta phase with firing rate increases information about position (Jensen and Lisman, 2000; Reifenstein et al., 2012). Analogous to precession, phase rolling is a position dependent signal. Thus, we hypothesize that rolling by itself may reduce uncertainly about animal position. Distinctly from precession, a given phase rolling phase corresponds to multiple possible regions in space. Restricting the potential position of the animal to several short ranges may aid in estimating the position of the animal.

Following the observation of between-cycle phase rolling, we asked whether and what is unique in fields undergoing phase rolling. Phase rolling occurrence was found to be correlated with field size (Table 4; Fig. 6). Previous work showed that place fields are smaller and precession is more informative in cue-poor environments (Sharif et al., 2021). The association between larger field size and phase rolling may suggest that precession and rolling are distinct phase changes, exhibited in different environments. In small and/or cue-rich environments, place fields are smaller, and spike theta phase accurately encodes position. However, in larger and/or cue-poor environments, which typically exhibit larger place fields, rolling might increase spatial information by effectively partitioning the field into smaller subregions.

While the current study showed that phase rolling is correlated with animal position, we did not demonstrate that phase rolling contains spatial information or plays a causal role in spatial navigation. In future work, a functional role for rolling may be supported by a decoding analysis of animal position based on the instantaneous phase of phase-rolling fields, as has been previously done for precession (Jensen and Lisman, 2000). The functional role of rolling may be directly tested by manipulating spontaneous phase changes, or by inducing artificial phase changes and quantifying the influence on spatial behavior.

## Footnotes

This work was supported by the Israel Science Foundation Grant #638/16, the European Research Council Grant #679253, and by the Canadian Institutes of Health Research (CIHR), the International Development Research Centre (IDRC), the Israel Science Foundation (ISF) and the Azrieli Foundation Grant #2558/18. We thank Ortal Amber-Vitos for establishing and managing the animal colony and Leore Heim for help collecting preliminary data. We also thank Kamran Diba, Shaked Palgi, and Meishar Shahoha for constructive comments.

The authors declare no competing financial interests.

- Correspondence should be addressed to Eran Stark at eranstark{at}tauex.tau.ac.il