Abstract
Running direction in the hippocampus is encoded by rate modulations of place field activity but also by spike timing correlations known as theta sequences. Whether directional rate codes and the directionality of place field correlations are related, however, has so far not been explored, and therefore the nature of how directional information is encoded in the cornu ammonis remains unresolved. Here, using a previously published dataset that contains the spike activity of rat hippocampal place cells in the CA1, CA2, and CA3 subregions during free foraging of male Long-Evans rats in a 2D environment, we found that rate and spike timing codes are related. Opposite to a preferred firing rate direction of a place field, spikes are more likely to undergo theta phase precession and, hence, more strongly affect paired correlations. Furthermore, we identified a subset of field pairs whose theta correlations are intrinsic in that they maintain the same firing order when the running direction is reversed. Both effects are associated with differences in theta phase distributions and are more prominent in CA3 than in CA1. We thus hypothesize that intrinsic spiking is most prominent when the directionally modulated sensory-motor drive of hippocampal firing rates is minimal, suggesting that extrinsic and intrinsic sequences contribute to phase precession as two distinct mechanisms.
SIGNIFICANCE STATEMENT Hippocampal theta sequences, on the one hand, are thought to reflect the running trajectory of an animal, connecting past and future locations. On the other hand, sequences have been proposed to reflect the rich, recursive hippocampal connectivity, related to memories of previous trajectories or even to experience-independent prestructure. Such intrinsic sequences are inherently one dimensional and cannot be easily reconciled with running trajectories in two dimensions as place fields can be approached on multiple one-dimensional paths. In this article, we dissect phase precession along different directions in all hippocampal subareas and find that CA3 in particular shows a high level of direction-independent correlations that are inconsistent with the notion of representing running trajectories. These intrinsic correlations are associated with later spike phases.
Introduction
Hippocampal place cells establish a neuronal representation of space by exhibiting elevated firing rates at only few locations in an environment called place fields (O'Keefe and Dostrovsky, 1971). Place field firing is thought to underlie the capacity of an animal to navigate in space and to form spatial memories (Morris et al., 1982; Moser et al., 1993; Nakazawa et al., 2002). Place field firing also includes a temporal code associated with the theta oscillation (4–12 Hz) of the local field potential (O'Keefe and Recce, 1993): As a rat passes through a place field, the spikes phase precess, that is, they occur at successively earlier theta phases thereby encoding the relative location of the animal within a place field (O'Keefe and Recce, 1993; Harris et al., 2002). Phase precession is generally thought to implement a compression of behavioral sequences to the theta time scale so that in a 100 ms time window, spikes of multiple place cells are elicited in the same order as the activation of the associated place fields along the trajectory of the animal over the time scale of seconds (Melamed et al., 2004; Dragoi and Buzsáki, 2006; Foster and Wilson, 2007; Jaramillo and Kempter, 2017).
In one-dimensional mazes, the trajectory of an animal can be uniquely mapped to a sequence of place fields; thus, spike sequences on the theta time scale (Foster and Wilson, 2007) can be easily interpreted as reflecting memories of previous locations or planning of future actions (Feng et al., 2015). In two-dimensional environments, place fields can be entered from multiple directions, and hence place cells generally take part in encoding multiple trajectories. Thus, place cells could either be linked to multiple sequences, or place field sequences could be directional. In the former case, sequential structure would be imposed by sensory-motor (extrinsic) inputs, whereas the latter case would render sequences of intrinsic origin supported by recurrent circuits or associative loops. Previous reports revealed that pair correlation lags of place fields in the CA1 subregion depend on running direction (Huxter et al., 2008) and thus support the extrinsic hypothesis, but similar analyses for the CA2 and CA3 subregions, which differ substantially in their cytoarchitecture, plasticity, and protein chemistry, are missing.
In addition to sequence order, directional information is also available to the entorhinal-hippocampal circuits from head direction cells of the postsubiculum (Taube et al., 1990a,b) and the medial entorhinal cortex (Sargolini et al., 2006; Giocomo et al., 2014) and, to a limited extent, from within the cornu ammonis itself (Leutgeb et al., 2000). These inputs might explain observed directionality of some place fields (Leutgeb et al., 2004; Acharya et al., 2016; Mankin et al., 2019); however, it is unclear to which extent this rate directionality is related to directionality of theta sequences.
Past studies have shown that compared with the CA1 region, place cells in CA3 demonstrate a more persistent and consistent activity pattern over an extended period of time (Mankin et al., 2012) and across cue-altered environments (Lee et al., 2004), as well as a higher stability of place field dynamics across multiple recording sessions in the same environment (Mizuseki et al., 2012). Consistently, CA3 place representations stabilized more slowly after a change of environment than in CA1 (Leutgeb et al., 2004). The stability, consistency, and slower stabilization of CA3 place fields, in addition to the classical anatomy indicating strong recurrent connectivity (Amaral and Witter, 1989; Ishizuka et al., 1990), lead to a belief that CA3 activity patterns are more reliant on internal network dynamics and less influenced by external sensory inputs. Therefore, we hypothesized that the theta sequence activity in CA3 should be less dependent on the directionality in the behaviors of the animal than CA1.
Our results show that although the extrinsic contribution to theta scale firing is dominant in all subareas, this is indeed least visible in CA3. Moreover, we observe directionality in the phase precession properties so that CA3 displays later spike phases in the running direction opposite to the best firing rate direction.
Materials and Methods
Experimental design and statistical analysis
We reanalyzed a previously published dataset in Mankin et al. (2012, 2015). For a detailed description of the data collection, we refer to the original work. In brief, the dataset involves eight male Long-Evans rats that were trained to forage for randomly scattered cereal crumbs in either a square (80 cm × 80 cm) or a 16-sided polygon (50 cm radius, also referred to as circular) enclosure. The experiment began after animals were trained 9–20 d in the enclosure. Single units were recorded simultaneously from the CA1, CA2, and CA3 subregions for the course of the experiment, which lasted 2 d. On each day the rats completed two blocks of four 10 min sessions, with two sessions in the square enclosure and two sessions in the circular enclosure assigned in random order. The enclosures contain a 20 cm-wide white cue card on an inside wall, and the cue card maintained a constant angle with the cues outside the room.
Data analysis and statistical tests were performed on Python using SciPy (Virtanen et al., 2020), NumPy (Harris et al., 2020) packages, and custom routines based on CircStat toolbox (Philipp, 2009). We used nonparametric Kruskal–Wallis tests with a post hoc Dunn's test for statistical comparisons. Normally distributed data were tested using Student's t test. The Watson–Williams test was used for circular data. For categorical data, we used χ2 and Fisher's exact tests. The p-values were adjusted for multiple comparisons using the Benjamini–Hochberg procedure (Benjamini and Hochberg, 1995). We used two-tailed tests throughout except for binomial tests, and p = 0.05 is chosen as the significance level.
Place field detection and delimination
For each place cell, we computed a spatial map of firing rates from each 10 min session by dividing the spike counts by the occupancy time in each space bin (1 cm × 1 cm). Spikes and occupancy were smoothed by a Gaussian filter with an SD of 3 cm. Place fields were segmented along the closed contour line located at 20% of the maximum rate and accepted if the peak rate exceeds 1 Hz, the field area is larger than 25 cm2, and the average firing rate within the field is larger than outside the field. In the analysis, we separated place fields into a border and a nonborder group. The groups are distinguished by the 20% line touching the boundary of the enclosure. Each place cell could have multiple place fields that were analyzed separately.
Two place fields are said to be a pair if both field areas intersect and contain at least 16 spikes of one field that occur next to a spike of the other field within a time window of ±0.08 s (see Fig. 4B). We defined pairs as border pairs if at least one of the place fields touches the border. In nonborder pairs, both place fields do not touch the border.
Directionality of place fields
To obtain the directional tuning curve of place fields, we used a maximum likelihood maximization (MLM) model (Cacucci et al., 2004). In brief, the model assumes an independent relation between positional and directional firing probability distributions, whose product is the firing probability. The product assumption furthermore ensures that directionality tuning is restricted to the place field and that the tuning depends on a weighted sum of position and direction. The solution of the directional term, which is also the estimated directional tuning curve, can be fit by iteratively maximizing firing likelihood to the observation of spikes. This MLM model has an advantage of reducing sampling bias, which usually arises at the enclosure borders where certain heading directions are severely undersampled.
Significance of directionality is determined by comparing the mean resultant vector length (R) of the directional tuning curve to a shuffling distribution, which is obtained by randomly shifting the spike times in a cyclic fashion along a trajectory concatenated from all path segments traversing the place field. The random time shift was repeated 200 times for each place field. A place field is classified as significantly directional if the R value exceeds the 95th percentile of its shuffling distribution. For field pairs, the same method is used to determine the significance of directional selectivity, except that the R is calculated from spike pairs.
Significance of the preferred precession direction is also determined by comparison of R to a distribution obtained from shuffling spike times over all passes through the field. In this case the R value for precession directionality is computed from the directions of field traversals that exhibit phase precession. A pass is labeled as precessing if the phase-position relation of the spikes has a negative slope derived from linear-circular regression (Kempter et al., 2012).
Correlation lags
We produced cross-correlograms for the spike trains of every field pair, with the resolution of 5 ms and a time window between −150 ms and 150 ms. The resultant cross-correlograms were then bandpass filtered (5–12 Hz) to derive the correlation lag as the phase at 0 time lag from the Hilbert transform of the theta filtered correlograms.
Classification of extrinsic and intrinsic pairs
To quantify the dependence of correlation lag on directionality, we devised the measures of “extrinsicity” and “intrinsicity” for each field pair. The extrinsicity is computed as the Pearson's correlation coefficient between the cross-correlogram for runs from one field to another and the cross-correlogram for runs in the opposite direction but with the sign of time-axis flipped. The Pearson's correlation coefficient (r) is then linearly transformed (
Inclusion criteria for analysis
The animal trajectory is split to passes entering and exiting the place fields or pairs. Intervals in which the speed of the animal is below 5 cm/s were excluded. The passes are then chunked to smaller segments in which the speed was always above 5 cm/s.
The pass segments and their spikes are only included if the pass duration is longer than 0.4 s and satisfies a straightness threshold as in the following:
Pair-crossing passes are classified as A → B if they start from an area within field A but not field B, and end in an area within field B but not field A. The opposite criteria apply for B → A passes. Passes that do not satisfy the above criteria, or cross either one of the field boundaries more than once, are not assigned to any of the directional groups (A → B or B → A). These unassigned passes were also included in computing the firing rate directionality but excluded for cross-correlation analyses. As a result, field pairs that have no spike pairs along the A → B or B → A passes were further excluded in the cross-correlation analyses.
Model simulation
Romani and Tsodyks (2015) proposed a recurrent network model of the hippocampus, which we adapted and simulated using our trajectory data. Here, we briefly summarize the key equations of the model. The dynamics of firing rate mi(t) of place cell i at time t is given by the following:
The parameter 0 < U ≤ 1 denotes the release probability. The introduction of xi(t) penalizes the recurrent input from the highly activated place cells with a delay, and therefore it produces asymmetrical weight couplings that are stronger in the forward direction as the animal moves.
We adopted the model parameters from the toroidal environment described in Romani and Tsodyks (2015). The periodicity of the environment was removed by clipping the cosine function cos(x) at the value of −1 for |x| > π. Our simulation thus has 32 × 32 neurons equally spaced across a 2π × 2π unit squared environment. We randomly chose one recording session in the square arena from the experimental data as the trajectory, and rescaled the trajectory into the range of 0 and 2π and the average speed to be the same as in the original study (2π/5 unit per second). Simulation was implemented with 1 ms temporal resolution using the Euler method. Spikes from each place cell were then subsampled by a fraction so that the average spike count of all simulated place fields is the same as the experimental data.
Results
Directional tuning of place cells
Directionality of hippocampal place field activity has been reported in a number of previous studies (Leutgeb et al., 2004; Acharya et al., 2016; Mankin et al., 2019), but quantitative comparisons between those studies were hampered because of the use of different methods and behavioral paradigms. Here, we analyze the firing properties of simultaneously recorded CA1, CA2, and CA3 neural networks under identical experimental conditions. We thus first applied one established rate-based directionality analysis on the datasets used in this article (Mankin et al., 2012, 2015) for further comparison. Directional tuning for each place field was quantified using the mean resultant vector length (R) obtained from directionality fields derived by the MLM proposed in Cacucci et al. (2004); Fig. 1A, single examples, 1B populations. Because mean resultant vector lengths are strongly biased by the number of observations (Fig. 1B), we decided to analyze directional tuning as a function of the spike count threshold criterion for including place fields (Fig. 1C) and include only the fields with spike counts higher than 40 in our statistical analysis of place field directionality. We found that all CA regions contain a fraction of directional fields that is above chance level (Binomial test; CA1, 146/800 = 0.1825, p = 1.8e – 41; CA2, 111/521 = 0.2131, p = 2.4e – 38; CA3, 45/396 = 0.1136, p = 3.5e – 07). The amount of directionality in all regions not only depends on the overall spike count threshold of the place field (showing an initial increase that is expected from the gain in statistical power) but also on whether the place field is located at the boundary of the arena (Fig. 1C).
Directionality of place field firing rates. A, Examples of place fields in square and circular enclosures of a free-foraging experiment overlaid with the trajectory of the animal in one recording session and spike events color coded by heading direction (color wheel at right). Directional tuning curve is shown to the right of each spike position plot. Mean direction (black bar), mean resultant vector length R (top), peak rate in Hz (bottom). B, R and total within-field spike counts of all border (B, teal) and nonborder (N-B, gold) place fields, as well as marginal distribution of R (right). R values are strongly biased by sampling. Therefore, we excluded the fields with spike counts below 40 in our directionality analysis, indicated by the shaded region. C, Median R (top row) and fraction of significantly directionally tuned place fields (bottom) by different spike number thresholds for all (solid line), border (dotted line), and nonborder (dashed line) fields in each brain region as indicated. CA1 and CA2 directionality is strongly border driven. D, Fraction of all place fields (top) and border fields (bottom) by spike count thresholds.
Comparing the place field directionalities between CA regions, CA1 and CA3 have a higher median R than CA2 [Kruskal–Wallis test; CA1 (n = 800) vs CA2 (n = 521) vs CA3 (n = 396), H(2) = 22.18, p = 1.5e – 05; post hoc Dunn's test with Benjamini–Hochberg correction; CA1 vs CA2, p = 0.0003, CA2 vs CA3, p = 3.8e – 05, CA1 vs CA3, p = 0.1920]. The fraction of significantly directional fields in CA3 is lower than in CA1 and CA2 (Fisher's exact test for independence of significant fractions with Benjamini–Hochberg correction; CA1 vs CA2, p = 0.2656; CA1 vs CA3, p = 2.2e – 05; CA2 vs CA3, p = 0.0074). However, as we increase the spike count threshold to admit only fields that are highly sampled, the fraction of significantly directionally tuned fields becomes similar in CA1 and CA3, as far as our data allow such a comparison because of the only very few CA1 and CA3 fields with high spike numbers (Fig. 1D).
To accurately interpret the above results, we looked into possible confounds. A major influence on directionality could arise from the presence of arena boundaries, both because they act as salient sensory landmarks and they introduce a behavioral bias. We thus further separated place fields into border and nonborder fields. Including all the place fields regardless of spike counts, we found that border fields in CA1 generally have higher directional selectivity than in the nonborder case [Kruskal–Wallis test; border (n = 537) vs nonborder (n = 263), H(1) = 62.33, p = 2.9e – 15; Fisher's exact test; p = 5.7e – 05], whereas CA3 exhibits significant border difference in median R [border (n = 257) vs nonborder (n = 139), H(1) = 15.31, p = 9.1e = 05] but not in significant fraction (Fisher's exact test; p = 0.0885). Similar to CA1, directionality measures in CA2 also exhibit a significant border effect [Kruskal–Wallis test; border (n = 359) vs nonborder (n = 162), H(1) = 31.95, p = 1.6e – 08; Fisher's exact test; p = 3.6e – 07].
We thus conclude that place field rates in all CA areas encode running direction, and directionality in CA1 and CA2 is more strongly induced by borders, whereas this is not the case for CA3 in which the directionality is more similar between border and nonborder fields. Assuming boundaries to induce a strong sensory-motor constraint, this is already a first hint that CA1 activity is more strongly influenced by extrinsic factors than CA3.
Preferred direction for phase precession in place fields
In addition to the firing rate code, because place field activity is also temporally organized on the theta scale by phase precession (Fig. 2A), we also asked to what extent directionality is also reflected in this temporal code. For each place field, we therefore identified a direction in which phase precession is more likely to occur using single pass phase precession analysis (Schmidt et al., 2009; Kempter et al., 2012; see above, Materials and Methods for inclusion criteria). In brief, we fit a linear-circular regression line for the phase-position relation in every single pass. Passes with negative regression slope between –2π and 0 (per pass length) are classified as phase precessing (Fig. 2B).
Phase precession per spike is most prevalent opposite to the direction of highest firing rate. A, Top row, Three examples of phase precession pooled over all passes in one recording session. Position is normalized to be between 0 and 1 with respect to the moments of entering and exiting the place field. Linear-circular regression line is indicated in black, characterizing phase precession by its slope and onset phase. Bottom rows, Phase precession in example passes. Passes with slope between –2π and 0 are classified as incidents of phase precession. B, Fraction of fitted negative slopes among all passes as a function of pass direction θpass for CA1, CA2, and CA3 (color as indicated) indicates a similar amount of precessing passes in all subregions. C, Distribution of phase precessing passes in the whole dataset (left) mirrors the elevated spike count along the preferred rate direction θrate of the field (middle) pooled over all precessing passes as a function of absolute angular deviation |d(θpass, θrate)| from pass direction θpass. We therefore normalized the spike distribution by precession occurrences (ratio of left and middle graph) to obtain a distribution of phase precession per spike (right), which increases with angular distance |d(θpass, θrate)| (p-values from Spearman's correlations; CA1, rs = 0.83; CA2, rs = 0.39; CA3, rs = 0.62) indicating an excess of precession at |d(θpass, θrate)| = π that cannot just be explained by increased firing. D, The same analysis on a fieldwise level shows preferred precession directions (per spike) θprecess of single fields and best firing rate direction θrate to be offset by about π (marked by black line). E, Top, Same as D shown as normalized polar histograms of θprecess relative to θrate. Bottom, Only low-spike passes (25th percentile) are admitted (p-values are derived from V-test vs the null hypothesis of a circular mean at π) to control for high rate bias. Arrow marks direction of the mean resultant vector of the distribution, with the best rate direction pointing to the right. F, Cumulative distribution of R of all place fields in CA1, CA2, and CA3. Kruskal–Wallis test indicates strongest directionality in CA3.
First, we computed the density of phase precession occurrences from all passes in all fields as a function of pass direction relative to the preferred firing rate direction of the respective field (Fig. 2C). We find most precessing passes along the best rate direction, reflecting the fact that more spikes should give rise to more detectable phase precession. However, not all spikes may contribute to phase precession to the same degree, either because of different levels of phase noise or because they occur outside theta sequences. If phase precession is directional beyond a simple spike count effect, it needs to show in an analysis per spike. Thus, computing the density of precessing passes per spike (Fig. 2C, right), we found that phase precession is more likely to occur the more the pass direction differs from the preferred rate direction of the field (Spearman's correlation; CA1, rs(6719) = 0.83, p = 2.0e – 12; CA2, rs(4837) = 0.39, p = 0.0085; CA3, rs(2625) = 0.62, p = 7.5e – 06), indicating that spike rate and phase precession differentially contribute to rate directionality.
In addition to this population-wide analysis, we also identified the direction in which precession is most probable per spike for each field separately and call it the best precession angle θprecess. Histograms of preferred precession angles from all place fields separately (Fig. 2E, top row) also demonstrate a significant π shift from their preferred firing direction (V-test vs π direction; CA1, V(731) = 155.54, p = 2.2e – 16; CA2, V(466) = 51.75, p = 0.0003; CA3, V(341) = 47.12, p = 0.0002). However, by comparing the R values of precessing passes to a shuffling distribution (see above, Materials and Methods, shuffling procedures), we found that only 24/829 (2.9%), 11/560 (2.0%), and 20/441 (4.5%) of place fields in CA1, CA2, and CA3, respectively, exhibit significant preferred precession direction, which is not significant under binomial tests (CA1, p = 0.9990; CA2, p = 0.9999; CA3, p = 0.7034). We thus conclude that although on the level of the single field the antiphase relation between spike count and phase precession is weak and does not reach significance (and therefore has likely not been identified previously), there is a strong indication of such a relation on the population level.
To further rule out that the π shift between best rate and best precession direction might arise as an epiphenomenon of different spike counts, with the opposite of the preferred rate direction being overrepresented by the normalization process, we recomputed the histograms of preferred precession angles by limiting the passes to only those with low spike counts (<25% quantile of all precessing passes in each CA region) so that there is no firing rate directionality left in the used data. Our results show that the place fields in all CA regions still demonstrate a significant π shift from the preferred rate direction (Fig. 2E, bottom row; V-test vs π direction; CA1, V(294) = 46.58, p = 6.1e – 05; CA2, V(226) = 20.44, p = 0.0273; CA3, V(120) = 17.04, p = 0.0139). Thus, on the population level, the direction of best phase precession displays a consistent and significant bias toward the opposite of the direction of best firing rate, corroborating the hypothesis of distinct coding schemes and, hence, input streams, for spike timing and rate (Huxter et al., 2003).
Finally, we also compared distributions of Rayleigh vector lengths R for phase precession tuning and found that although CA1 and CA2 seem to have similar directionality, CA3 exhibits a significantly higher directional selectivity (Fig. 2F; Kruskal–Wallis test; CA1 (n = 753) vs CA2 (n = 485) vs CA3 (n = 363), H(2) = 9.37, p = 0.0092; post hoc Dunn's test with Benjamini–Hochberg correction; CA1 vs CA2, p = 0.3351; CA2 vs CA3, p = 0.0084; CA1 vs CA3, p = 0.0267), further suggesting that directional information in CA3 place field activity is distinct from CA1 and CA2.
Phase precession properties show dependence on pass direction
To further support the existence of directional effects on phase precession and to elucidate the underlying processes, we searched for directional modulations of phase precession properties by quantifying single pass onset phase and precession slope and plotted their occurrence density (Fig. 3A) for different relative pass directions |d| (defined as the absolute circular difference between pass direction and preferred rate direction of the place field). Although there was no significant correlation between precession slope and |d| in any CA region (CA1, r(5020) = –0.008, p = 0.5768; CA2, r(3514) = 0.031, p = 0.0564; CA3, r(2044) = 0.013, p = 0.5558), we found that phase onsets slightly but significantly increase as the heading deviates more from the preferred rate direction in CA3 (r(2625) = 0.078, p = 5.4e – 05), also, but barely significantly, in CA1 (r(6719) = 0.024, p = 0.0459), but not significantly in CA2 (r(4837) = –0.023, p = 0.0952). Thus, the opposite directions of best firing rate are signified by later phases, which could reflect more prospective parts of theta sequence activity (Foster and Wilson, 2007).
Directional dependence of phase precession properties. A, Marginal distributions of precession density as a function of onset phase and slope, color coded by the difference between pass angle θpass and best rate angle θrate (top, illustration of the color code). Black dot indicates the circular mean of marginal density. B, Average phase-position relations for cases where the animal is running along (blue) and against (green) θrate. Schematic illustration above. Against-θrate condition has higher onset than along-θrate condition; po denotes the p-value from Watson–Williams test for onset difference (along-θrate vs against-θrate, mean ± SEM in radians; CA1, 3.66 ± 0.04 vs 3.77 ± 0.05, F(1,2073) = 2.93, p = 0.0871; CA2, 3.78 ± 0.04 vs 3.78 ± 0.05, F(1,1301) = 0.00, p = 0.9978; CA3, 3.44 ± 0.06 vs 3.94 ± 0.07, F(1,848) = 25.34, p = 5.9e – 07). There is no difference in slopes between both cases; ps denotes the p-value from Kruskal–Wallis test for slope difference (mean ± SEM in radians per unit position; CA1, −4.62 ± 0.08 vs −4.54 ± 0.10, H(1) = 0.84, p = 0.3584; CA2, −4.48 ± 0.10 versus −4.57 ± 0.13, H(1) = 0.06, p = 0.8117; CA3, −4.44 ± 0.12 vs −4.20 ± 0.14, H(1) = 1.69, p = 0.1940). C, Histograms of spike phases from precession samples show higher spike phase for passes against θrate; p-values are derived from Watson–Williams test for the difference between the circular means (shown as vertical bars) between the two cases (along-θrate vs against-θrate, mean ± SEM in radians; CA1, 3.66 ± 0.04 vs 3.77 ± 0.05, F(1,2073) = 2.93, p = 0.0871; CA2, 2.50 ± 0.02 vs 2.50 ± 0.02, F(1,19271) < 0.01, p = 0.9867; CA3, 1.61 ± 0.02 vs 2.22 ± 0.03, F(1,17151) = 285.63, p = 1.5e – 63).
For illustration, we plotted the typical characteristics of phase precession for cases when the animal runs along (|d| < 30°) or against (|d| > 150°) the preferred rate direction (θrate) by separately fitting a regression line for the phase-position relation to passes from all fields in these two cases (Fig. 3B) and found that passes aligned to the opposite of best rate direction indeed have on average a significantly higher onset of precession than those with different directions in CA3 region but not in CA1 and CA2. The average precession slopes do not differ between the two groups of passes.
The difference in onset phases between along-θrate and against-θrate passes is also corroborated by the phase histograms (Fig. 3C), where in CA1 and CA3 the against-θrate group exhibits a significant shift to later phase as compared with along-θrate group.
Thus, at least in CA3, phase precession tends to start from a higher phase when the running direction of the rat aligns with the opposite of preferred rate direction, corroborating that phase precession exhibits directional modulations.
Directional selectivity in paired place fields
Phase precession is often considered a single-cell reflection of place cell sequences during theta (Dragoi and Buzsáki, 2006; Foster and Wilson, 2007; Feng et al., 2015; Leibold, 2020), and as such it should show up in peak lags of pair correlation functions too (Dragoi and Buzsáki, 2006; Huxter et al., 2008; Geisler et al., 2010; Schlesiger et al., 2015). In two-dimensional environments, such correlation lags have been shown to flip signs depending on the order in which a trajectory samples the place fields (Huxter et al., 2008), arguing for strong external (behavioral/sensory) drive of sequence structure. We therefore hypothesized that if phase precession reflects sequence firing, correlation lags should also be tuned to certain directions. To test our assertion, we first identified overlapping pairs of place fields (Fig. 4A) and included only spike pairs that are spaced, at most, 80 ms in time (Fig. 4B). This criterion allowed us to admit only the spike pairs that form part of a putative theta firing sequence.
Pair correlations are highly directional in CA3. A, Examples of pairs of place fields (black contour lines) and the directional tuning curves of paired spikes. Number below the contour plot indicates the amount of field overlap (calculated as 1 – Dks, where Dks is 2D Kolmogorov–Smirnov distance between two place fields). Resultant vector length (R) of the directional distribution is printed below the tuning curve. B, Raster plots of spike times in a pair of overlapping fields during traversal from field B to A. Admitted paired spikes with time difference <0.08 s between field A and B are in orange. C, R and spike-pair count for the whole populations of border and nonborder pairs in CA1, CA2, and CA3. Pairs with spike-pair counts below 40 are excluded in our statistical comparisons, indicated by the shaded region. D, Median R (top) and fraction of significantly directional pairs (bottom) by different spike-pair number thresholds. Spike pairs in all CA regions were significantly directional [CA1, 31/258 = 12.02% (Binomial test, p = 6.9e – 06); CA2, 17/181 = 9.39% (p = 0.0097); CA3, 9/88 = 10.23% (p = 0.0319)]. Comparing regions in terms of R, we found that CA3 has higher directional selectivity than CA1 [Kruskal–Wallis test; CA1 (n = 258) vs CA2 (n = 181) vs CA3 (n = 88), H(2) = 10.82, p = 0.0045; post hoc Dunn's test with Benjamini–Hochberg correction; CA1 vs CA2, p = 0.1413; CA2 vs CA3, p = 0.0030; CA1 vs CA3, p = 0.0315]. In CA1, and different from the single spike results in Figure 1, also in CA3 directionality is induced by the proximity to the border in terms of resultant vector lengths [Kruskal–Wallis test; CA1, border (n = 217) vs nonborder (n = 41), H(1) = 4.39, p = 0.0361; CA3, border (n = 56) versus nonborder (n = 32), H(1) = 9.64, p = 0.0019]. E, Fraction of all place fields (top) and border fields (bottom) by spike-pair count thresholds.
Overall, directionality results are very comparable between spike pair and single spike analysis from Figure 1. We observe significant directionality in all subregions and a higher median R in CA3 (Fig. 4D; statistics in legend). Most importantly, however, directional tuning of pairs is much higher (in terms of median R) than of single spikes (Kruskal–Wallis test for median R difference between single spikes and spike pairs at spike count threshold 40; CA1, single (n = 800) vs pair (n = 258), H(1) = 153.84, p = 2.5e – 35; CA2, single (n = 521) vs pair (n = 181), H(1) = 113.50, p = 1.7e – 26; CA3, single (n = 396) vs pair (n = 88), H(1) = 78.70, p = 7.2e – 19), and thus we conclude that spike correlations in theta sequences induce additional directionality in line with our initial hypothesis.
Again, separating field pairs into border and nonborder, we find that pair directionality in CA3 is also border sensitive (Fig. 4, statistics in legend) in contrary to single spike directionality. Particularly, border-sensitive CA3 pairs extend to high R values even for spike-pair counts >100 (Fig. 4C). These findings suggests that border sensitivity in CA3 is specifically tied to the correlation structure, whereas in CA1 it is mostly inherited from the directional firing rates.
Directionality of pair correlation
Because the pair firing rate showed region-dependent differences, we hypothesized that these differences should also transfer to spike timing correlations (Fig. 5A). Huxter et al. (2008) reported that spike correlation lags in CA1 depend on path direction and thus show a strong extrinsic (sensory/behavioral) dependence, and so we followed their approach and confirmed their main results for all three CA regions, that is, correlation lags decrease as the overlaps of the field pairs increase, and the sign of the lags flip if the path direction is reversed (Fig. 5B).
Strong extrinsic effect on correlation lag in all CA regions. A, Top row, Examples of a field pair with low overlap (left) and its cross-correlograms when the animal runs from field A to field B (middle) and in the opposite direction from B to A (right). The spike correlation lag ϕ, indicated by the distance of the closest peak to midline, generally flips sign when the direction is reversed. Bottom row, Example of a field pair with high overlap. Note that the magnitude of correlation lag is smaller than the high-overlap example. B, The higher the field overlap, the larger the spike correlation lag. Relation between spike correlation lags and the amount of field overlap for all pairs in CA1, CA2, and CA3 when the animal runs from field A to B (top) and from field B to A (bottom). Black straight line shows linear-circular regression line (Direction A to B, CA1, r(412) = –0.43, p < 1.0e – 63; CA2, r(260) = –0.39, p = 2.9e – 09; CA3, r(134) = –0.48, p = 6.6e – 08; Direction B to A, CA1, r(414) = 0.49, p < 1.0e –63; CA2, r(280) = 0.17, p = 0.0030; CA3, r(121) = 0.39, p = 1.3e – 05). Vertical and horizontal curves are the marginal distribution of spike correlation lags and field overlap, respectively.
However, a closer inspection of the correlation functions reveals that they often do not have a clear single peak, and thus we devised a new approach, taking into account the symmetries of the full correlation function. Field pairs whose activities rely on intrinsic dynamics and are insensitive to sensory stimulus should show a similar shape of the correlation function regardless of reversing the path direction. In contrast, the correlation lags of extrinsic pairs should flip sign as an effect of direction reversal (Fig. 6A,C, single pass example). Based on this principle, we were able to quantify intrinsicity and extrinsicity of a pair, using the overlap of the correlation functions of both pass directions (original and flipped, extrinsicity; original and original, intrinsicity; Fig. 6B,D; also see above, Materials and Methods). Pairs with higher extrinsicity than intrinsicity are classified as extrinsic and vice versa.
Place field correlations reveal region-specific extrinsic and intrinsic contributions. A, Top, Illustration of theta sequence from an example extrinsic pair during a single pass when the rat runs from field A to B (A → B, left). Place cell A (sky blue) fires ahead of B (purple) in each theta cycle (gray and white shaded intervals). When the trajectory direction is reversed (B → A, right), place cell B fires ahead of A. Bottom, Phase-position relations of spikes for the example pair in A, pooling from all passes in condition A → B (left) and B → A (right). B, Intrinsicity (blue text) is computed by Pearson's correlation (see above, Materials and Methods) between the cross-correlograms of two running directions (A → B and B → A, left), whereas the computation of extrinsicity (red text) takes the cross-correlogram of A → B and flipped cross-correlogram of B → A (right). C, Same as A, but an intrinsic pair. Cell B fires ahead of A even in the direction A → B. D, same as B, but with higher intrinsicity than extrinsicity. E, Top, Scatter plots of intrinsicity versus extrinsicity for all field pairs. The diagonal line represents the decision boundary; pairs above it are classified as intrinsic pairs and pairs below it as extrinsic pairs. Ratio of extrinsic pairs to intrinsic pairs (red number) and p-values (one-way χ2 test) indicate significant bias toward extrinsicity for CA1. Bottom, Histograms of differences between extrinsicity and intrinsicity. Mean μ (black bar) and p-value of the t test of mean versus zero suggest a bias toward extrinsicity for all brain regions. F, Intrinsicity versus extrinsicity for pairs with similar (difference <90°, top row) and dissimilar (>90°, bottom row) best rate angles. Detailed statistics are reported in the text.
The ratios of extrinsic to intrinsic field pairs in CA1 and CA2 are significantly different from the expected equality, whereas that of CA3 is not (Fig. 6E, top, one-way χ2 test for equal proportion of extrinsic and intrinsic pairs: CA1,
The most parsimonious explanation for intrinsic pair correlation structure would be to assume that place cells are bound into rigid sequences that play out independent of running direction. In such a scenario, place field firing should be highly directional, and pairs with similar preferred direction should reveal more of the intrinsic correlation structure than pairs with opposite preferred direction. We thus further separated pairs into those with similar (angle difference <90°) and dissimilar (>90°) best rate angles (Fig. 6F) and found that pooling over all CA regions, the similar pairs indeed have a lower extrinsic-intrinsic ratio than dissimilar pairs (similar = 203:158 = 1.28; dissimilar = 131:83 = 1.58, Fisher's exact test; p = 0.0353). Resolving for the different CA regions, the effect was significant in CA1 (similar = 108:84 = 1.29, dissimilar = 76:44 = 1.73, p = 0.0441) and CA3 (similar = 20:23 = 0.87, dissimilar = 17:7 = 2.43, p = 0.0333) but not in CA2 (similar = 75:51 = 1.47, dissimilar = 38:32 = 1.19, p = 0.0930). Thus place field pairs with similar directional tuning may contribute more to the activation of intrinsic sequence activation at least in CA1 and CA3.
Contrary to the extrinsic pairs, the intrinsic pairs are invariant to sensory inputs and maintain their firing order even when they are reversely sampled (Fig. 6C, single pass example). A possible mechanism could be that intrinsic pairs are asymmetrically connected and bias the generation of theta sequences in one direction. We thus hypothesized that the pair correlation structure reflects two contributions, extrinsic sequences that are produced by sensory inputs and flip order with the running directions and, in addition, intrinsic sequences produced by the internal connections in a fixed order.
To dissociate the intrinsic sequences from the extrinsic ones, we focused on intrinsic pairs, and for those we separated the running trajectories into two groups; either they align with (“Same” group) or are opposite (“Opposite” group) to the intrinsic directional bias that is identified via their cross-correlation signals. A directional bias of A → B would mean the correlation signal always has a peak at negative time lag even when the animal traverses from B to A (Fig. 7A). Similarly, a directional bias of B → A is determined by having both correlation peaks at positive time lag (Fig. 7B). In the Same passes, both extrinsic and intrinsic sequences could be present and hardly distinguishable as they have the same firing order. Whereas in the Opposite passes, the occurrence of extrinsic sequences is less, as indicated by the lack of reversal of the cross-correlation signal, and hence the theta sequences would be more representative of the underlying intrinsic firing order. Therefore, a comparison between Same and Opposite passes could reveal the difference between extrinsic and intrinsic sequences. More specifically, we asked whether phase precession, as a potential single-cell reflection of a theta sequence, would differ in onsets, slopes, and phase distributions between Same and Opposite running directions.
Prospective representation is revealed in CA3 intrinsic pairs. A, Top, Correlogram of an example pair with directional bias A → B (A always fires before B). Bottom, Passes that align with the directional bias A → B (left, brown, denoted as Same), and those are opposite to the directional bias (right, blue, denoted as Opposite). Origins of passes are marked as black crosses. B, An example pair with directional bias B → A. C, Fraction of passes that are precessing for Same and Opposite (Opp) conditions; p-values are derived from Fisher's exact test. D, Top, Phase-position relation for all precessing Same passes. Black solid curve marks the average of all individual linear-circular fits from each precession, and black bar marks the marginal mean phase. Bottom, Same as top but for Opposite condition instead. Linear-circular fit and marginal mean phase for Same passes from the top are also shown in dashed lines here for comparison. E, Top, Cumulative density of onsets for all precessing Same and Opposite passes. Colors as in C; p-values are derived from Watson–Williams test for differences in onset phases. Bottom, Cumulative density of slopes; p-values are derived from Kruskal–Wallis test for differences in slopes. F, Distribution of spike phases for all precessing Same and Opposite passes; p-values are derived from Watson–Williams test for differences in mean phases (vertical bars). See text for detailed statistics.
In line with our initial hypothesis, we found that precession is more likely to occur when passes align with the directional bias (Fig. 7C). The Same condition has a higher fraction of precession than Opposite condition in all CA regions, but only reaching significance in CA1 and CA2 (Fisher's exact test; CA1, p = 0.001; CA2, p = 0.0284; CA3, p = 0.0550). Considering only the precessing passes, there is no difference in precession slopes between Same and Opposite [Fig. 7E; bottom; Kruskal–Wallis test; CA1, Same (n = 236) vs Opposite (n = 174), H(1) = 0.69, p = 0.4048; CA2, Same (n = 213) vs Opposite (n = 164), H(1) = 0.15, p = 0.7015; CA3, Same (n = 61) vs Opposite (n = 39), H(1) = 2.52, p = 0.1126], but the phase onset of precession exhibits clear regional differences. CA3 has higher phase onset when the passes are opposite to the directional bias (Fig. 7E; top, Watson–Williams test; CA3, F(1,98) = 5.24, p = 0.0243), whereas such difference is not significant in CA1 (F(1,408) = 2.54, p = 0.1115) and CA2 (F(1,375) = 0.35, p = 0.5552). Consistently, pooling the spikes from all precessing passes, the spike phase distribution from Opposite passes shows a significant shift to later phases as compared with Same passes only in CA3 (Fig. 7F; Watson–Williams test; CA3, F(1,2002) = 18.76, p = 1.6e – 05) but not in CA1 (F(1,6017) = 1.86, p = 0.1730) and CA2 (F(1,6272) = 0.48, p = 0.4890).
According to our hypothesis that extrinsic and intrinsic sequences are played out in parallel and that intrinsic pairs result from a suppression of extrinsic sequences in one direction, the association between intrinsic sequences and later spike phases should be more strongly visible in intrinsic pairs than in extrinsic pairs. We thus also inspected the phase distributions in extrinsic pairs. Because by definition, extrinsic pairs do not have a directional bias, we resorted to compare the spike phase distributions from all precessing passes between extrinsic and intrinsic pairs, regardless of their pass directions. The upward shift of spike phases in intrinsic sequences (if they include more intrinsic pairs), should then show up as a higher marginal spike phase than that of extrinsic pairs. Statistical analysis confirms this prediction and shows that precessing passes from intrinsic pairs indeed exhibit higher spike phases than extrinsic pairs in CA2 and CA3, but not in CA1 (Watson–Williams test, mean ± SEM in radians; CA1, Ex 2.07 ± 0.018 vs In 2.03 ± 0.026, F(1,15816) = 0.82, p = 0.366; CA2, Ex 2.30 ± 0.020 vs In 2.59 ± 0.032, F(1,9227) = 45.28, p = 1.8e – 11; CA3, Ex 1.20 ± 0.035 vs In 1.44 ± 0.044, F(1,4232) = 12.68, p = 0.0002).
Our results support our hypothesis that intrinsic pairs display two types of sequences, extrinsic and intrinsic in the Same direction and predominantly intrinsic sequences in the Opposite direction. Trajectories Opposite to the directional bias of spike pair lags display later spike phases and higher onsets.
It might seem counterintuitive that when the intrinsic pairs are traversed in the opposite direction, sequences play out in a backward order, and yet there is phase precession. We therefore propose a possible explanation for the coexistence of intrinsic sequence and phase precession that is in line with our analysis results by assuming that intrinsic pairs consist of a leading cell and an enslaved cell. The enslaved cell only fires with a delay after the leading cell was active. Consequently, their correlation structure is fixed in both directions. As the leading cell undergoes phase precession, the enslaved cell fires after the precessing spikes of the leading cell, and, hence, also precesses with a phase shift (Fig. 8A).
Co-occurrence of phase precession and intrinsic sequences. A, Schematic illustration of how intrinsic sequences arise from enslaved spiking. The enslaved cells only fire at a fixed delay after the leading cell, causing a fixed temporal order even if the animal runs in the opposite (Opp) direction (bottom row). This relation predicts that enslaved cells mostly fire in the overlap region of the two place fields (shaded area). B, Phase precession of leading (cyan dots) and enslaved (dark purple) cells plotted for Same, Opp, within, and outside overlap region, respectively, from left to right. Data were generated artificially for illustration purposes. Marginal distribution of spike phase (solid curve) and its circular mean (solid black bar) are shown on the right of each plot. Light violet dots and dashed lines represent an alternative dataset of enslaved spikes with larger phase shifts from the leading spikes. Note that the phase means of Opp passes and in nonoverlap regions are higher than Same and overlap, respectively. C, Schematics of theta sequences produced by Normal pairs without enslaved spikes. Note that the firing order is reversed when the pair is traversed in the opposite direction and both cells fire symmetrically in the nonoverlap region. D, Phase precession in normal pairs. Note that there is no difference in marginal phase means across all conditions. E, Phase distribution of precession spikes in Overlap and Nonoverlap regions for intrinsic (top) and extrinsic (bottom) pairs in CA1, CA2, and CA3; p-values are derived from Watson–Williams test comparing the circular difference of phases between nonoverlapping regions and overlaps. CA3 shows a larger effect in the phase difference as compared with CA1 and CA2. See text for detailed statistics.
The first prediction of the enslavement hypothesis is that the enslaved cell only fires within the overlap region between the two fields. As a result, the Opposite condition would have a higher spike phase than the Same condition because both leading and enslaved spikes are present in the prospective cycles at the beginning of the field (Fig. 8B, left). Whereas in Normal pairs, without enslaved spikes (Fig. 8C), there would be no phase difference between Same and Opposite (Fig. 8D, left). The prediction is compatible with our experimental findings in Figure 7, where the marginal phase mean and onset in the Opposite condition are higher than in the Same condition, indicating a possible contribution of enslaved spikes to the intrinsic structure of CA3.
Another assumption in the enslavement hypothesis illustrated in Figure 8 is that the precession of leading cells in the Opposite condition would start at a lower onset and last for a shorter number of cycles than in Same because of the diminished forward recurrent connection when the pair is traversed in the opposite order. It would lead to the second prediction that the mean phase in the region where place fields are nonoverlapping should be higher than within the field overlap (Fig. 8B, far right). On the contrary, the “Normal” pairs would show no phase difference between overlap and nonoverlap regions as phase precession is symmetrical in both directions (Fig. 8D, far right).
The prediction of higher overlap phase in nonoverlap regions from the enslavement hypothesis can be tested by comparing the phase distributions between overlap and nonoverlap regions in the experimental data (Fig. 8E, top row). We found that in CA3 intrinsic pairs, nonoverlap spiking indeed has a higher mean phase than in the overlap (mean ± SEM in radians; Nonoverlap 2.02 ± 0.07 vs Overlap 1.35 ± 0.04, Watson–Williams test, F(1,2745) = 50.69, p = 1.4e – 12). A similar trend is also observed in CA1 intrinsic pairs, but the difference is smaller than in CA3 (Nonoverlap 2.41 ± 0.03 vs Overlap 2.18 ± 0.02, F(1,9905) = 26.03, p = 3.4e – 07). These results support that in CA1 and CA3, the intrinsic structure is generated by enslaved spikes. As a control, we also inspected the nonoverlap-overlap phase differences in the extrinsic populations (Fig. 8E, bottom row). In CA3, the phase difference has decreased as compared with its intrinsic counterpart (CA3 extrinsic, Nonoverlap 1.57 ± 0.05 vs Overlap 1.14 ± 0.04, F(1,3813) = 33.53, p = 7.6e – 09), whereas in CA1, there is even a higher mean phase in the Overlap region (CA1 extrinsic, Nonoverlap 1.97 ± 0.03 vs Overlap 2.06 ± 0.02, F(1,17 523) = 7.48, p = 0.0062), indicating that enslaved spiking is less involved in the extrinsic pairs.
Intrinsic sequences and phase precession are thus not contradicting each other, if we assume that one field exhibits only enslaved spikes. The existence of enslaved spikes can explain both higher spike phase in the opposite traversal direction as well as in the nonoverlap area, which are both in agreement with the data.
Relation between pair correlation and phase precession
Following up the hypothesis from Figure 8 that intrinsic pairs are resulting, at least partly, from enslaved spikes, we suggested that both phase precession and pair correlations contribute different aspects to theta sequences, as proposed previously (Middleton and McHugh, 2016). We thus finally asked whether we can identify a direct relation between the intrinsicity/extrinsicity property and the directionality of phase precession in the data to further corroborate our enslavement hypothesis. To this end, we distinguished trajectories that are either parallel (<90°) or opposite (>90°) to the best rate angles of the fields of a pair with parallel best rate directions (Fig. 9A) and included only passes with phase precession in both fields. Because the late spike phases appear to be associated with the passes opposite to the best rate angles in single fields (Fig. 2), as well as the intrinsicity in CA3 (Fig. 7), we hypothesized that the association should also transfer to field pairs. Our analyses show that indeed when the trajectories oppose both of the best rate angles, there is a higher proportion of intrinsic pairs contributing to phase precession (Fig. 9B) and a higher marginal spike phase of co-occurring phase precession in CA3 (Fig. 9C) than in the case when the trajectory runs in parallel to both best rate angles. These directional differences are not observed in CA1 pairs. CA2, surprisingly, shows the opposite trend that higher spike phases are more associated with the best rate direction, but there is no difference between the contributions of extrinsic and intrinsic pairs. The findings on CA3 corroborate that both correlation structure and phase precession are direction dependent in CA3, where the intrinsic field pairs seem to lack strong extrinsic drive in one running direction and thus would exhibit enslaved spikes at later spike phases when the running direction opposes their best rate direction. In extrinsic pairs the extrinsic drive seems strong in all running directions to overrule the intrinsic structure and allows cells to fire at earlier phases.
Field pairs in CA3 exhibit higher intrinsicity and fire at later spike phases in the direction opposite to the best rate angles. A, Schematic illustration of how trajectories are divided into two groups, Both-parallel (left), if trajectories across a field pair deviate <90∘ from their best rate angles, and Both-opposite (right), if trajectories deviate >90∘ from their best rate angles. B, Numbers of extrinsic and intrinsic pairs exhibiting phase precession in Both-parallel and Both-opposite cases. In CA3, there is a significantly higher contribution from intrinsic pairs in Both-opposite case than in Both-parallel case; p-values are derived from Fisher's exact test. C, Phase distributions of phase precessing passes. The Both-opposite case in CA3 shows a significantly later spike phase than the Both-parallel case. Vertical bars mark the circular means; p-values are derived from Watson–Williams test (Both-parallel vs Both-opposite, mean ± SEM in radians; CA1, 2.15 ± 0.0252 vs 2.18 ± 0.0369, F(1,7917) = 0.53, p = 0.4654; CA2, 2.57 ± 0.0321 vs 2.22 ± 0.0369, F(1,4638) = 40.31, p = 2.4e – 10; CA3, 1.26 ± 0.0561 vs 1.67 ± 0.0609, F(1,1841) = 18.55, p = 1.7e – 05).
Computational model as control
To further corroborate that the observed directionality effects in place field activity are not just an artifact of running trajectories or data analysis, we applied our analysis to simulations of a spiking model that does not include directional information as input (Romani and Tsodyks, 2015) but can be fed with the animal trajectories of our experiments (see above, Materials and Methods for simulation details). In brief, the model neurons integrate place-specific inputs, theta-periodic inputs, and symmetric recurrent connections with short-term synaptic depression. The recurrent connections give rise to omnidirectional phase precession but, because of their symmetry, do not impose preferred intrinsic sequential activity. Also, because the place-specific inputs are not directionally modulated, we reasoned that any directionality and any intrinsicity we find in the model must be artificial. Indeed, as expected, directionality analyses on the single field firing rates showed that the fraction of significantly directional simulated place fields does not exceed chance level of 5% (61/1024 = 5.96%, p = 0.0936, Binomial test; Fig. 10A). Similarly, by pooling over all fields again, we did not find a significant correlation of phase precession per spike with the distance to preferred rate direction (Fig. 10B; Spearman's correlation coefficient rs(1386) = –0.27, p = 0.0710). A further inspection on the distribution of preferred precession angles on a single field basis did not reveal a significant π shift with respect to the preferred rate angles (Fig. 10D; V-test against π; for all passes, V(383) = –26.21, p = 0.9709; Fig. 10E; for low-spike passes, V(174) = –30.61, p = 0.9995). After again separating the passes into groups of against-θrate (n = 105) and along-θrate (n = 361), with >150° and <30° angular difference in radians from the best rate angles, respectively (Fig. 10F), we consistently found no significant difference of their average precession curves (Kruskal–Wallis test for difference in slopes; H(1) = 0.09, p = 0.7591; Watson–Williams test for difference in onsets; F(1464) = 0.01, p = 0.9415), and there was also no significant difference in the means of spike phases between groups against-θrate and along-θrate (Watson–Williams test; F(1,7639) = 2.78, p = 0.0953; Fig. 10G). Thus, as expected by model design, properties of phase precession in the model generally do not show a significant dependence on the heading directions of the animal.
Simulation results using the model by Romani and Tsodyks (2015). A, Directionality of simulated single fields. Median R (left) and fraction of significantly (Sig. Frac.) directional place fields (middle) by spike count thresholds for all (solid line), border (dotted), and nonborder (dashed) fields. Right, the fraction of all border and nonborder place fields by spike count thresholds. Frac., Fraction. B, Distribution of precession incidences (dark blue) and spikes (orange) as a function of |d(θpass,θrate)|, the difference between pass direction and best rate direction of the place field. Ratio of blue and orange line in green shows no significant trend (Spearman's correlation, rs(1386) = –0.27, p = 0.0710). C, Scatter plot for the relation between best precession angle θprecess and rate angle θrate of all place fields reveals no obvious structures. D, Distribution of θprecess directions of all precessing passes corrected to θrate shows no significant 180° difference of two best angles. E, Distribution of θprecess directions of low-spike passes. F, Average precession slopes (phase-position curves) from precession samples against (green) and along θrate (blue); ps and po are derived from Kruskal–Wallis test for slope difference and Watson–Williams test for onset difference respectively. G, Distribution of spike phases for precession samples against and along θrate. Vertical bars denote the circular means of the distributions; p-value is derived from Watson–Williams test comparing the difference of two circular means. Detailed statistics are reported in the text.
A virtue of the investigated Romani and Tsodyks (2015) model is that it contains recurrent synaptic connections, and we asked to what extent they might explain intrinsic correlation structure observed in the data. We thus repeated directionality analysis at the level of field pairs on the simulated data. Similar to the single field results, the fraction of significantly directional pairs never exceed chance level in the model (169/3197 = 5.29%, p = 0.2394, Binomial test; Fig. 11A); however, the dependence of correlation lags on place field distance can be well reproduced (Fig. 11B), demonstrating that phase precession induces direction-dependent sequential firing of model place fields.
Directionality of pair correlations in model simulations. A, Median R (left) and fraction of significantly directional pairs (middle) by spike count thresholds for all (solid line), border (dotted), and nonborder (dashed) pairs. Right, Fraction of all, border, and nonborder field pairs. B, Relation of spike correlation lags and field overlap of all pairs for direction A → B (left) and B → A (right). Linear-circular regression line is in black (A → B, r(2659) = –0.38, p <1.0e –90; B → A, r(2663) = 0.35, p < 1.0e – 90). Vertical and horizontal curves show the marginal distributions. C, Intrinsicity versus extrinsicity for all pairs (one-way χ2 test of equal extrinsic-intrinsic ratio). D, Density of differences between extrinsicity and intrinsicity (mean µ, black bar; p-value from Student's t test of mean vs 0) show a significant trend toward extrinsicity.
We next computed extrinsicity and intrinsicity of the model correlations and found a fraction of extrinsic pairs similar to CA1 (Fig. 11C; Ex:In = 1313:851 = 1.54; one-way χ2 test for equal proportion of extrinsic and intrinsic pairs;
To conclude, although a network model with asymmetric weight couplings can reproduce the results of phase precession and correlation lags in a two-dimensional environment including the extrinsic bias prominently observed in CA1 ensembles, further extension is still required to account for directional selectivity in place field networks and their effect on phase precession, specifically on the onset phase, and the less extrinsic dynamics of CA3.
Discussion
We showed that theta-scale timing of place cell activity in all regions of the hippocampus proper exhibits directional modulations. This first applies to phase precession, which is most prominent (per spike) in the direction opposite to the one with the largest firing rate, and this effect is more prevalent and acute in CA3 than in CA1 and CA2. In addition, best precession tuning in CA3 is associated with higher onset phases. Second, by inspecting directionality of spike pairs from two overlapping place fields, we showed that pair correlation exhibits higher directionality than single spikes from individual place fields. Furthermore, using cross-correlation analysis, we demonstrated that CA1 pair correlations are better explained by external sensory inputs induced by overt movement than CA3 pairs. This suggests that CA3 place field correlations are more strongly intrinsically determined. In addition, a closer inspection into the intrinsic pairs revealed that the intrinsic sequences, which are invariant to the sampling order of place fields, are also associated with higher onset and spike phases in CA3 phase precession. Finally, we used the model by Romani and Tsodyks (2015) and showed that the interaction between directionality and theta activity (described above, Results) cannot be solely explained by an omnidirectional generation mechanism of theta sequences. Nonetheless, the model shares a similar dependence on extrinsic information, thus being more analogous to CA1 rather than CA3.
The directionality of place field firing in 2 d open environments has a long history of debate. Although classical studies on the influence of landmark cues (Kubie and Muller, 1991), as well as early modeling work (Touretzky and Redish, 1996), clearly acknowledge the availability of directional information to the place cell system, there was disagreement about the degree to which this information becomes overt in place field activity (Muller et al., 1994; Anderson et al., 2006; Acharya et al., 2016). In the present study, we have found 10−20% of significantly directional place fields, which is similar to previous reports (Markus et al., 1995; Acharya et al., 2016; Mankin et al., 2019), although directionality seems to strongly depend on the behavioral setting. Markus et al. (1995) found 80% of place cells with significant directionality in eight-arm mazes and 20% in an open circular platform. Acharya et al. (2016) also showed that the significant fraction of heading-direction modulated cells could vary as a function of the width of visual cue on the wall as more place cells were head-direction modulated near the border.
Most interestingly, and consistent with these previous reports, our analysis revealed that directionality of CA1 in familiar environments is heavily induced by proximity to the border, even at a high spike count threshold where the sampling bias is relatively small. However, this border effect is less prominent in simultaneously recorded CA3 place cells. One possible explanation is that there is no border-related directionality in novel environments but develops via learning processes, which is more difficult and, hence, slower in CA3 because of its recurrence. This idea has been suggested in a previous modeling study from Brunel and Trullier (1998), which demonstrated that place cells become less directionally selective through synaptic plasticity when they are visited in all directions (as in the nonborder case) as compared with the border case when fields are visited only in a subset of directions. Also a second model by Navratilova et al. (2012) suggested that the firing rate of place cells is initially independent of directions but develops to be directionally selective through experience and synaptic plasticity. Assuming that the same learning mechanism further differentiates between border and nonborder directionality, the border difference would thus develop slower in CA3 because of its larger degree of recurrent connectivity. Both models support the idea that the differential directionality between border and nonborder cells in familiar environments could be explained via plasticity induced by early experiences in novel environments. Our finding, however, that borders induce strong directionality of CA3 spike pairs (Fig. 4), questions the assumption of these models that directionality is learned but rather supports the idea that CA3 pair directionality reflects the directional imbalance of extrinsic and intrinsic activity (Fig. 8).
The main conjecture of the present work was that inherent network dynamics could induce directionality of place cell spike timing. Through multiple lines of analysis, we repeatedly found that in our dataset, CA3 exhibits more hints to intrinsic dynamics than CA1. This was, for example, suggested by their theta-scale correlations being more invariant to the trajectory. These findings are very much in accordance with anatomy. Pyramidal cells in CA3 project extensive recurrent collaterals in rodents and primates (Amaral et al., 1984; Anderson et al., 2006), whereas such recurrence is not so much obvious in the CA1 region (but see Deuchars and Thomson, 1996, for rodents). However, a further possibility to explain CA3 bias toward intrinsic dynamics could be its unique position in the mammalian hippocampal-entorhinal circuitry. Sensory information from entorhinal cortex layer II is projected to CA3 not only via the direct perforant pathway but also via mossy fiber pathway of dentate gyrus (DG), contributing to an extra layer of internal processing through the granule and mossy cell loop. Indeed, CA3 activity in rats with DG lesions showed reduced prospective firing in an eight-arm maze (Sasaki et al., 2018), which is strikingly consistent with our observed direction dependence of theta onset phases. The anatomic consistency of the enthorinal-hippocampal formation across mammals, as well as consistent reports on hippocampal phase precession in several mammalian species including Chiroptera (Eliav et al., 2018) and primates (Qasim et al., 2021), let us speculate that the bias toward stronger intrinsic sequences in CA3 is a common feature within this animal class as a whole and specifically also expected to be seen for primates and humans despite their less precise spatial firing patterns (Ekstrom et al., 2003).
The directional selectivity of phase precession provides widespread support of our initial hypothesis that directional information should also express itself in hippocampal theta sequences, yet the anatomic foundations of this directional modulation can only be speculated on so far. One possible explanation would be to assume that it is the direct entorhinal inputs to CA3 that induce the directional dependence of firing rate, which would be consistent with directional drive from postsubiculum and entorhinal head direction and conjunctive cells (Sargolini et al., 2006). According to Figures 2D and 7, this input would be reduced in passes opposite to best rate directions. Consequently, intrinsic CA3 activity would then mostly be explained by the indirect DG mossy fiber pathway, which was shown to induce prospective out-of-field spiking (Sasaki et al., 2018), potentially indicating enslaved spikes. In the best rate direction, according to our hypothesis, strong direct entorhinal input should induce earlier theta spike phases in addition to the DG-induced theta sequences. Conversely, against the best rate direction, DG-induced theta sequences would still play out invariantly despite the weaker entorhinal directional drive. In such a scenario, we would predict that intrinsic pairs would dominate phase precession properties over those extrinsically induced. A direct test of these predictions is provided in Figure 9. Although the outcomes are consistent with our predictions, the conclusiveness is limited because of small sample sizes, particularly in CA3. Yet the above hypotheses could be tested even more directly by acute differential suppression of entorhinal and DG pathways, which we would predict to selectively alter phase precession in the different running directions.
Intrinsic hippocampal correlations are incompatible with the interpretation of the hippocampus as a pure spatial map but rather imply the existence of prestructure (Dragoi and Tonegawa, 2011). The degree to which intrinsic activity is expressed in familiar environments is, however, relatively small and clearly visible only in CA3, thus arguing for a relatively small role of intrinsic structure in supporting spatial navigation in familiar environments. The situation in novel environments, however, might be different with a stronger effect of prestructure on organizing hippocampal representations. This is in line with the finding that theta sequences only develop with time (Feng et al., 2015) as initially some pairs may predominantly fire in the wrong intrinsic order. The idea—contrary to the models of Brunel and Trullier (1998) and Navratilova et al. (2012)—that directional invariance develops in an experience-dependent way on top of prestructured sequences has important implications for theories of spatial memory formation because it would favor the notion of sensory integration into existing temporal structure (cf. Buzsáki and Tingley, 2018; Leibold, 2020).
Footnotes
This work was supported by the German Research Association under Grant LE2250/13-1 and the National Institute of Mental Health under Grant R01MH119179. We thank Geoffrey W. Diehl for help with data transfer.
The authors declare no competing financial interests.
- Correspondence should be addressed to Christian Leibold at christian.leibold{at}biologie.uni-freiburg.de
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
