Uncovering spatial representations from spatiotemporal patterns of rodent hippocampal field potentials.

Spatiotemporal patterns of large-scale spiking and field potentials of the rodent hippocampus encode spatial representations during maze runs, immobility, and sleep. Here, we show that multisite hippocampal field potential amplitude at ultra-high-frequency band (FPAuhf), a generalized form of multiunit activity, provides not only a fast and reliable reconstruction of the rodent's position when awake, but also a readout of replay content during sharp-wave ripples. This FPAuhf feature may serve as a robust real-time decoding strategy from large-scale recordings in closed-loop experiments. Furthermore, we develop unsupervised learning approaches to extract low-dimensional spatiotemporal FPAuhf features during run and ripple periods and to infer latent dynamical structures from lower-rank FPAuhf features. We also develop an optical flow-based method to identify propagating spatiotemporal LFP patterns from multisite array recordings, which can be used as a decoding application. Finally, we develop a prospective decoding strategy to predict an animal's future decision in goal-directed navigation.

Introduction

Cutting-edge, large-scale, high-density electrode arrays have enabled us to record spatially distributed neural activity within local or multiple brain circuits (Berényi et al., 2014; Jun et al., 2017; Chung et al., 2019), yet it remains challenging to reliably uncover spatial representations of large-scale neural activity. The hippocampus contributes to a wide range of brain functions responsible for spatial and episodic memories, learning, and planning (Pfeiffer and Foster, 2013). Population spike activities from hippocampal place cell assemblies provide a readout of the rat's spatial location (Zhang et al., 1998; Davidson et al., 2009; Kloosterman et al., 2014; Grosmark and Buzsáki, 2016; Hu et al., 2018; Ciliberti et al., 2018). However, direct use of spike information for neural decoding remains challenging due to practical issues of spike sorting and unit instability, let alone the prohibitive computational cost in the era of big data (Hu et al., 2018; Ciliberti et al., 2018; Rey et al., 2015; Rossant et al., 2016; Gridchyn et al., 2020). In contrast, hippocampal field potentials consist of collective local synaptic potentials around the recording site, serving as an alternative information source for spatial representation (Buzsáki et al., 2012; Agarwal et al., 2014). In addition to the animal studies, the high-density electrocorticography (ECoG) grid has allowed us to access large-scale field potentials from the human brain (Khodagholy et al., 2015; Chang, 2015; Zhang and Jacobs, 2015). Therefore, developing effective and scalable statistical methods for uncovering neural representations of these signals during memory tasks or sleep remains a central goal in computational neuroscience.

Hippocampal sharp-wave ripples (SPW-Rs) are important hallmarks for memory reactivations during immobility and non-rapid eye movement (NREM) sleep (Roumis and Frank, 2015; Buzsáki, 2015; Fernández-Ruiz et al., 2019). Decoding the content of hippocampal replays during ripple events can help dissect circuit mechanisms of memory consolidation, planning, and decision making (Pfeiffer and Foster, 2013; Davidson et al., 2009; Foster, 2017; Grosmark and Buzsáki, 2016; Ólafsdóttir et al., 2018). Evidence has suggested that hippocampal local field potentials (LFPs) may encode neuronal ensemble activations during SPW-Rs (Taxidis et al., 2015), but it remains unknown how the information embedded in field potentials is related to spatial representations and how these spatially distributed hippocampal field potential representations relate to ensemble spike representations at different brain states.

To address these questions and challenges, we systematically investigate spatiotemporal patterns derived from large-scale multisite recorded rodent hippocampal field potentials and their spatial representations during maze running, memory replay decoding, and spatial decision-making. In addition to previously identified demodulated theta (amplitude/phase) features (Agarwal et al., 2014), we propose a new set of independent and complementary spatiotemporal features derived from continuous hippocampal field potentials and develop innovative supervised and unsupervised learning methods while validating them in position decoding during a maze run (on rats and mice) and during hippocampal ripple events (on rats). We demonstrated their excellent decoding performance compared with the current state-of-the-art decoding analysis methods. Together, these robust spatiotemporal features and simple decoding methods provide not only direct and robust solutions to position-decoding and decision-prediction problems in real-time brain-machine interface (BMI) applications, but also means for visualizing high-dimensional neural data or finding latent structures of these high-dimensional data in the absence of behavioral measure.

Results

Spatially distributed hippocampal field potentials encode the position in a maze run

Rats and mice were implanted with large-scale or high-density silicon probes (varying from 64 to 512 channels; see Table S1) in the dorsal hippocampus or hippocampi while freely foraging in a circular track, linear track, T maze or open-field environment (see Figure S1). We processed the extracellular raw voltage traces of hippocampal recordings and extracted independent neural signals at different frequency bands, including unclustered and clustered spikes (Figure 1A). During run epochs, we found that both demodulated theta (4–12 Hz) activity—or in short dLFPθ (amplitude or phase or both) (Agarwal et al., 2014) and the instantaneous amplitude at ultra-high-frequency (>300 Hz) band, or FPAuhf —a proxy of the unclustered multiunit activity (Stark and Abeles, 2007; Bansal et al., 2012; Smith et al., 2013) recorded from multiple recording sites could be used to reliably decode an animal's position in the unseen held-out data (Figures 1B and 1C). Therefore, simple linear decoding methods (see STAR Methods) based on these continuous and un-thresholded activities provided comparable or improved decoding accuracy compared with advanced Bayesian decoding methods using sorted-spike or unsorted spikes (Kloosterman et al., 2014; Hu et al., 2018; Ciliberti et al., 2018; see STAR Methods). Notably, the same decoding method derived from different features may produce different position-decoding distribution statistics (Figure S2A). Examining the joint probability distribution function of the decoding error and position or behavioral sampling (occupancy) can provide a guideline in how to design the experiments and data collection (Figure S2B). Augmenting multiple independent and complementary features, such as combining clustered spike, dLFPθ, or FPAuhf with their own activity history (Figure 1) or combining FPAuhf with dLFPθ (Figure S1B), further improved the 10-fold cross-validated decoding accuracy (Figure 1C; two-sample Kolmogorov-Smirnov test, maximal absolute difference (D) and p value: D = 134 cm and p = 1.6 × 10−4, D = 136 cm and p = 5.1 × 10−5, D = 142 cm and p = 1.5 × 10−2 for three respective panels). We systematically assessed the decoding performance using instantaneous amplitude or phase features at various frequency bands and temporal bin sizes and found that >300 Hz amplitude was most effective (Figure S3A). In temporal windows of <300 ms the various features and their combinations were comparable, whereas at longer time intervals decoding errors increased more rapidly for LFPθ than for spike-containing combinations (Figure S3B). Based on the same features, different decoding methods (Table S2) produced varying degrees of decoding accuracy (Figures S3C, S3D, and S4), which depended on the data dimensionality, electrode implant, channel layout, and behavioral sampling of the animals. Unless stated otherwise, we used 100 ms temporal bin size, 2 cm position bin size, and an optimal linear estimation (OLE) decoding method for all supervised decoding analyses, and all decoding error statistics were reported based on 10-fold cross-validation.

Figure 1

Spatiotemporal features of the hippocampus encode the rat’s position during maze running

(A) Illustration of processing raw voltage traces of hippocampal recordings into four different signals: unsorted spikes (green), clustered spikes (red), filtered LFP at theta band- (blue), and high-pass-filtered signal (>300 Hz, light green).

(B) Decoded rat's trajectory on a circular track (dataset 1 in the key resources table) based on four features derived from (A): unsorted spikes (orange), clustered spikes (red), demodulated LFP theta (referred to as dLFPθ, blue), >300 Hz LFP amplitude (referred to as FPAuhf, green). The rat's position is marked by the black line the and decoded position shown as scatterplot.

(C) 2D histograms (left) and error cumulative distribution function (CDF; right) curves of 10-fold cross-validated decoding errors during maze run (numbers indicate the median error). Augmenting the covariate with their respective history in the OLE decoder (second to fourth rows) further improved the decoding accuracy in each case by ∼10% (dotted line in the CDF plot, numbers in brackets indicate the median decoding error of held-out data).

(D) Illustration of position decoding in a 2D open-field environment (dataset 5 in the key resources table) based on the same decoding strategy derived from clustered spikes, dLFPθ, FPAuhf, and joint (dLFPθ + FPAuhf) features.

(E) CDFs of decoding error derived from the four different signals shown in (D). Numbers indicate the median decoding error of held-out data.

(F) Decomposition of a broadband raw voltage trace (black) into the sum of a spike-free component (magenta) and a spike-only component (orange). To investigate the effect of spike waveform and amplitude (dataset 1 in the key resources table), we normalized individual spike waveforms by their respective peak amplitudes (brown). Ticks in the top represent identified spikes from clustered units and their averaged spike waveforms.

(G) Down-sampled (to 1,250 Hz) and high-pass filtered (>300 Hz) version of the four signals in (F).

(H) CDFs of decoding error based on amplitude features derived from the four signals shown in (G). Numbers indicate the median decoding error of held-out data (dataset 1 in the key resources table).

See also Figures S1–S4; Tables S1 and S2; Videos S1.

Spatial coverage of the hippocampal place fields is crucial to the decoding accuracy in a larger environment. In the open field, we found that 64-channel joint dLFPθ + FPAuhf features yielded better decoding accuracy than clustered spike-based decoding (dataset 5 in the key resources table; mean ± SD of bootstrapped median decoding error for held-out data: 10.81 ± 0.25 cm for spikes and 8.85 ± 0.21 cm for dLFPθ + FPAuhf; rank-sum test, p = 1.95 × 10−21; Figures 1D and 1E; see also Video S1).

Video S1. Demonstration of Bayesian filter (likelihood plus a temporal prior) decoding in an open-field environment (dataset 5) based on clustered spikes, LFPθ, MUA, joint (LFPθ + MUA) features, and combining all features, related to Figure 1

Video frame: 3.3 Hz.

It is noted that the >300 Hz filtered signal also contains spike activity (Zanos et al., 2011; Buzsáki et al., 2012; Waldert et al., 2013). To remove the contamination of spike activity, we decomposed the broadband signal into spike-free (by removing detected spikes) and spike-only components (by retaining only detected spikes; Figure 1F). Following down-sampling (1,250 Hz) and high-pass filtering (>300 Hz) operations (Figure 1G) and repeating decoding analyses, we showed that the spike-only component produced the best decoding accuracy (Figure 1H), suggesting that the position decoding contribution of FPAuhf was mainly derived from the spiking-related activity, although decoding the animal's position does not require a preprocessed clustering step. To further investigate the effect of spike waveform and amplitude, we normalized individual spike waveforms by their respective peak amplitudes and repeated the same decoding analysis (Figure 1F). These results show that, in our FPAuhf decoding strategy, the most important information came from spike amplitude feature; this is not surprising considering our high-pass filtering and Hilbert transforming operation (Figure 1H). It is noteworthy to point out that, even after the removal of the above-threshold spike waveforms, the remaining spike-free high-frequency component can still yield decent decoding performance (Figure 1H, median error 9.1 cm derived from the spike-free component and 5.3 cm from the raw LFP)—which might be due to the contamination of undetected spike activities on the spike-free component. The unclustered spikes (discrete-valued) and FPAuhf (continuous-valued) can be viewed as different forms of multiunit activity, but the information between them is overlapped but not identical. As FPAuhf contains both above-threshold and sub-threshold unit activities from pyramidal cells and interneurons, it therefore represents richer information beyond the sorted or detected spiking activity. Put together, these results not only confirmed the previous findings (Hu et al., 2018; Agarwal et al., 2014), but also indicate the advantages of the new proposed or combined features based on a simple linear decoder.

Since spiking activity is derived mainly from the somatic layer, we further examined the impact of anatomical location of high-density silicon probe electrodes on decoding performance (Figure 2; dataset 3 in the key resources table). As expected, recording sites in the CA1 and CA3 pyramidal layers contributed most to overall FPAuhf decoding accuracy, since these two layers are known to have higher cell density. Furthermore, the decoding accuracies from all test conditions were significantly better than the chance level (derived from the randomly shuffled data, Figure 2B).

Figure 2

Position decoding using recording sites in different hippocampal layers

(A) Sharp-wave ripple-triggered current source density analysis revealed distinct hippocampal LFP sinks and sources in different anatomical layers (LFP traces from different layer [Pyr., Rad., Lm., DG, and CA3p] indicated by different colors, 256 channels, 8 shanks, dataset 3 in the key resources table).

(B) Demodulated dLFPθ (top) and FPAuhf (bottom) median position decoding error using recording sites in different layers (blue) compared with median position decoding error (error bar: SD) derived from 1,000 random selections of matched number of recording sites (red) and 1,000 random shuffle of their position data (orange, decoding baseline; see STAR Methods) (Pyr, stratum pyramidale; Rad, stratum radiatum; Lm, stratum lacunosum moleculare; DG, dental gyrus; CA3p, CA3 pyramidal layer). The numbers in parentheses indicate the channel count within that layer/region. CA1 layers with changing layer-by-layer theta-phase shift gave best position decoding by dLFPθ decoding strategy, while FPAuhf was most effective in CA1 and CA3 pyramidal layer. p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001; two-sample t test.

Hippocampal FPAuhf features decode replay events

During immobility and NREM sleep, putative hippocampal memory replay events occur during SPW-Rs, where population firing bursts coincide with a high ripple band amplitude (Figure 3A). A central task of statistical analysis is to identify the content of these replays and assess their significance in the absence of behavioral ground truth. However, population spike activity may be sparse, especially during NREM sleep, bringing a challenge in decoding analysis. Therefore, a field potential-based features may be very appealing in replay decoding analysis. As an illustration, here we used three different decoding approaches (clustered spikes with traditional Bayesian replay decoder, clustered spikes with the OLE decoder, and FPAuhf with the OLE decoder; see STAR Methods for details) and found that all three approaches could reliably reconstruct the replay trajectories (Figures 3B and S5A). To evaluate the significance of the reconstructed trajectories, we further computed the distance correlation and the shuffled statistics (see STAR Methods). We found that FPAuhf and clustered spike-based decoding strategies produced slightly different replay trajectories (Figure 3C) and Z score or Monte Carlo p value statistics (Figures 3D, S5B, and S5C). This is interesting since all three different approaches could produce excellent decoding results in the maze run state but slightly different replay trajectories in the NREM sleep state. In addition, the derived Z score statistics were also different since spike-based and LFP-based decoding approaches used different data-shuffle strategies. To gain more insight into how varying degrees of spiking information influence the replay significance result, we analyzed the total number of firing units in every candidate event and further compared their replay significances. The results showed that the significance of replay events was dependent on the total number of firing units during that event (Figure S5D). This suggests that not only the selection of feature and decoder but also the data-shuffle method may affect the ripple replay analysis and assessment. Despite its demonstration purpose in this preliminary replay analysis (due to limited data recordings), our proposed approach suggests an alternative means for future study of hippocampal replay analyses.

Figure 3

Spatial distributed hippocampal LFP patterns encode rat position during ripple events

(A) Illustration of a broadband hippocampal LFP (black) during a candidate replay event. The candidate event was determined by a joint criterion (see STAR Methods) using the ripple band amplitude (blue) and total unclustered spike count (red, FPAuhf). The red and blue horizontal dotted lines indicate their respective thresholds.

(B) Three representative example trajectories of virtual maze positions during hippocampal SPW-Rs events, derived from a traditional Bayesian decoder using clustered spikes (top), from an OLE using clustered spikes (middle), and from an OLE using FPAuhf (bottom). Dark pixels show high probability values. Red circles, green rectangles, and green triangles mark the estimated replay trajectories using different approaches.

(C) 2D histogram comparisons between trajectories from three estimators shown in (B).

(D) Comparison of Z-scored statistics of distance correlation (see STAR Methods) for ripple replay candidate events (n = 727 ripples during waking and NREM sleep). Decoding methods (OLE or Bayesian) and features (FPAuhf or clustered spikes) are labeled in each paired comparison. Dashed lines mark the Z score value of 2.33 (p = 0.01). Four distinct colors are used to illustrate the event membership within divided regions. The top right region indicates the significant replay events detected by both methods.

See also Figure S5.

Real-time linear decoding strategy for large-scale hippocampal recordings

Our results directly point to an efficient hippocampal FPAuhf- or dLFPθ-based decoding strategy with two key advantages: low sampling rate requirement and easy implementation, as well as the stability and longevity of signals (in contrast to unreliable spike detection and sorting of unstable units). First, we validated the robustness of the FPAuhf- or dLFPθ-based linear decoding strategy. In multiple consecutive recording sessions (separated by days in between) from the same animal and the same environment (dataset 6 in the key resources table), we trained the linear decoder using a one-run session and tested it on remaining run sessions. Remarkably, the FPAuhf and dLFPθ decoding strategies produced outstanding performance across several consecutive run sessions/days (Figure 4A). Furthermore, we assessed the average cross-session decoding performance when training and testing sessions were separated by different intervals, and found a V-shaped performance drop with increasing time intervals (Figure 4B). When the intervals between training and testing session were short, the FPAuhf features produced better decoding performance, yet the performance of LFPθ features decayed slower. Since FPAuhf and dLFPθ features were complementary, combining them further improved the decoding accuracy, especially in the case of chronic recordings where the number and quality of recorded units are unstable.

Figure 4

Robustness and scalability of linear decoding strategies based on different features

(A) Cross-session decoding error matrices summarize within-session (cross-validated) and between-session assessment using dLFPθ, FPAuhf, and dLFPθ + FPAuhf features in five consecutive run sessions (dataset 6 in the key resources table, mouse 1). The number in the matrix shows the median decoding error. Sessions (labeled by A, B, C, D, and E) were separated by ∼24 h.

(B) Average statistics of the kth (−4 ≤ k ≤ 4) diagonal of cross-session decoding error matrices in (A). k = 0 implies within-session cross-validated decoding; k < 0 (lower diagonal) implies training on historical sessions and testing the subsequent sessions in time; k > 0 (upper diagonal) implies the reverse order.

(C) OLE decoding using FPAuhf (green, 20 ms; blue, 100 ms bin size) with respect to the number of channels, N. In our computer simulations, we assumed a linear mapping with dimensionality of N × 145 (290 cm track, 2 cm position bin size). Signals were down-sampled (to 1,250 Hz) and high-pass filtered (>300 Hz) before computation. Shaded areas represent SD from 1,000 realizations.

(D) Illustration of online detection and assessment of hippocampal memory replay during NREM sleep (for details, see STAR Methods). In brief, the computational procedure is as follows. First, from the raw extracellular voltage trace (first row) of a pre-selected electrode channel, the onset of a candidate event was identified using a threshold criterion based on the ripple band amplitude (second row). Next, a “spatial trajectory” was reconstructed using an OLE decoder derived from hippocampal FPAuhf (20 ms bin size; third row). Once the duration of candidate event became more than three bins, the significance of the decoded spatial trajectory of 60 ms or longer was assessed based on online shuffling statistics (fourth row), reporting the Monte Carlo p value (red curve). Furthermore, an accumulative score (blue curve) was continuously updated if p < 0.05 (red horizontal dashed line): . A significant replay was determined when the accumulative score was above a threshold (blue horizontal dashed line). The accumulative score was set to 0 at the detection onset (marked by arrow) and reset to 0 when the cumulative score threshold was reached. The computation time for online evaluation is shown in the last row. As shown, all computation times per bin were below a sampling interval of 20 ms, supporting the real-time analysis.

See STAR Methods for detailed descriptions. See also Table S3.

Furthermore, we ran computer simulations to test the scalability of our decoding algorithm. Since the training and decoding of our algorithms are fully automatic, we could update our linear decoder at any time. For the sake of simplicity, we choose to use FPAuhf features alone to achieve best decoding accuracy and to reduce computational complexity. Running an optimized C/C++ code on a multi-core desktop computer, the readout could accommodate a real-time speed at a scale up to hundreds of thousands of (>100,000) recording channels during both wake and sleep (Figures 4C; 100 ms temporal bin size for wake and 20 ms for ripple replay). Next, it is important to assess the statistical significance of online decoded replay events during NREM sleep. We adapted a real-time spike-based replay assessment strategy (Hu et al., 2018) to accommodate FPAuhf decoding (see STAR Methods). By presetting the pseudorandom shuffle operations (n = 2,000 shuffles), we could accommodate a real-time speed (<20 ms) for both decoding and significance assessment for 128 channels recoding with 1,250 Hz sampling rate in a 1D environment (Figure 4D), with comparable results derived from offline significance assessment (Table S3). In addition, all these real-time scalability tests were performed on a regular desktop computer with CPU only; however, the state-of-the-art clusterless online decoding algorithm (Hu et al., 2018) has to employ GPU hardware and special software that may increase the complexity and compromise the reliability in practical applications.

Put together, in comparison with the previous clusterless decoding algorithm (Kloosterman et al., 2014; Hu et al., 2018; Ciliberti et al., 2018), the low sampling rate requirement, simple feature extraction (high-pass filtering and Hilbert transform), efficiency decoding method (OLE, linear mapping/matrix multiplication), and the stability and longevity of signals make our approach more suitable for online closed-loop neural manipulation experiments. The scalable computation also enables the feasibility of many hippocampal replay studies based on large-scale neural recordings, especially for closed-loop neuroscience experiments (Chen and Pesaran, 2021).

Unsupervised learning reveals consistent representations between low-rank structures of hippocampal FPAuhf with spikes of neuronal ensembles

Establishing a mapping between the animal's position and neural activity requires training samples during maze runs. However, from the perspective of an internal brain observer, it is important for downstream structures of the hippocampus to quickly infer spatial representations without direct behavioral measures. Unsupervised learning and topological representations provide an alternative perspective for hippocampal neural codes (Chen et al., 2014; Linderman et al., 2016; Maboudi et al., 2018). In parallel to supervised decoding strategies, here we employed two latent variable models and unsupervised learning methods, nonnegative matrix factorization (NMF) and reconstruction-independent component analysis (RICA), to extract lower-rank features from multisite recorded hippocampal FPAuhf (Figures 5A and 5B; see STAR Methods for details). Position-averaging on the lower-rank feature matrix derived from NMF or RICA revealed localized structures in the latent state space.

Figure 5

Unsupervised learning methods for extracting multichannel FPAuhf features

(A) Illustration of NMF (top, blue) and RICA (bottom, red) methods for extracting lower-rank features from FPAuhf during the maze run (dataset 1 in the key resources table). The greyscale color bar illustrates the strength of feature amplitude. After sorting, the position-average patterns (from the rows of NMF's matrix H or RICA's matrix Z) yielded localized position-tuned features covering the entire track (similar to “place fields”). The colored boxes (matrices H and Z, derived from NMF and RICA, respectively) were used to derive the position-tuned features.

(B) Same as (A), during ripple events. The virtual position was decoded from clustered spikes.

(C) A Bayesian hidden Markov model (HMM) was used for inferring latent state sequences {S1(t)} and {S2(t)}, based on the spike ensembles and FPAuhf features, respectively. The sorted ensemble spikes are assumed to be Poisson distributed. The lower-rank FPAuhf features were obtained from the NMF's matrix H, followed by a feature resampling procedure (see STAR Methods). The matrix frame colors (green, blue, red) represent the source of the feature matrix (clustered spike count, H, Z).

(D) During maze running, the Bayesian HMM analysis revealed a strong one-to-one correspondence map between the animal's position and the inferred latent states {S1(t), S2(t)}, as well as between {S1(t)} and {S2(t)} (upon sorting). The grayscale color indicates the normalized count or occupancy statistics (with dark color representing higher occupancy). These consistency results indicated that the most valued information in FPAuhf is spike activities (top row: based on NMF; bottom row: based on RICA). Label color (S1 in green and S2 in blue or red) represents the latent state sequence derived from clustered spikes and lower-rank features (H or Z).

(E) During ripples in the absence of behavior measures, the Bayesian HMM analysis also revealed a correspondence map between {S1(t)} and {S2(t)}, which were derived from the clustered spikes and FPAuhf features (top row: based on NMF; bottom row: based on RICA), respectively. Label color same as (D).

See also Figures S6–S8.

From the lower-rank FPAuhf features (derived from NMF or RICA, followed by feature resampling; see STAR Methods and Figure S6), we further trained an unsupervised Bayesian hidden Markov model (Chen et al., 2014; Linderman et al., 2016) and inferred the latent state trajectories (Figure 5C). During run, the lower-rank FPAuhf-inferred state trajectories matched well with the animal's position as well as spike-inferred state trajectories (Figure 5D). During ripples, we also observed similar correspondence between two latent state trajectories derived independently from lower-rank FPAuhf features and clustered spikes (Figure 5E). In both maze run and ripple events, the preprocessed lower-rank feature extraction was critical for learning the latent state trajectories; our analyses showed that direct use of LFP or FPAuhf features alone yield poor performance.

Furthermore, projecting high-dimensional FPAuhf features onto a 2D embedding space revealed coherent structures in agreement with the animal's position or spatial topology of the environment (Figure S7). In addition, in order to visualize the individual extracted NMF or RICA features in the recording channel space, we used a deflation approach to project the features back to the LFP channel space to account for their contributions (Figure S8; see STAR Methods for details). Since the sorted NMF/RICA features correspond to a “position code,” the channel space maps reveal the information where the individual “place code” is located in the channel space. Together, these results confirmed that unsupervised learning methods can uncover neural representations of large-scale hippocampal field potentials without a priori measurement of the animal's behavior, which has been a challenge for inferring memory replays of hippocampal non-spatial representations or hippocampal-neocortical representations during sleep (Chen et al., 2016; Chen and Wilson, 2017; Allen et al., 2016).

Independent, parallel, and complementary hippocampal spatiotemporal patterns for spatial representation

Propagating waves in the rat hippocampus have been reported at the theta and ripple frequency bands (Agarwal et al., 2014; Lubenov and Siapas, 2009; Patel et al., 2012, 2013). In many cases, it is difficult to identify the spatiotemporal patterns of electrode array recordings because of the lack of regular patterns. To fully exploit spatially distributed features of high-density electrophysiological recordings, here we introduce a method, motivated from optical flow used in computer vision, to extract spatiotemporal “wave” patterns (see STAR Methods for details). Specifically, we investigate how these spatiotemporal patterns across multisite hippocampal recordings were coordinated in space and time while animals ran in the maze. We mapped the band-pass-filtered hippocampal field potentials (at 4–12 Hz or >300 Hz band) onto the 2D electrode space in time to obtain a video sequence, and then employed an optical flow estimation approach to compute the vector field between two consecutive image frames (see STAR Methods and Figure 6A). We found that the spatiotemporal optical flow revealed position-tuned wave patterns in the FPAuhf- or other LFP-derived features, suggesting that they can serve as a proxy for spatial readout. Therefore, the optical flow method can help reveal the intrinsic spatiotemporal patterns of band-limited signals (beyond theta and ultra-high-frequency bands) during run or other task behaviors. In principle, this method can also be used for other high-density array recordings, such as ECoG, EEG, and MEG, despite the difference in their spatial coverage.

Figure 6

Spatiotemporal patterns of multisite recorded hippocampal LFP encode run trajectories

(A) Illustration of optical flow estimation and visualization by 2D vector field. High-density silicon probe recordings (8 shank, 32 sites each shank) in the same hippocampus in the transversal (T) and longitudinal (L) axes of a rat are shown separately (dataset 3 in the key resources table). Each 2D vector was visualized as an arrow (characterized by size and direction). (i) Background color in the raw map represents the voltage of filtered LFP traces. (ii) Background was color coded according to the arrow direction. The color map was chosen to represent a continuous circular angle variable (0°–360°).

(B) Optical flow estimated from band-pass-filtered (4–12 Hz) 2 × 256-site recorded LFP signals show position-dependent patterns on a T maze (dataset 3 in the key resources table). The vector field patterns show trial-by-trial consistency at the same position on the maze, suggesting that low-frequency hippocampal spatiotemporal activity carries spatial information of an animal's position. See Video S2 for the demonstration of clear propagating wave patterns.

(C) Schematic of extracting parallel and independent spatiotemporal features from hippocampal field potentials at theta (4–12 Hz) and ultra-high (>300 Hz) frequency bands. For notation simplicity, features are labeled as 1, 2, 3, and 4.

(D) Comparison of median position decoding error during the maze run derived from four different spatiotemporal features in (C) and their feature combinations (dataset 3 in the key resources table). Error bars show the bootstrapped SD. See also Figure S9 and Table S4; Video S2 and S3.

In addition to exploratory data visualization, we further investigated whether the high-dimensional optical flow features were useful for decoding the animals’ positions. Remarkably, these spatiotemporal patterns estimated directly from band-pass-filtered (4–12 Hz) hippocampal field potentials provided sufficient information for reconstructing the animals' positions (referred to as LFPθ-flow in Figure 6B; see also STAR Methods, Figure S9, and Video S2). This result suggested that low-frequency hippocampal spatiotemporal activity (without demodulation) carries spatial information of the animal's position. In fact, we discovered that a set of hippocampal field potential-derived spatiotemporal patterns (Figure 6C) contained complementary information for position decoding (Figure 6D for dataset 3; see also Figure S9C for dataset 1). The variability in contributions of individual spatiotemporal patterns in different datasets might be ascribed to different spatial sampling as well as different layout in the implanted probes (Figure S1C). Together, these hippocampal spatiotemporal features or their combinations may provide a rich repertoire for position decoding across different brain states (Table S4).

Video S2. Demonstration of location-tuned spatiotemporal patterns in 512-channel rat hippocampal LFP signals filtered at 4–12 Hz (dataset 3), related to Figure 6

The animal's actual position (gray) and decoded position (blue) are shown. The 2D vector field illustrates the optical flow patterns across two Buzsaki256 probes. The background is color coded according to the arrow direction. Video frame rate: 10 Hz.

Video S3. Demonstration of location-tuned spatiotemporal patterns in 120-channel rat hippocampal MUA features while running in a circular maze (dataset 1), related to Figure 6

The animal's actual position (gray) and decoded position (blue) are shown. The 2D vector field illustrates the optical flow patterns across two Buzsaki64SPL probes (6 × 10 layout). The background shows the MUA amplitude (warm colors represent high value). Video frame rate: 10 Hz.

It is noteworthy to point out that, while using hippocampal theta-band LFP demodulation, it requires that all recorded sites should share a strong common oscillation pattern (such as theta), such that the useful information can be extracted based on a global demodulation operation (Hu et al., 2018). Otherwise, there is no global pattern available for demodulation. Unfortunately, this requirement is rarely satisfied in LFP recordings across many brain states and many frequencies. In contrast, the proposed optical flow method focuses on a local spatiotemporal pattern and does not have any prerequisite (expect the regularity of the layout of recoding sites supported by large-scale high-density electrode arrays). Therefore, the optical flow method has a great potential in many exploratory high-density LFP or ECoG studies, such as position/item decoding, decision-prediction, mental stimulation, and virtual reality experiments.

Prospective decoding strategy for predicting goal-directed navigation behavior

Hippocampal place cells exhibit trajectory-dependent or prospective/retrospective memory coding in spatial navigation (Frank et al., 2000; Ferbinteanu and Shapiro, 2003; Ji and Wilson, 2008; Wood et al., 1999; Catanese et al., 2014; Grieves et al., 2016). To demonstrate our proposed decoding strategies, we further investigated if a prospective decoding strategy based on hippocampal field potentials could predict an animal's decision in a spatial alternating task (dataset 6 in the key resources table). We adapted the standard decoding strategy and attempted to predict the animal's prospective L/R arm position, while the animal ran in the central arm toward the choice point, with forward time lag ranging from 0.5 s to 1.5 s (5–15 temporal bins, 100 ms bin size; Figure 7A). Remarkably, both hippocampal spikes and joint (dLFPθ + FPAuhf) features (Figure 7B; see also Figure S10) carried predictive representations of upcoming maze locations in the L/R arm. We found that the prediction accuracy in prospective coding improved as animals moved closer to the choice point, where the predictability was far beyond the chance level (see shuffled statistics in Figures S10A and S10B). In two tested animals, the joint features achieved a high (86%–97%) prediction peak accuracy in correct trials (n = 4 sessions for each mouse; 50% chance level; see Figures S10A and S10B for incorrect trials). Overall, we found that the FPAuhf feature produced the overall best results, followed by the dLFPθ feature. Since individual hippocampal place cells can remap their place fields based on the behavioral context, it is likely that an ensemble of such trajectory-dependent place cells carries sufficient information for MUA-based L/R arm prediction. The exact mechanism of demodulated theta contributing to this phenomenon remains to be determined. According to Buzsáki et al. (2012), local LFP “snapshot” features (such as the demodulated theta) can reflect unique constellations of cell assemblies responsible for the discharge of neuron. Therefore, time-evolving constellation of the LFP map derived from high-density LFP recordings may provide information of the evolving cell assemblies. Together, joint features derived from both LFP amplitude and phase information can be used to predict an animal's goal-directed navigation behavior.

Figure 7

LFP, FPAuhf, and spike-based features predict the mouse's decision-making

(A) Schematic of prospective decoding, with a time lag of 0.5–1.5 s (5–15 temporal bins, 100 ms bin size) before the choice point in the T maze.

(B) Comparison of cross-validated prediction accuracy between clustered spike-based and joint (dLFPθ + FPAuhf) decoding strategies (dataset 6 in the key resources table). Mean and SD results were pooled over multiple sessions. Left panel: mouse 1 (n = 4 sessions, n = 120 correct trials per session); right panel: mouse 2 (n = 4 sessions, n = 120 correct trials per session).

(C) Top panel: cross-validated support vector machine classification accuracy (for mouse 1). The peak correct classification rate from cross-validation was ∼90% and the AUROC statistic was ∼0.94. Bottom panel: boxplot statistics of the animal's actual location at specific temporal bins across all tested trials.

See also Figure S10.

In another independent test, we further trained a linear support vector machine classifier to predict the L/R decision (for correct trials only) by pooling dLFPθ and FPAuhf features at five consecutive temporal bins before the choice point (Figure 7C; see STAR Methods for details). The cross-validation classification accuracy and trend were also comparable with Figure 7B in two tested animals (peak accuracy, correct trials: 86%–98% for spikes, 77%–95% for joint [FPAuhf + dLFPθ] features; incorrect trials: 66%–94% for spikes, 64%–89% for joint features). The gradual decay in classification accuracy away from the choice point was partially due to the increased trial variability in the animal's actual position because of the variability in run speed. Notably, the degraded accuracy trends were slightly different between forward and backward directions, when triggering temporal bins from either the end (“backward”) or start (“forward”) position of the central arm (Figures S10C and S10D).

Discussion

Rodent hippocampal clustered spikes and spatially distributed theta waves are known to encode an animal's position. Reading out spatial representations from unclustered spikes or field potentials can provides a mechanistic insight into hippocampal circuits (Kloosterman et al., 2014; Hu et al., 2018; Ciliberti et al., 2018; Agarwal et al., 2014). Here, we further demonstrate that unthresholded and unclustered hippocampal spatiotemporal patterns at different frequency bands contain rich spatial representations in wake and sleep. Based on large-scale high-density electrophysical recordings of the rodent hippocampus, we propose a set of independent spatiotemporal features derived from hippocampal field potentials; we further develop efficient statistical decoding methods to achieve ultrafast readout of spatial representations. These features are complementary; the decoding performance derived from each feature depends heavily on the layout of the silicon probe and the implanted location of hippocampus layer. Based on the implanted electrode location and spatial coverage at the hippocampus, we can assess the differential contributions of spatial representation from different layers of hippocampal subregions. We have tested the proposed decoding methods and compared them with several standard methods across six datasets (and a few more, but the results not shown due to the similarity in conclusion). The decoding accuracy comparison between LFP-based with spiking-based methods varied, which may depend on multiple factors, such as behavioral sampling, unit yields, the number of electrode channels, electrode array layout, and the implant location. In general, electrode arrays with large spatial spread will likely favor demodulated theta- or traveling wave-type decoding methods. Electrode implants in the layers with the densest CA1 or CA3 units will likely favor spike- or multiunit-based decoding methods. The demodulated theta feature (dLFPθ) and optical flow-theta feature are somewhat correlated, as combining these two features did not significantly improve the decoding accuracy. Although these two features aim to capture the relative phase variation across electrodes, the optical flow feature is spatially local and requires no global demodulation as in the originally proposed demodulated theta information (Agarwal et al., 2014).

The instantaneous readout of multisite spatiotemporal patterns of hippocampal field potentials provides a readout of integration of spatiotemporal information along the septotemporal axis of the hippocampus, forming the traveling waves (Lubenov and Siapas, 2009; Patel et al. 2012, 2013). Such traveling waves can also propagate the coded information to the downstream structures of the hippocampus. Our results on FPAuhf decoding were also consistent with previous findings that unsorted hippocampal ensemble spikes can reliably encode a rodent's position (Kloosterman et al., 2014; Hu et al., 2018; Deng et al., 2015). However, our current method is simpler and different from the previous methods in that the multiunit activity in FPAuhf contains spiking information from both pyramidal cells and interneurons, whereas clusterless decoding methods focus on unsorted spikes from hippocampal pyramidal cells (Kloosterman et al., 2014; Hu et al., 2018). Recently, the information encoded at the wide frequency band of hippocampal field potentials has been independently demonstrated based on deep learning (Frey et al., 2021). High-frequency field oscillations are known to represent multiunit activity or the spectral leakage of spiking activity (Scheffer-Teixeira et al., 2013). The similar ultra-high-frequency band (300–1,000 Hz) of unsorted spiking activity can also predict movement in BMIs (Nason et al., 2020). Based on such ultra-high-frequency band information, our analysis also demonstrated that the prospective decoding strategy is robust and reliable to predict an animal's decision-making in goal-directed navigation. The most computationally efficient FPAuhf decoding strategy is particularly useful for content-based, real-time closed-loop circuit manipulation in rodents with a chronic implant, which may prove critical for causal circuit dissection of hippocampal-cortical/subcortical coordination in decision-making (Singer et al., 2013; Schmidt et al., 2019). In addition, the FPAuhf feature can potentially overcome the issues of gradually degraded unit loss during chronic electrophysiological recordings. The degraded decoding accuracy based on LFP phase or amplitude information (Figure 4B) might be due to the drift of electrodes due to brain movement, or due to gradual remapping and slow drift of hippocampal representation (Frank et al., 2004; Dong et al., 2021; Rule et al., 2019). In another independent investigation (Tu et al., 2020), we applied the OLE decoding method to decode a rodent's position based on calcium-imaging recordings, and demonstrated that population decoding accuracy degraded across days even when a large number (>500) of CA1 units were used. Easy acquisition, efficient sampling, processing, and transmission of band-limited signals make it ideal for wireless-empowered implanted devices.

A central task in computational neuroscience is to link neural representations with an animal's behavior or internal cognitive state. Therefore, development of unsupervised learning methods is important for discovering intrinsic latent structures of high-dimensional spatiotemporal neural data, especially in the absence of behavioral measures (Linderman et al., 2016; Cunningham and Yu, 2014; Townsend and Gong, 2018; Chaudhuri et al., 2019). Since high-density LFP recordings produce a high degree of input correlation between neighboring channels, low-rank feature extraction or dimensionality reduction (e.g., ICA, NMF, and embedding) can help visualization of exploratory data and subsequent decoding analysis. Extending the previous report (Agarwal et al., 2014), our unsupervised learning results present another demonstration to uncover spatial correlations directly from rodent hippocampal field potentials at the ultra-high-frequency band, especially for replay analysis during sharp-wave ripples. Notably, our approach is rather general and may be easily extended to other neural circuits as well as other behavioral contexts. Our preliminary analyses have shown that the FPAuhf feature can be useful for identifying structures from multisite LFP recordings during quiet wakefulness and sleep replay. However, because of the difference in data-shuffle operation, it remains to be determined whether spike- or LFP-based features are more reliable for identifying significant memory replay events. Systematic comparisons or improved detection methods are further required in future investigations.

Hippocampal spatiotemporal patterns are often displayed in the form of traveling waves (Lubenov and Siapas, 2009; Patel et al. 2012, 2013) and sequential structures (“hippocampal sequences”) in space and time during various behavioral states (Buzsáki and Tingley, 2018). In addition to ripple sequences, hippocampal theta sequences have been known to encode look-ahead trajectories or current goals in planning (Foster and Wilson, 2007; Gupta et al., 2012; Wikenheiser and Redish, 2015). Furthermore, internally generated sequences may predict an animal's choices in a memory task or guide navigation when external spatial cues are reduced (Pastalkova et al., 2008; Wang et al., 2015; Villette et al., 2015). To date, rodent hippocampal memory replays have been widely studied in NREM sleep. During REM sleep, rodent hippocampal LFP theta oscillations are pronounced, yet their spatiotemporal patterns are less well understood (Louie and Wilson, 2001; Zielinski et al., 2021). Investigation of REM sleep hippocampal theta waves may reveal new biological insight into learning and memory, as well as dreams. The spatiotemporal patterns of multisite hippocampal field potentials and decoding methods proposed here may provide new opportunities to investigate these REM sleep-associated spatiotemporal patterns.

Although we only investigated the rodent hippocampal circuit in this study, our “place” decoding analysis and methods can be readily applied to other rodent neocortical areas that encode spatial information, such as the entorhinal cortex, retrosplenial cortex, parietal cortex, primary somatosensory cortex, and primary and secondary visual cortices (Hafting et al., 2005; Whitlock et al., 2008; Mao et al., 2017; Ji and Wilson, 2007; Haggerty and Ji, 2015; Hu et al., 2018; Long and Zhang, 2021; Long et al., 2021). Moreover, it has been shown that these sensory “place cells” are also modulated with local cortical theta rhythms during run behaviors and preserved phase precession (Long and Zhang, 2021; Long et al., 2021). To date, experimental evidence of spatial modulated responses across many neocortical areas beyond the traditional hippocampal formation has been accumulating. It remains to be determined whether our proposed LFP-based decoding methodologies can be applied to readout the spatial representations from high-density cortical recordings. Examination of the coordinated representations of large-scale hippocampal-neocortical spatiotemporal activity during various behavioral states would help dissect the mechanisms of memory, learning, planning, and decision-making (Chung et al., 2019).

Finally, the superior representational power of hippocampal field potentials points to potential applications for predicting human hippocampal memory task outcomes based on ECoG recordings (Zhang and Jacobs, 2015) or for investigating human memory replays based on high-density MEG or EEG recordings (Kurth-Nelson et al., 2016; Liu et al., 2019; Huang et al., 2018). Although non-invasive MEG/EEG signals do not carry ultra-high-frequency spiking information, MEG/EEG source localization can still reveal interesting spatiotemporal patterns of the hippocampal-entorhinal network or medial temporal lobe related to spatial memories or other task variables. Our proposed optical flow methods based on the derived band-pass-filtered spatiotemporal features may play a subserving role in exploratory data analysis.

Limitations of the study

In this study, we have not thoroughly investigated the application of our proposed methods in NREM/REM sleep-associated hippocampal recordings due to the lack of sleep recordings with both high-unit yields and high channel counts. Systematic verification of our methods on high-density hippocampal electrophysiological recordings from independent experimental labs will be the subject of future investigation.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Zhe S. Chen (zhe.chen@nyulangone.org).

Materials availability

This study did not generate new materials.

Experimental model and subject details

All experimental studies were performed in accordance with the National Institutes of Health (NIH) Guide for the Care and Use of Laboratory Animals to ensure minimal animal use and discomfort, and were approved by the New York University School of Medicine (NYUSOM) Institutional Animal Care and Use Committee (IACUC). Long-Evans rats (n=3) and mice (n=2) were used in the behavior tasks, and were maintained on a 12:12-h light-dark cycle and housed individually with access to food and water.

Electrophysiological recordings

Datasets 1 and 2: Circular track and linear track (rat 1)

Male Long-Evans rats were bilaterally implanted with two 6-shank silicon probes (128 channels in total) parallel to the septo-temporal axis of the dorsal hippocampus. The silicon probe is a custom Buzsaki64SPL probe (NeuroNexus). This 6-shank silicon probe had 10 sites at each shank, and there were 4 extra channels spaced every 1.25 mm starting 1.25 mm from the tip of Shank 4 (6×10+4=64 channels). All sites were vertically staggered along the shank with 20 μm spacing between sites (Figure S1C). We selected one rat (‘Achilles’) in the analysis. The recording session consisted of a long (∼4 hour) pre-RUN sleep epoch in a familiar room, followed by a RUN epoch (∼45 minutes) in a novel circular maze (1 m diameter, Dataset 1 in the key resources table) or a linear track (1.6 m long, Dataset 2 in the key resources table). After the RUN epoch the animal was transferred back to its home cage in the familiar room where another long (∼4 hour) post-RUN sleep was recorded. Details of experimental protocols and data have been published (Grosmark and Buzsáki, 2016). The electrophysiological data are publicly available (https://crcns.org/data-sets/hc/hc-11/).

Datasets 3 and 4: T-maze and linear track (rat 2)

Two 256-channel (512 channels in total) custom-made silicon electrodes were implanted to the right hippocampus of the rat. The silicon probe is a custom Buzsaki256 probe with 32×8 array layout (Figure S1C). The two probes were perpendicularly aligned to each other in order to record along both the septotemporal and subiculo-fimbrial axis of the hippocampus. In the T-maze, rat was trained to perform a delayed alteration task, in which the animal had to choose either the left or the right arm at the decision point. After returning to the start area, the rat was confined for 10 seconds. In the following trial, the rat had to choose the opposite direction and would obtain water reward with a correct choice. In the linear track, animal simply foraged back and forth to collect water reward. Details of experimental protocols have been published (Berényi et al., 2014), and data are available at http://buzsakilab.com/wp/datasets/.

Dataset 5: Open field (rat 3)

Male Long-Evans rats foraged and chased randomly dispersed drops of water or food on an elevated square platform (120 ×120 cm2). The silicon probe consists of two 32-channel (4×8) array (Figure S1C). The probe was implanted in the rat’s left and right dorsal hippocampi. Data are available at https://crcns.org/data-sets/hc/hc-3.

Dataset 6: Circular T-maze (mouse 1 and mouse 2)

Male mice were trained to perform a continuous spatial alternation task in a circular T-maze. The silicon probe consists of a 64-channel (4×16) poly2 layout (Figure S1C). Each animal’s recordings consisted of multiple sessions during 4-6 consecutive days, and each session lasted 50-60 min. In each session, animal was able to run 120-160 trials, with a minimum of 120 correct trials. The probe position was not changed across sessions in order to assure the stability of LFP channels. Spikes were sorted separately for each session, yielding varying number of units (Table S1).

Method details

Feature extraction methods

Spike sorting

Extracellular representations of action potentials were extracted from recorded broadband (0.3 Hz-10 kHz) signals followed by threshold-based spike detection algorithm. Individual spikes were automatically clustered using the KlustaKwik algorithm (http://klustakwik.sourceforge.net/) or Kilosort algorithm (https://github.com/cortex-lab/KiloSort). The clusters were manually refined by discarding multiunit clusters showing lack of clear refractories in the autocorrelogram, and groups with unstable firing patterns over time.

Theta-band LFP demodulation

We applied a complex Morlet wavelet (with central frequency at 8 Hz) filter to the broadband LFP through (non-causal) zero-phase filtering (MATLAB function ‘filtfilt’). Then, we applied a complex-valued principal component analysis (PCA) (Agarwal et al., 2014) to the filtered multi-channel LFP signal, where produced the first principal component (PC) that reflects the oscillatory component common to all channels. At each channel, from the theta-band filtered analytic LFP signal, we applied demodulation using the following equation to obtain theta-demodulated LFP signals where , denotes the band-pass filtered signals at the theta band, denotes the phase of the first principal component (PC) of the complex-valued, filtered multi-channel LFP signal within the theta frequency band (Agarwal et al., 2014).

Field potential amplitude (FPAuhf) at ultra-high frequency band

The multi-channel LFP signals were first high-pass filtered by a 4th order Butterworth filter at 300 Hz cutoff frequency. Then we applied the Hilbert transform to the filtered signal, and computed the amplitude feature using the absolute value of analytical signal. This FPAuhf can be viewed as a continuous version of the proxy for the multi-unit activity, which contains both sub-threshold and supra-threshold spiking activity from both pyramidal neurons and interneurons.

Decoding analysis

Several spike-based or LFP-based decoding methods were investigated (see Table S2 and Figure S2C). In all decoding analyses, spikes or LFP features were Z-scored across features (neurons or channels), respectively.

Cross-validation

In all supervised decoding analysis, we used 10-fold cross-validation to separate training and testing data. For a specific session, we concatenated all data together and used the first 90% data to train the decoder, and decoding performance was only assessed on remaining 10% data. We repeated this process for another 9 times to make sure all data were used as testing data in different repetition.

Optimal linear estimation (OLE) decoding

We used the standard OLE method previously described before (Agarwal et al., 2014). Let denote the decoded variable, and denote the C-dimensional observation (LFP feature or Spike). In one-dimensional environment, we mapped the linearized position (with length L) via K equally spaced von Mises functions characterized by a circular variable : for k=1,…,K. The OLE tries to solve a least-squared (LS) estimation problem:Where denotes the unknown parameters, and denotes the vectors of spike count or multichannel LFP features at time t, and denotes a K-dimensional vector that expands the position in the basis. Unless stated otherwise, we used K=75, and a temporal bin size of 100 ms (maze run) or 20 ms (ripple events). In all decoding analyses, all LFP features were first averaged within each time bin and then Z-scored across all channels. We rescaled the OLE output between 0 and 1 via a softmax operation.

In addition, we considered incorporating the prior history of LFP activity (such as ) as additional covariance in the linear regression estimator. However, this was at the cost of doubling computational complexity and potential overfitting.

In the open field, we mapped the two-dimensional position by replacing K von Mises basis functions with a set of two-dimensional Gaussian tiling (Agarwal et al., 2014):Where the vectors represents the animal’s 2D position, and denotes the (diagonal or isotropic) covariance, and denotes the mean vectors of K Gaussians that tile the field. We used K=144 and temporal bin size of 300 ms for the open field environment.

OLE decoding by imposing a sparsity constraint

In the OLE method, the traditional linear least-squared (LS) regression method is subject to overfitting in the presence of high-dimensional features. To improve generalization and assist variable, we used a sparse Bayesian linear regression method based on variational Bayes (VB) and automatic relevance determination (ARD) (Drugowitsch, 2013). Briefly, we assumed that the scalar output conditional on a multivariate input followed a Gaussian likelihood model, and the prior of regression coefficients and the precision variable has a conjugate normal inverse-gamma distribution. The goal of VB is to maximize the lower bound of the marginal likelihood (evidence) to obtain the variational posterior distribution of the regression coefficients. Because of the use of conjugate priors in Gaussian likelihood, the inference of variational posteriors can be analytically derived (Drugowitsch, 2013). The introduction of ARD assigned a separate shrinkage prior to each element of the regression coefficients (to impose the sparsity), which was in turn adjusted by a hyper-prior. Software implementation can be found in https://github.com/DrugowitschLab/VBLinLogit. Since the output was mutually independent, we extended the multi-input single-output (MISO) regression problem to a multi-input multi-output (MIMO) regression problem. Finally, the VB inference produced the posterior mean and posterior variance of the individual parameter in . For variable selection, the parameter with a small mean and a small variance would be discarded. See Figure S4 for an illustration of encoding and decoding application with VB-ARD.

LFP feature likelihood-based decoding

We used a Gaussian likelihood-based decoder (i.e., noninformative prior) based on LFP featureswhere denotes the variance, and the mean is represented by a sum of one or two-dimensional basis functions

The Gaussian log-likelihood function, denoted as , is written aswhere . The unknown parameters were estimated from the maximum likelihood estimation.

Bayesian spike decoding

For decoding with sorted spikes, one method we used is traditional replay reconstruction algorithm (Zhang et al., 1998). According to Bayes’ rule, the probability of animal at position given the spiking activities in a short window is:

By assuming a Poisson firing for each unit, the likelihood function written as:where donates the average firing rate of unit at position , which can be inferred from the maze run data (ripple event analysis, Figure 3) or training data (position decoding analysis, Figure 1D). Δ is the temporal bin size, and is the spike count of unit within that window during the ripple event. Combining equations above yields posterior probability:

where is a normalization factor that can be calculated by the normalization condition . denote the prior probability, which can be assumed as a uniform distribution (Figure 3). In the case of uniform prior, the Bayesian decoder reduces to a likelihood-based decoder. Alternatively, can depend on the previously decoded position (Gaussian filter, Figure 1D) (Agarwal et al., 2014):

Finally, the reconstructed position is derived from the maximum a posteriori (MAP) estimate:

Ensemble unsorted spike decoding

For performance comparison, we also used a previously developed decoding method based on unsorted hippocampal spikes (Kloosterman et al., 2014; Hu et al., 2018; Ciliberti et al., 2018). Briefly, we estimated a joint probability density function of animal’s position and spike waveform features using kernel density estimation (KDE) methods. To reduce the dimensionality of , we used the spike waveform’s principal components (PCs) of each channel. At each shank, we also ran a kernel compression algorithm to reduce redundancy in the training samples by progressively merging samples based on a compression threshold. In the decoding phase, we computed the likelihood by accumulating the spikes collected from K independent shanks:where n denotes the number of spikes observed at the k-th shank during the interval [t, t+Δt). The generalized rate functions and were defined by the density ratioswhere and denote the joint and marginal densities derived from KDE at the k-th shank, respectively; denotes the mean firing rate at the k-th shank; and denotes the spatial occupancy probability distribution estimated from KDE. Finally, we sought the decoded estimate of among all candidate positions that produced the maximum likelihood. Details of the methods are referred to (Hu et al., 2018).

Spike removal from raw voltage signals

To remove putative spikes from the broadband signal (up to 20 kHz sampling rate), we first detect the spikes by using Kilosort algorithm. Once the spikes were identified in each channel, we removed the spikes, by using a linear interpolation from the start to end points of spike waveform, and further obtained the resulting spike-free signal (see Figure 1F).

Region-specific assessment of decoding

To assess the region-specific contribution to specific decoding strategy, we first split the channels according to their implanted anatomical regions by using current source density map during sharp wave ripple (Figure 2). To control the imbalance of channel counts between regions, we randomly selected the matched number of channels from other regions and repeated the same decoding analysis. We repeated the random selection for 1,000 times and compared their Monte Carlo statistics.

SPW-Rs Replay Analysis

Identification of ripples and awake and sleep replay candidate events

We first used the electromyography (EMG) and LFP for sleep staging. NREM sleep was primarily determined by the low EMG amplitude, high delta/theta power ratio, the presence of slow waves and sleep spindles in LFP activity. REM sleep was determined by the low EMG amplitude and high theta/delta power ratio. SPW-Rs were identified based on a previously described method (Grosmark and Buzsáki, 2016). The integrated power of the filtered LFP signal was calculated in a sliding window for each electrode. To identify the candidate events of memory replay during immobility and NREM sleep (see Figure 3A), we used a combined criterion of hippocampal LFP amplitude at ripple band (140-250 Hz) and total spike count (threshold >mean+3 SD). For visual inspection, we computed the spectrogram to assist identification. We selected subsets of candidate events during immobility (in the maze or sleep box) and post-NREM periods for FPAuhf or clustered spike-based decoding analysis.

Statistical assessment of ripple decoding

During maze run (velocity threshold: 5-15 cm/s depending on the spatial environment), we evaluated the decoding accuracy by the absolute error between the animal’s actual position and the decoded position: (where was derived from a specific decoding strategy using either clustered spikes, FPAuhf or dLFPθ features). We used 10-fold cross-validation to assess the median decoding error, and computed the bootstrapped standard deviation (SD).

During ripple events, we computed the distance correlation (Liu et al., 2018), and then computed 2,000 random shuffles for each event to derive the Z-score statistic or Monte Carlo P-value. For spike-based likelihood decoding, we used cell identity shuffling and receptive field shuffling; for FPAuhf decoding, we used channel order shuffling and linear map shuffling (with respect to the rows of matrix V from OLE). Among different types of shuffling operations, we used the worst Z-score or Monte Carlo p-value as the final shuffle statistics. Monte Carlo P<0.05 was used as the statistical significance criterion.

Real-time online decoding

Real-time FPAuhf decoding based on the OLE method

Our real-time FPAuhf decoding operation consisted of three steps. Step 1: Broadband voltage signals were filtered (>300 Hz) in data acquisition hardware. Step 2: In the memory buffer (20 ms during ripples and 100 ms during run), multichannel LFP signals were processed by a Hilbert transform in parallel (using multi-core CPU), and instantaneous LFP amplitude features were computed in raw 1250 Hz sampling rate and then averaged. The memory was constantly updated in time. Step 3: Position was reconstructed based on a pretrained linear map V from previous session. The pipeline was initially implemented in MATLAB R2019a (MathWorks), and optimized by C/C++ implementation (with Intel C++ compiler). We tested all real-time operations on a desktop (3.6 GHz 4-core Intel i7700 CPU, with Windows10 OS, Visual Studio 2017 and Intel Parallel Studio XE 2018 environment).

Once the LFP features were computed across channels and averaged within each time bin, the remaining computational complexity of OLE operation was O(N∗L) per time bin (where N denotes the number of channels and L denotes the number of spatial bins), which was split into multiple threads in parallel for computational speedup. The computational bottleneck was the Hilbert transform.

Real-time significance assessment of decoded memory replay

To assess the significance of the FPAuhf decoding result, we ran shuffle analyses. To accommodate the real-time computation, we performed and saved two types of pseudorandom shuffles (i.e., channel order shuffle by permutating columns and receptive field shuffle by circular-shifting each column) in advance and directly applied those in online decoding.

We used Hilbert transform to computed the LFP amplitude at ripple band (140-250 Hz) for a preselected channel that shows the largest ripple amplitude, and then compared ripple band amplitude with a predetermined threshold (e.g., mean+SD) to detect the onset of candidate event (Figure 4D). Since the decoding speed is ultrafast, we used a “decode-as-you-go” strategy. Specifically, we calculated >300 Hz high-pass filtered amplitudes features (FPAuhf) for all recording channels and continuously decoded the content for each time bin. But the significance assessment of the candidate event was only executed after the length of event above 3 bins. Once the significance assessment is started, we ran 2,000 (1,000 for each type of shuffle operation) random shuffle decoding analysis for the whole event at each time step and computed the Monte Carlo P-value (based on the distance correlation measurement). After the P-value is calculated, we update the cumulative score as followswhere i denotes the time bin index, and at the time of event onset. Once the online cumulative score was above a predetermined threshold (e.g., used in our current analyses), the event was deemed statistically significant and we reset the cumulative score to be 0.

Unsupervised learning methods

Nonnegative matrix factorization (NMF)

Let Y denote the N-by-T matrix consisting of nonnegative N-channel FPAuhf features. NMF is aimed at finding an approximate low-rank matrix factorization form:where W is a low-rank (m-by-N, m

Deflation of NMF features

Upon the completion of NMF, we reconstructed the specific source of interest in the channel space (i.e. the signal contributed merely to the i-th feature) (Chen et al., 2007):where denotes the matrix pseudoinverse of . By projecting to the original channel’s positions (e.g., channel layout in Dataset 1 in the key resources table), we could yield the channel space map of i-th extracted feature. In addition, to evaluate the relative contribution of every channel to the extracted feature. we need to consider the joint effect of and . For this purpose, we further calculated the weighted estimate of channel space map as follows:where denotes the Hadamard (elementwise) product, denotes the i-th row vector of the matrix , and is the back-projected channel space map from i-th feature. As the result, is the final N-by-T weighted channel space map of i-th extracted feature. To visualize the of all m features, we further calculated the temporal variance of every map and project it to the channel space:where (N-by-1) is the Variance of weighted map of the i-th feature in the channel space, which reflects the contribution of every channel to the i-th extracted feature (Figure S8).

Reconstruction independent component analysis (ICA)

Let Y denote the N-by-T matrix consisting of real or complex-valued N-channel LFP features at a specific frequency band. The ICA assumes that Y is a linear combination of independent source signals: , where is a square mixing matrix, and denotes the N-by-T source signal matrix consisting of mutually statistically independent source signals. ICA is aimed to seek an optimal demixing matrix, , operated on the , such that recovers the independent source signals . For real-valued FPAuhf features, we employed a reconstruction ICA (RICA) algorithm (Le et al., 2011). In contrast to ICA, RICA replaces the orthonormal constraint with a soft reconstruction penalty and allows an overcomplete representation (i.e., is a non-square matrix).

Bayesian hidden Markov model (HMM) for unsupervised learning analysis

To discover latent structures of large-scale hippocampal population codes, we used an HMM for analyzing hippocampal ensemble spike activity during spatial navigation and sleep (Chen et al., 2014, 2016; Linderman et al., 2016; Chen and Wilson, 2017). In a basic HMM, we assumed that the latent state process follows a first-order discrete-state Markov chain {} {1,2,···,N}, and the observations of neural activity at discrete time index t, follow a conditional probability distribution (conditioned on the latent state ). In the case of nonnegative features (either in the form of neuronal spike counts or multi-site FPAuhf features), we assumed a Poisson probability with associated tuning curve functions . In the case of amplitude features, we used a feature resample technique (see the subsection below) to convert the chi-squared distributed amplitude features into Poisson-distributed observations. The joint probability distribution of observed and latent variables is given bywhere denotes an N-by-N state transition matrix, with P representing the transition probability from state i to j; denotes a probability vector for the initial state S1; y denotes the number of spike counts from cell c within the t-th temporal bin (bin size: 100 ms during wake and 20 ms during ripples) and denotes the time series of C-dimensional neural response vector. In the case of spike counts, we assumed the conditional probability distribution had a factorial form: , which defined the products of conditionally independent Poisson distributions. We used a Bayesian inference procedure to identify the unknown parameters and the latent state sequences (Chen et al., 2014). We further used a Bayesian nonparametric version of the HMM, the hierarchical Dirichlet process-HMM (HDP-HMM), which extends the finite-state HMM with a nonparametric HDP prior, and inherits a great flexibility for modeling complex data (Linderman et al., 2016). The associated Markov chain Monte Carlo (MCMC)-based inference algorithm allowed us to infer the number of hidden states N from the inferred posterior distribution.

LFP feature resampling

In the HDP-HMM, the assumed observations were Poisson-distributed spike counts (nonnegative integers). Therefore, we need to adapt the likelihood model of the HMM to accommodate hippocampal field potential observations. In the case of nonnegative features (e.g., derived from NMF operated on nonnegative FPAuhf or dLFPθ power), we generated the same size of random samples from a desired Poisson distribution based on a rank-invariant resampling procedure (Honey et al., 2009; Tu et al., 2020), and replaced the original samples with the ordered new samples (by keeping their rank or order unchanged). We treated as the pseudo spike count observations and then repeated the same Bayesian inference procedure assuming a new Poisson likelihood . For the real-valued Z-scored features, we employed a similar resampling procedure. A schematic illustration of the resampling procedure is shown in Figure S6. The choice of the mean statistic in the new Poisson distribution was ad hoc, yet the final performance was robust with respect to a wide range of the Poisson mean statistic.

Stochastic neighbor embedding (SNE) of multichannel FPAuhf features

The t-distributed SNE algorithm is a robust probabilistic dimensionality reduction method for visualization of high-dimensional data (Van der Maaten and Hinton, 2008). The goal of embedding was to approximate the probability distribution of high-dimensional FPAuhf features with a low-dimensional (n=2 or 3) map, which reflect the similarity between individual samples. The similarity measures can be defined by the user (default: Euclidean distance). During animal’s maze run, we color coded the embedded FPAuhf features in a two-dimensional embedding space according to the animal’s linearized position. We used the ‘tsne’ function (with the default setup) in MATLAB implementation.

Optical flow analysis

Characterization of LFP traveling waves

To characterize the traveling wave patterns of multichannel LFP signals, we used the optical flow (vector field) method. Optical flow is commonly referred to the pattern of apparent motion of objects, and edges in a visual scene. In computer vision, the optical flow methods try to calculate the motion between two image frames which are taken at times t and at every voxel or pixel position. These methods are called differential since they are based on local Taylor series approximations of the frame images; that is, they use partial derivatives with respect to the spatial and temporal coordinates. For a 2D+t dimensional case, a voxel at location () with intensity will be moved by {, , } between the two image frames, and the following brightness constancy constraint can be given:Based on the first-order Taylor series expansion, we derivedwhere and . In the vector form, the above equation is:

We projected the filtered (4-12 Hz or >300 Hz) multichannel LFP signals onto a 2D image according to the recording sites array layout. We then used the Horn-Schunck estimation method (Horn and Schunck, 1981) to estimate the optical flow based on neighboring frames. The optical flow was represented and visualized by a 2D vector field, with arrows indicating the direction, and the size of arrow proportional to the scale. To decoding animal’s position, we further temporally averaged the vector field movie at a frame rate of 10 Hz (i.e., bin size: 100 ms). These high-dimensional flow patterns were used in a simple linear decoding method (such as OLE). In real-time applications, estimation of optical flow can be implemented via graphic processing unit (GPU), such as NVIDIA optical flow SDK (https://developer.nvidia.com/opticalflow-sdk).

Prospective decoding

Prospective decoding and prediction of animal’s future spatial decision-making

In goal-directed navigation tasks, we developed a look-ahead decoding strategy to predict animal’s future decision (L vs. R turns) based on either ensemble spikes, dLFPθ, or FPAuhf features collected in the central arm of T-maze. During the training phase, we correlated the animal’s future position at the L/R arm with neural activity with a time lag (0.5-1.5 s, 100 ms temporal bin size) while animal ran in the central arm (see Figure 7A for schematic illustration). We assessed the “mode” of decoded L vs. R-turn trajectories and determined the prediction outcome based on their sums of weighted scores---the one with a higher cumulative sum was deemed the winner. The score at each time point was weighted with a slow decay from the future to the past. For instance, we used a linear decay weight vector [1, 0.9, 0.8, 0.7, 0.6], by imposing a larger weight on the future predicted outcome (yet the ultimate prediction outcome was insensitive to the exact values in the decay weight vector). To assess the prediction accuracy, we conducted 10-fold cross-validation. In addition, we ran 1,000 Monte Carlo runs (by randomly selecting 90% training trials in each run) and reported the mean±SD prediction accuracy. As control for each strategy, we also performed 1,000 random decoder channel order shuffles and computed the chance-level prediction accuracy.

Furthermore, we trained a linear support vector machine (SVM) classifier (MATLAB Machine Learning toolbox) to predict animal’s future decision based on different features. For the ease of interpretation, we used a linear kernel in the SVM. In addition to the classification accuracy, we also reported the AUROC (area under ROC curve) statistic. An AUROC value of 0.5 indicates a chance level, whereas a high AUROC value (close to 1) implies excellent performance.

For all the prediction analyses, we only used the trials that animal start from the beginning of center arm and run directly to the next left/right arm without hesitation or stopped before the choice point. By excluding these invalid trials, we reduce the possibility that animal will change its decision throughout center arm run.

Quantification and statistical analysis

Position decoding errors were quantified by the probability distribution function as well as the median statistic. Additional quantification of decoding error was computed with respect to the position or behavioral sampling (occupancy). To assess the statistical significance, Monte Carlo shuffles were conducted to calculate Monte Carlo P-values.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
Datasets 1, 2	https://crcns.org/data-sets/hc/hc-11	N/A
Datasets 3, 4	http://buzsakilab.com/wp/datasets/	N/A
Dataset 5	https://crcns.org/data-sets/hc/hc-3	N/A
Dataset 6	http://buzsakilab.com/wp/datasets/	N/A
Experimental models: organisms/strains
Long-Evans rats	Charles River Labs	RRID: RGD_2308852
Mice C57BL/6J	The Jackson Laboratory	Stock No.: 000664
Software and algorithms
MATLAB	MathWorks	N/A
PYTHON	Open source	N/A
VB-ARD linear regression	Drugowitsch, 2013	https://github.com/DrugowitschLab/VBLinLogit
Custom code for various decoding methods	This paper	https://doi.org/10.5281/zenodo.5527940