Abstract
Neural activity at the large-scale population level has been suggested to be consistent with a sequence of brief, quasistable spatial patterns. These “microstates” and their temporal dynamics have been linked to myriad cognitive functions and brain diseases. Most of this research has been performed using EEG, leaving many questions, such as the existence, dynamics, and behavioral relevance of microstates at the level of local field potentials (LFPs), unaddressed. Here, we adapted the standard EEG microstate analysis to triple-area LFP recordings from 192 electrodes in rats to investigate the mesoscopic dynamics of neural microstates within and across brain regions during novelty exploration. We performed simultaneous recordings from the prefrontal cortex, striatum, and ventral tegmental area in male rats during awake behavior (object novelty and exploration). We found that the LFP data can be accounted for by multiple, recurring microstates that were stable for ∼60–100 ms. The simultaneous microstate activity across brain regions revealed rhythmic patterns of coactivations, which we interpret as a novel indicator of inter-regional, mesoscale synchronization. Furthermore, these rhythmic coactivation patterns across microstates were modulated by behavioral states such as movement and exploration of a novel object. These results support the existence of a functional mesoscopic organization across multiple brain areas and present a possible link of the origin of macroscopic EEG microstates to zero-lag neuronal synchronization within and between brain areas, which is of particular interest to the human research community.
SIGNIFICANCE STATEMENT The coordination of neural activity across the entire brain has remained elusive. Here we combine large-scale neural recordings at fine spatial resolution with the analysis of microstates (i.e., short-lived, recurring spatial patterns of neural activity). We demonstrate that the local activity in different brain areas can be accounted for by only a few microstates per region. These microstates exhibited temporal dynamics that were correlated across regions in rhythmic patterns. We demonstrate that these microstates are linked to behavior and exhibit different properties in the frequency domain during different behavioral states. In summary, LFP microstates provide an insightful approach to studying both mesoscopic and large-scale brain activation within and across regions.
Introduction
The coordination of neural activity across the entire brain has remained incompletely understood. At the same time, recent research has demonstrated the far-reaching influence of behavior and other state changes on the activity in distant and substantial parts of the brain (McGinley et al., 2015; Musall et al., 2019). Previous work has demonstrated that neural activity is structured into short-lived epochs of coordinated activity on the local scale (e.g., of small ensembles of cells; Luczak et al., 2009, 2015) and on the large scale (e.g., measured noninvasively; Lehmann et al., 1987; Pascual-Marqui et al., 1995).
On the local level, transient spatiotemporal patterns have been observed in visual cortex (Arieli et al., 1995; Chiu and Weliky, 2001; Kenet et al., 2003; Fiser et al., 2004; Omer et al., 2019), auditory cortex (Luczak et al., 2009; Sakata and Harris, 2009), and other cortical areas in the awake, anesthetized, and sleep state in various species using different recording techniques (Mao et al., 2001; Petersen et al., 2003; Battaglia et al., 2004; Massimini et al., 2004).
On the large-scale level, electroencephalography (EEG) has been a valuable tool to study whole-brain activity in humans. Here, transient spatiotemporal patterns have been identified and termed “EEG microstates.” EEG microstates are defined by clustering topographical maps around peaks in global field power (GFP; Lehmann et al., 1987) using modified k-means clustering (Murray et al., 2008). A large number of studies have shown that four clusters account for a large fraction of variance and are observed consistently across myriad studies spanning nearly 4 decades (Michel and Koenig, 2018). Microstates are used to study human cognition, including visual processing (Britz and Michel, 2011), perceptual awareness (Britz et al., 2014), neuropsychiatric disorders including schizophrenia (Dierks et al., 1997; Lehmann et al., 2005; Kindler et al., 2011; Tomescu et al., 2015; da Cruz et al., 2020), and resting-state networks (Britz et al., 2010; Musso et al., 2010; Yuan et al., 2012; Khanna et al., 2015; Bréchet et al., 2019; Zanesco et al., 2020).
The relation between large-scale EEG microstates and small-scale neural population spike correlations is poorly understood, in part because different analysis methods are applied at different scales (Alishbayli et al., 2019). Here we combined multichannel recordings (up to 192 electrodes across three brain areas) with the analysis techniques developed in the context of EEG microstates to investigate the dynamics of neural activity on the meso-scale [i.e., local field potentials (LFPs)]. We recorded LFPs simultaneously from the prefrontal cortex (PFC), striatum (STR), and ventral tegmental area (VTA) in rats, while they moved freely in an arena in which novel objects were introduced.
We observed reliable LFP microstates in each region that shared characteristics previously identified in EEG microstates, including robustness, temporal dynamics, temporal continuity, and relation to global field power (Khanna et al., 2015; Milz et al., 2016; Mishra et al., 2020). Using concepts from EEG microstate analysis, we verified the validity of these meso-scale, LFP microstates. For example, the LFP microstates explain a large proportion of the variance of the signal, in contrast to their temporally shuffled surrogates. Furthermore, we linked different microstate properties like occurrence rate, temporal coverage, and state duration to behavioral states. A striking finding was that LFP microstates appear to be coordinated across regions: microstate time series from distant brain regions exhibited distinct cross-correlation patterns that depended on the behavioral state, suggesting functional interactions of these brain regions, manifested by the microstate activation patterns. Functional coupling of these regions is a well established result (Gao et al., 2007; Zhang et al., 2016), which can be analyzed from a larger perspective using LFP microstates. Thus, we propose that the EEG microstate method is an insightful approach to identify, define, and characterize more local, transient states of neural activation. As proposed by Luczak et al. (2015), these microstates might form a vocabulary of neural activation that could help structure our understanding of whole-brain dynamics.
Materials and Methods
All experimental procedures were performed in accordance with the European Union directive on animal experimentation (2010/63/EU), and the Dutch nationally approved ethics project 2015–0129. All recordings were performed in the laboratory of M.X.C. We included four male Long–Evans TH:Cre rats (age, ∼3 months; weight, 350–450 g at the time of recordings). Rats were housed singly in Makrolon type III cages (40 × 30 × 35 cm; UNO) in a room with controlled temperature (21 ± 2°C) and humidity (60 ± 15%). Animals were kept on 12 h light/dark cycles, and the experiments were performed under the light phase.
Electrode implants
Each silicon probe was loaded on a plastic 3D-printed drive purchased from 3Dneuro (https://www.3dneuro.com). The 3Drives were lightweight (0.43 g), with a small footprint and base dimensions of 5 × 6 mm and height of 16.7 mm, allowing the travel distance of 9 mm and resolution of 62.5 µm per one-quarter turn. Rats were anesthetized with 5% isoflurane (2–3% maintenance) and placed in a stereotaxic frame. They were given a dose of Rimadyl/Carprofen (0.5 mg/kg) and Xylocaine (maximum, 0.4 ml) subcutaneously at the incision site. Small craniotomies on the skull were made above the VTA, STR, and PFC. The animals were implanted with custom-designed high-density silicon probes (NeuroNexus) with specific electrode layout of densely packed recording sites designed for each target region (Fig. 1A, red dots inside the coronal slices).
In PFC and STR, each probe contained four laminar shanks, having 16 recording electrodes per shank in forms of tetrodes for spike recordings and single sites for LFP recordings, thus creating 64 recording sites per probe. Since the focus of the current investigation is on LFP signals, spike recordings are not analyzed here. Sixty-four electrodes covered an area of 1 × 2 mm with typical spacing of 225 µm in each shank and 330 µm between shanks in PFC.
STR electrodes also covered an area of 1 × 2 mm with the same shank distance (330 µm). However, two shanks contained only tetrodes, and two shanks had only single sites with typical spacing of 130 µm between single sites and 660 µm between tetrodes.
VTA contained eight shanks of eight electrodes each and covered an area of 1.5 × 0.14 mm. All electrodes were referenced to a reference screw implanted on the skull over the cerebellum.
Implants were surrounded by and grounded to a copper mesh Faraday cage built around the implants and was fixed to the skull surface using dental cement (Vandecasteele et al., 2012). Rats were monitored daily for 1 week following surgery and recordings began at least 1 week after surgery. After recovery (minimum, 7 d), the electrodes were advanced into the brain toward a final target location of 4 mm below the skull for PFC, 5 mm for striatum, and 8 mm for the VTA.
Experiment design and apparatus
Animals were handled and trained for 1 week to freely move and get habitualized in a black plastic square box (length × base × height, 60 × 40 × 40 cm) open field covered with bedding in the bottom of the enclosure. The square enclosure was placed on an electrically grounded Faraday cage inside the recording room. In each session, rats were allowed to move freely in the box. A novel object (e.g., a cup or a toy) was presented in the middle of the box in the second session, and the same object was repeated in the fourth session (see Fig. 4A). Rats were allowed to explore the novel object freely. A camera was placed above the box to track movement. Camera resolution was 320 × 240 (width × height) and captured at a frame rate of 30 Hz. The alignment between camera and OpenEphys acquisition system was implemented using event markers. All sessions within a trial were recorded without any intersession delay. A maximum of one trial per animal was recorded on a single day. There were 27 recording trials in total.
Recording and cleaning of electrophysiology data
Electrophysiological data were sampled at 30 kHz via OpenEphys hardware and software solutions (Siegle et al., 2017). Offline, data were downsampled to 1 kHz and stored in the EEGLAB (Delorme and Makeig, 2004) format in MATLAB (MathWorks). Excessively noisy data segments were removed manually, and excessively noisy channels were removed from the data based on visual inspection of raw signals and variance thresholds (6.6 ± 5.2 channels per dataset per region). Per recording, removed segments comprised 6.2 ± 5.0% (median ± median absolute deviation) of all movement states and 2.0 ± 1.5% of all stationary states. Similarly, removed segments comprised 2.8 ± 2.8% of all object interaction states and 2.7 ± 2% of nonobject interaction states. Thus, removed noisy segments are significantly smaller than the rest of the datasets (p < 1e-14, Wilcoxon signed-rank test) in each behavioral state, ensuring that sufficient data are available for subsequent behavioral state-specific analyses. More importantly, several microstate features that we analyzed were normalized to the amount of data (e.g., the rate of occurrence), and the surrogate-based permutation tests shuffle the retained data and are thus robust for the amount of data. Thereafter, average referencing was performed separately for the set of electrodes located in each region. The datasets were then bandpass filtered in the range of 1–30 Hz, using a Hamming window sinc FIR (finite impulse response) filter, to be consistent with the common practice in EEG microstate analysis (Milz et al., 2017).
Microstate analysis
The microstate analysis was performed using the Microstate EEGLAB Toolbox (Poulsen et al., 2018). Before running the microstate analysis, data from all sessions were concatenated to create one continuous epoch of ∼20 min. Extracting microstates on the aggregated data improves the quality of the decomposition and facilitates a direct comparison of microstate dynamics across sessions.
Detailed descriptions of microstate extraction can be found in previous publications (Pascual-Marqui et al., 1995; Murray et al., 2008; Milz, 2016; Michel and Koenig, 2018; Poulsen et al., 2018; Mishra et al., 2020). Here we provide a brief overview; the following procedure was performed separately for each animal and each recording session. The first step is to compute GFP time series, which is the SD across channels. Microstates are extracted at local maxima of GFP, with the constraint that two GFP peaks must be separated by at least 10 ms. The second step is to cluster the spatial maps at GFP peaks according to a modified k-means method that has been adapted for EEG microstate analysis. We selected four clusters for subsequent analysis. The number of clusters was chosen based on the global explained variance (GEV) criterion (Eq. 2; Fig. 1C1–C3). We noticed that the incremental change in GEV saturated after three to five clusters. Moreover, cluster centers with similar spatial patterns started emerging after four clusters in most cases. Although the “correct” number of clusters is difficult to determine exactly, a consistent number of clusters across all recordings facilitates analyses and comparisons. Thus, we decided to extract four clusters for all datasets. Because k-means involves random initializations, we repeated the clustering algorithm 10 times using 10 different random starting values; the iteration with the highest total GEV was retained for the final analyses. The third step in microstate analysis is to “backfit” the microstates onto the data, which provides a time series containing the dominant microstate map at each time point; thus, successive time points with the same map are termed “microstates.” This is achieved by correlating the spatial map at each time point with the four microstate maps, and labeling each time point according to the map with the largest absolute correlation (absolute correlation is used to ensure polarity invariance of the backfitting). We set the minimum duration of a microstate to be 25 ms, and any microstates shorter than this are absorbed into the neighboring microstate.
Microstate properties estimation
We quantified the temporal pattern of occurrence of microstates using the following five properties.
Mean spatial correlation.
Defined as the mean Pearson correlation between LFP signals and their corresponding microstate map. We used the absolute value of the correlation for polarity-independent analyses, as follows:
GEV.
Global explained variance is a measure of similarity between the data and the microstates. GEV for the nth sample of the data is defined as follows:
Occurrence rate.
The mean number of appearances of a microstate per second.
Temporal coverage.
The fraction of total time covered by all occurrences of a given microstate.
State duration.
Defined as the median duration of each microstate occurrence.
Although GEV and spatial correlation are correlated, they are separate measures, and we therefore report both.
Animal movement tracking and behavioral state determination
Animal movement tracking was performed using DeepLabCut (Mathis et al., 2018). We used this deep neural network-based method to track five body markers (i.e., snout, left ear, right ear, head center, and tail start; see Fig. 4A1,A2). We also tracked the location of the object center, for sessions in which an object was present. This was done because the objects could be moved around slightly as the animals explored it. A moving median filter with a window length of 13 frames was used to postprocess the movement-tracking results. This improved accuracy by removing isolated video frames with incorrect location of body markers and objects.
To determine the video frames where the animal was not moving, a threshold on the speed of the animal was manually set based on visual inspection (see Fig. 4D2).
We defined a radius around the object center, based on the size and shape of the object and determined separately for each object. For the extraction of frames where the rats were interacting with the object, the distance between the snout and the object center was computed. “Object interaction” was defined as frames in which this distance was less than the interaction radius. Two binary vectors, associated with two types of behavioral states–object interaction and movement, were upsampled to 1 kHz to match with the LFP data.
Microstate activation cross-correlation and their power spectra
The ephys data were projected onto the extracted microstates to generate microstate activation time series. Cross-correlations were performed on all pairs (within and across regions) of these microstate activation time series, up to a maximum lag of ±2000 ms (using xcorr in MATLAB). The power spectra of the cross-correlograms [cross-spectral density (CSD)] were estimated by taking the square of the discrete Fourier transform (fft in MATLAB) divided by the length of cross-correlograms. To facilitate comparison and subsequent analyses, power spectra were z-normalized by taking the SD over frequencies.
Power-spectrum clustering
The cross-correlation power spectrum profiles were clustered using the k-means method. Squared Euclidean distance criteria was used as the distance metric for the clustering. The optimal number of clusters, determined using silhouette measure, was two. Clustering was performed separately for each dataset.
Surrogate data for permutation testing
To test the statistical significance of GEV by four microstates, and the effect on microstate properties by the change in behavioral states, we generated surrogate datasets in which the local temporal structure per channel was preserved, while the temporal alignment between channels was randomized. This was achieved in the LFP data by choosing a random point in time for each channel and swapping the time series around these time points. This surrogate technique was chosen because it not only destroys the spatial alignment of the activity, but preserves the local temporal structure (see Fig. 3A–B). A complete randomization of the time points would remove all spatial and temporal structures, and thus result in a “straw-man” null-hypothesis against which any empirical microstate would be highly significant. One hundred surrogates were produced, and GEV values were calculated for each surrogate to compare with the GEV of the empirical data. We also generated surrogate movement data using the same method. Again, 100 surrogates were produced. Microstate properties for each microstate were calculated for the surrogate movement data and compared against their empirical value to compute the significance (see Fig. 5C1–C3).
To determine a significance threshold for the cross-correlation analyses, we used 50 randomly shuffled surrogates per dataset (total of 50 × 27 = 1350 shufflings), and computed microstates and cross-correlation spectra from each shuffled dataset, as described above. From each of 27 datasets, we computed the 99th percentile of the shuffled CSD distribution. Finally, we took the median of the 99th percentile spectra across the 27 datasets as our final statistical threshold seen in Figure 6, E and F.
Statistical analysis involving behavior
Statistical analysis on microstate properties and their relationship with the behavioral states was performed without distributional assumptions, by using permutation tests. Statistical analysis on behavioral state durations and power spectra peak frequencies were performed using a two-sided Wilcoxon rank-sum test.
Results
We investigated the existence and dynamics of coordinated neural activity over three brain regions (PFC, STR, and VTA) in chronic recordings in behaving rats (N = 4; Fig. 1A). Recordings were collected simultaneously using up to 192 electrodes (64 per region) in conjunction with video recordings. We used microstate analysis, an established technique for detecting short-lived states of joint activation from the EEG literature (Lehmann, 1989; Michel and Koenig, 2018; Mishra et al., 2020) and deep learning-based video tracking (Mathis et al., 2018) to relate the occurrence of the detected microstates with ongoing behavioral states. We start by introducing the analysis technique to isolate microstates in multichannel LFP data.
LFP signals recorded from multiple sensors can be represented by LFP microstates. A, Neural recordings were concurrently collected from three brain regions (PFC (navy), STR (maroon), VTA (orange) in awake behaving rats. Large-scale electrode arrays (red dots) covered a coronal plane in each brain area (slices show approximate anteroposterior locations, taken from the Sprague Dawley rat atlas by Papp et al., 2014). The probe contained 16 recording electrodes per shank in PFC and STR. B, LFP data excerpt from all channels recorded in each brain region (3 × 64, top) and GFP (bottom). Microstate maps were extracted based on the local peaks of the GFP, using a modified k-means clustering algorithm (Michel and Koenig, 2018; Poulsen et al., 2018; Mishra et al., 2020; for details, see Materials and Methods). C1–C3, Cumulative global explained variance saturated after three to five microstates in all three regions. Gray lines indicate results from individual datasets, thick lines indicate the mean. D, The spatial patterns of LFP signals (top left) can be captured by a small number of microstate maps (top right; e.g., data and corresponding microstates at two GFP peaks, 331 and 751 ms). We focused on the four dominant microstates in each region (bottom, different colors). The current microstate is determined by correlation between microstates and the current LFP. E, Microstate maps extracted from each region while the animal was freely moving in an open field environment (A, top left). The left-most maps also show the electrode locations. Map values were spatially spline interpolated between electrode locations. Maps in D and E have the same color scale as that shown in D. D and E are from rat 1, recording day 1.
Large-scale LFP signals are accounted for by a small set of microstates
Multiarray, multiarea LFP signals were acquired simultaneously (Fig. 1B) and subjected to microstate analysis separately per region. The microelectrode array geometry enabled the acquisition of neural signals from a wide area (0.15−2 mm2). In each area, the spatial activity could be accounted for by a small number of microstates (i.e., joint activity patterns), whose residence time ranged between a few tens and a few hundreds of milliseconds (Fig. 1D, sample from PFC). Microstates are recurring spatial patterns of neural recordings. Since all similar maps can be grouped together, the representative structures from each group can be considered as the representation of the entire neural signal. Thus, the segmentation of the LFP signal into microstates is based on spatial correlation of neural activity (for details, see Materials and Methods). Qualitatively, this becomes apparent when the LFP data and corresponding microstate are visualized together (Fig. 1D, top; quantified in Table 1). The first four microstates accounted for 68% of the variance in the example shown (Table 1, details; see Fig. 8, region-specific details). The exact number of significant microstates from a clustering analysis can be difficult to ascertain (Murray et al., 2008; Michel and Koenig, 2018); we chose to use four microstates per region for each recording day, based on the total amount of explained variance (Fig. 1C1–C3), to facilitate analyses across days. Although other choices may yield slightly different results, our experience is that the patterns of findings are robust to the exact number of microstates.
Summary of microstate properties
Microstate topographies exhibited a mixture of distributed spatial features (Fig. 1E, left, Microstate1 for STR) and restricted local activity (Fig. 1E, middle left, Microstate2 for STR). Importantly, electrodes in the array were individually inspected for spurious recordings (e.g., caused by unreasonably high variance or large differences in impedance; for details, see Materials and Methods). Electrodes with excessive variance or artifacts were excluded from subsequent analysis, ensuring that the microstates reflected spatiotemporally coherent activation patterns based on well sampled neural activity.
Microstate maps are robust over multiple days
We repeated the experiment in the same animals over multiple days to determine the stability of microstate maps over time. The microstate analysis was performed separately on the recording of each day (between 5 and 9 d/animal). We found that spatial correlations of microstate maps across days were high (Fig. 2), indicating that many of the maps recurred on different days. Other microstates were less stable and were identified on some but not all days (Fig. 2A1, blue/green, Microstate3/4 in PFC). It is possible that these “missing” microstates were sorted into smaller clusters that were excluded from the primary analysis because of lower signal-to-noise ratios or infrequent occurrences.
Microstate maps are reproducible across different days. Simultaneous LFP signals were recorded on six different days (each day with a new object) in PFC, STR, and VTA. Microstate maps from each region on the first day were chosen as reference maps. Maps are rotated counterclockwise by 90° for PFC and STR and rescaled for visualization purposes. A1, In the PFC, the squared correlations among the four maps on each day and average microstate maps were computed, and are shown with red, black, blue, and green lines for microstates 1–4, respectively. Lower correlation values on certain days can be loosely interpreted as the absence of the Microstate on that day (in reference to the average Microstate). A2, Microstate (MS) maps are reordered based on their best match (A1). Most microstate maps retained the spatial features on each recording day. (A1 and A2 are from animal 1 in A3). A3, The fraction of maps that were repeated over time for each of four animals. B, C, These rows show the same results as described above, but for data from the striatum and VTA.
To assess the presence of microstate maps on experiment repeats, we averaged all highly similar maps (per region and per animal) together as their first principal component. We refer to these maps as subject average microstate maps. Squared spatial correlations between the subject average microstate maps and microstate maps on each experiment day were computed as an indicator of microstate recurrence (Fig. 2A1,B1,C1, example from one animal). We used a threshold of 0.5 (squared spatial correlation) to determine that a microstate repeated, although the exploration of a range of threshold values (0.4–0.7) yielded qualitatively similar results. In all three regions, the same set of microstates, within the same animal, showed a stronger repeatability over separate recording days (Fig. 2A3,B3,C3).
Spatiotemporal LFP structure underlies the explanatory power of microstate maps
The previous analyses demonstrate that microstates are rather stable over multiple days. We next investigated how much variance in the total signal is explained by the microstates. We quantified this using the GEV, which reflects the correlation between the signal at each time point and the microstate maps. Using four microstates per region accounted for 71 ± 6%,71 ± 9%, and 64 ± 10% (mean ± SD) of the global variance in PFC, STR, and VTA, respectively. We tested whether these GEV values were statistically significant by generating a null-hypothesis distribution, generated by shuffling the data in a way that preserved the local temporal structure but distorted global temporal structure (see Materials and Methods; Fig. 3A–B, illustration). Microstate analysis was rerun 100 times with different random surrogate datasets, and subsequently the GEV was recomputed.
Microstate maps represent the spatiotemporal structure of LFP activity. A1, Surrogate LFP datasets with conserved local structure were created by pivoting the data around a random time point (independently selected per trial; blue, before pivot point; red, after pivot point). A2, The data were swapped around these pivot points to generate a surrogate dataset. A1 and A2 illustrate the principle of the shuffling method; the actual recordings were ∼20 min long. The shuffling strategy disrupts the spatial structure while imposing minimal disruption to the time series within each channel (vertical line at 300 ms as an illustration). B, Microstate analysis was performed on these surrogate datasets, and GEV was computed. The corresponding GEV distribution for one of the recordings (C, red circle) showed a significant decrease in GEV in the surrogate data, indicating that the original temporal alignment between channels is underlying the microstate maps explanatory value. C, GEVs of the empirical data were always higher than the corresponding median GEV from surrogate datasets in all three brain regions, showing statistical significance of the microstate results.
We found that this removal of spatiotemporal LFP patterns in the surrogate data resulted in a significant drop in GEV (p < 1e-35; Fig. 3C). The GEV of the shuffled data were only 41 ± 9%, 37 ± 5%, and 31 ± 4% compared with ∼70% for the original data (see above), corresponding to a significant reduction by 43%, 47%, and 51% (p < 1e-35, permutation test) in PFC, STR, and VTA, respectively. GEV values differed considerably across regions (Fig. 3C), but in all cases using the spatial surrogates dropped the GEV significantly. In conclusion, microstates capture repeatedly occurring spatiotemporal activation patterns, a small number of which are sufficient to explain up to ∼80% of the total LFP variance.
Region-wise correlations between shuffled data GEVs and empirical data GEVs were not significant (Fig. 3C; PFC: r = 0.33, p = 0.09; STR: r = 0.23, p = 0.24; VTA: r = 0.058, p = 0.78).
Linking microstates to behavioral states
During the recordings, animals interacted with a novel or a familiar object in the arena. We subdivided the recording session along the following two behavioral dimensions: (1) “movement related” (i.e., stationarity or in motion); and (2) “object exploration related” (i.e., interacting or not interacting with the object). Tracking of the animal and the object was performed using DeepLabCut (version 2.1; Mathis et al., 2018), and each time point was assigned a binary value for each of the two dimensions based on empirically determined thresholds (Fig. 4; for details, see Materials and Methods). Each recording day comprised four sessions (Fig. 4B); movement tracking was performed separately for each session and then later pooled.
Automated movement tracking to segregate behavioral states. A1, An example video frame during the experiment. A2, Animal movement is accurately tracked using DeepLabCut, shown here using the same frame as in A1. Animal movement was tracked using five body markers (blue, black, red, green, and yellow dots represent, respectively, snout, tail start, left ear, right ear, and head center). B, Each experimental recording consisted of four sessions (∼5 min each session). Sessions 1 and 3 were open field recordings, while session 2 included a novel object (different each day) in the middle of the box. In session 4, the same object was present. Rats spend less time interacting with familiar objects than novel objects (each dot represents snout location in each frame; red dots show time points labeled as interaction). A small number of erroneous object locations were filtered out before analysis. C, The distribution of temporal coverage of each behavioral state (rest, movement, object interaction) across sessions was dominated by periods of rest. Object interaction can take place during rest as well as during movement, represented as an overlap in behavioral states. D1, Median temporal coverage from all recordings captures the trends over all sessions, recordings, and animals (N = 4). D2, Median distribution of animal speed movement from all the recordings. The dashed line shows the motion speed threshold separating resting from movement.
Animals spent most of their time in the stationary state (median ± SD, 88 ± 5%) followed by movement state (12 ± 5%). Next, we classified these two behavioral states further into object interaction and no object interaction and found that the fraction of time spent on object interaction mostly took place during stationary state (5 ± 4%) compared with interaction during moving state (0.9 ± 0.7%; Fig. 4C). Per recording, we removed 6.2 ± 5.0% (median ± median absolute deviation) of all movement states and 2.0 ± 1.5% of all stationary states. Similarly, removed segments consisted of 2.8 ± 2.8% of all object interaction states and 2.7 ± 2% of non-object interaction states. The rejected segment lengths are significantly lower than the remaining data in each behavioral state (p < 1e-14, Wilcoxon signed-rank test). Thus, overall only a small fraction of data from all the behavioral states is removed during data cleaning.
We next examined the evolution of behavioral states as time progressed. First, we compared the sessions in which the object was not present (sessions 1 and 3). We found that the animals spent more time in a stationary state in the later session (session 3, 93 ± 6%) compared with the earlier session (session 1, 85 ± 8%; p = 0.00,017, two-sided Wilcoxon signed-rank test). Next, we compared the time spent in the stationary state when the object was introduced (sessions 2 and 4). We again found that animals spent more time in the stationary state in the later session (90 ± 5%) compared with the earlier session (83 ± 6%; p = 0.00,017, two-sided Wilcoxon signed-rank test). A novel object was first presented in session 2, and the same object was again presented in session 4. The fraction of time spent on object exploration in session 4 (9 ± 12%) was significantly less (p = 0.0021, two-sided Wilcoxon signed-rank test) than in session 2 (15 ± 12%). There were no trends along experiment repeats, indicating that the animals remained interested in the novel objects each day.
Microstate dynamics were modulated by behavioral state
We quantified three microstate properties during each behavioral state—occurrence rate, temporal coverage, and state duration (for definitions, see Materials and Methods)—and reported their summaries in Table 1 (see Fig. 8). We separated these features for each microstate and brain region, and for the behavioral states “stationary” and “movement.” All three microstate properties in all three brain regions were modulated by behavioral state (Fig. 5A, example from one recording in PFC). Furthermore, these modulations were consistent across the different recording days (Fig. 5B1–B3, one animal in PFC). Microstates followed the same direction of state-related changes in most of the repeats (PFC, 76 ± 13%; STR, 82 ± 14%; VTA 82 ± 15%), indicating a systematic relationship between microstate properties and behavioral states.
Microstate parameters are affected by behavioral state. A, Example results from the recording in PFC (rat 1, recording day 3). A1, Temporal coverage fraction for all four microstates differed in the two behavioral states (statistics are in C1). A2, Occurrence rate during the two behavioral states (solid bars on the left refer to stationary periods, and lighter bars on the right refer to movement periods; the four numbers refer to the different microstates). A3, Similarly, microstate durations were modulated by behavioral state. B1–B3, The modulation of microstate properties by behavioral state (expressed as stationary minus movement) was consistently in the same direction (PFC, 76 ± 13%) across 6 different recording days (each dot corresponds to a day; only microstates with r2 > 0.5 repeatability are shown here). C1–C3, Similar to B but for data pooled across 27 novel objects in four rats. Significance was assessed by comparison against surrogate datasets (thin horizontal lines indicate a significance threshold of z = ±1.96). There was a high fraction of significant results for all three microstate properties in all three brain regions.
We assessed the statistical significance of the behavioral state dependency by generating a null-hypothesis distribution based on surrogate behavioral states, obtained by random temporal cuts of the time series of behavioral labels (Fig. 3A; also see Materials and Methods). Temporal coverage changed significantly between behavioral states, as all three regions showed a high fraction of significant changes (|z| > 1.96) exceeding 0.05 (fraction of significant tests: PFC, 0.76; STR, 0.79; VTA, 0.66). The occurrence rate showed a similar trend. The changes in the occurrence rate were also significant in all three regions (PFC, 0.35; STR, 0.44; VTA, 0.52). State duration was less affected by behavior (PFC, 0.22; STR, 0.23; VTA, 0.08; Fig. 5C). We also repeated these analyses for object exploration-related behavioral states and observed qualitatively similar results.
In summary, temporal dynamics of microstates were significantly affected by changes in behavioral state. The most significant changes were observed for temporal coverage.
Relationship of microstate dynamics across brain regions
We speculated that if microstates reflect a mesoscale organizational principle, then functional connectivity across brain regions might manifest as correlated microstate temporal dynamics. To investigate this possibility, we computed cross-correlations across all pairs of microstate activation time series, defined as the projection of sensor-level signals at each time point onto the microstate basis vectors, in the three regions (144 cross-correlation pairs in total, of which 12 are autocorrelations and 66 are inter-regional, omitting pairs that only differ by sequence; Fig. 6A). All cross-correlation and autocorrelation curves followed a damped oscillatory nature with an average peak of 9 ± 6 ms. The peak amplitude of cross-correlations varied substantially across microstate combinations and across animals (mean ± SD; absolute correlation value, 0.17 ± 0.07).
Microstate activation time series are rhythmically correlated within and between brain areas. A–E2, Illustrative results from one recording session (rat 1, recording day 2). A, Microstate time series cross-correlations show damped oscillatory dynamics across microstates both within regions and across regions. B, Example cross-correlation between activation time series of microstate2 from STR and microstate2 from VTA (this is a zoomed-in version of the boxed plot in A). The typical cross-correlation structure is absent in shuffled data (red line). C, The z-normalized power spectrum of the cross-correlation (cross-spectral density) shows a peak at ∼8 Hz. The shuffled data do not exhibit a dominant peak (red line). D, The z-normalized peak power from all combinations of microstate activation cross-correlations illustrates varying degrees of correlation within and between areas. E1, E2, The z-normalized power spectra from all combinations of microstates within and between areas were clustered into two groups, one with a lower-frequency peak (∼2-4 Hz; E1, black lines,) and one with a higher-frequency peak (∼8 Hz; E2, blue lines). Thick lines indicate cluster centers, and thin lines are individual power spectra. The red lines indicate CSD power for statistical significance calculated from permutation testing (p < 0.01). F1–F2, Results from all recordings in all animals show a consistent grouping of power spectra based on the frequency of the peak in power. The red lines indicate CSD power for statistical significance calculated from permutation testing (p < 0.01).
Remarkably, the cross-correlations exhibited robust oscillations, suggesting rhythmic patterns of temporal microstatic activation coordination across microstates and regions. The z-normalized power spectral profiles (Fig. 6C) contained a distinguishable set of peaks at different frequencies. Typically, one or multiple peaks were <10 Hz, along with smaller peaks in higher frequencies (10–30 Hz). These higher-frequency peaks appeared to be harmonics of the lower-frequency peaks, and we therefore focused on the lower frequencies. As expected from the cross-correlation profiles, peak power varied for different microstate combinations, suggesting the varying strength of cross-region interaction in different microstate combinations (Fig. 6D). To test the statistical reliability of these cross-correlation profiles, we repeated the analysis on shuffled surrogate data (see Materials and Methods). The damping oscillatory nature of cross-correlation structure was absent in the shuffled data (Fig. 6B, red line), indicating that the cross-correlation magnitude reflects microstate interactive temporal dynamics and not random fluctuations.
The power spectral profiles of the cross-correlations were grouped into two clusters, based on k-means and the silhouette method (Rousseeuw, 1987). The first cluster had a dominant peak at 2–4 Hz, while the second was dominated by peaks at ∼8–10 Hz (Fig. 6E1,E2). These spectral profiles were consistent both within animals (across experiment repeats) and across animals. All individual spectral profiles showed a consistent shape (Fig. 6F1,F2).
The occurrence of large-scale EEG patterns could be enabled by underlying co-occurrences of microstates with near-zero lag. While volume conduction prevents such an analysis in EEG recordings, our data allowed investigations into near zero-lag connectivity. For this purpose, we identified the peak in the cross-correlation with a lag closest to zero (“near-zero peak lag”). We found that the peak lag values were largely repeatable over experiment days within each animal. Computing confidence intervals (95%; Jackknife method) for each inter-region microstate pair on each experiment day led to the finding that >92% of inter-region microstate pairs across all experimental days in each animal included 0 ms in the 95% confidence interval around the peak, suggesting a substantial fraction of co-occurrence of microstates across regions.
In summary, microstate activations were temporally synchronized across brain regions both in slow (2–4 Hz) and theta/alpha (8–10 Hz) ranges. Moreover, most of these synchronizations indicated correlation peaks at zero or near-zero lag. This suggests that mesoscopic-scale networks interact in a synchronized rhythmic fashion across the brain areas.
Inter-regional interactions and behavioral states
We investigated the possible behavioral relevance of the microstate activation correlations by separating them according to the following two behavioral states: object interaction and long stationary periods. We defined a “long stationary period” as the trial segments where the animals remained stationary for >90th percentile of all stationary segment lengths (2.3 ± 0.9 s; calculated separately for each trial). While it is easier to distinguish moments of object interaction and no interaction, shorter segments of stationarity can be indistinguishable with movement preparation. Hence, to compare against the dynamics during the “true” stationarity, we chose to focus on longer stationary periods. Next, we defined dominant microstates during these two behavioral states in terms of their temporal coverage (Fig. 7A1–A3). We focused on temporal coverage because it showed the most robust dependence on object exploration in all three regions. All three microstate properties exhibited dependence on the behavioral state (Fig. 5), and their results were qualitatively similar. However, temporal coverage is proportional to both state duration and occurrence rate and, hence, was considered as the most suitable property for the subsequent analyses. The cross-correlation power spectra during long resting periods dominated the 8–10 Hz peak, whereas the power spectra during object interaction dominated the 2–4 Hz peak (Fig. 7B1–C2).
Power spectra of microstate activation cross-correlations are grouped into two clusters that are differentially modulated by behavior. A1–C2, Illustrative data from one recording (rat 1, recording day 2). A1–A3, Temporal coverage by each microstate during longer resting periods (left bars) and object interaction (right bars shown in light colors) indicates the dominant microstate during each behavioral state. B1, Dominant microstates (yellow) during longer stationary periods are shown across areas. All dominant microstates were grouped for further cross-correlation spectrum analysis. B2, Power spectra of cross-correlation of relevant microstates (shown in B1) during longer rest periods show a peak at ∼8 Hz (thick line shows the mean power spectrum). C1, Same as B1 for object interaction microstates. C2, Cross-correlation power spectra of dominant microstates (C1, yellow) during object interaction show an earlier peak at ∼2–4 Hz (thick line shows the mean power spectrum). D, Mean power spectra from all recordings show higher affinity to one of the spectral cluster centers (Fig. 6F1,F2), except very few spectra from object interaction microstate combinations that show relatively greater similarity to the cluster more populated by rest microstates.
Furthermore, power spectra associated with long resting periods were found to have a higher squared correlation with the 2–4 Hz cluster (r2 = 0.87 ± 0.22, mean + SD) compared with the 8–10 Hz cluster (r2 = 0.36 ± 0.32, p = 7.4e-05, two-sided Wilcoxon signed-rank test; Fig. 7D). On the other hand, the power spectra associated with object interaction correlated better with the higher-peak frequency cluster centers (r2 = 0.76 ± 0.24) compared with the low-peak frequency cluster (r2 = 0.50 ± 0.32, p = 0.024, two-sided Wilcoxon signed-rank test) in terms of squared correlation (Fig. 7D). Thus, each behavioral state exhibited a distinct power spectrum profile. The two clusters of power spectral profiles (Fig. 6F) were modulated by object exploration and movement.
In summary, the dependence of microstate properties on two contrasting behavioral states (object exploration and longer stationary periods) is supported by both the significant changes in microstate properties (Fig. 5) and the distinct nature of interaction patterns (cross-correlation spectra) across regions. In fact, these two types of interaction indicators relate to the two behavioral states. This suggests the specificity of the underlying neural mechanism of the two behavioral states (Figs. 8, 9).
Distributions of microstate properties in all three regions in all recording sessions provide insights on temporal and spatial features of microstates (Table 1). A, Temporal coverage fraction shows temporal area covered by each microstate. B, Occurrence rate shows that microstates in all three areas appeared approximately three to four times per second on average. C, The top four microstates exhibited a high global explained variance (60–80%) in all three brain regions. D, The top four microstates exhibited high spatial correlation (Pearson) with the electrode topography (∼0.5). E, State duration of microstates followed power law distributions in PFC and STR, and had a characteristic peak of ∼100 ms in the VTA.
LFP microstates are a small set of recurring spatial patterns of LFP signals recorded in a single brain region [here PFC (navy), STR (maroon), VTA (orange)]. These microstates, originating from different brain regions coordinate during different behavioral states. This is exhibited through their distinct behavioral state-specific cross-spectral density profiles. During the stationary state, a subset of microstates coordinates with 8–10 Hz cross-spectral density (csd) profile (pink), while during object interaction another subset of microstates shows a peak of ∼2–4 Hz in their csd profile (gray). Data from one recording session (rat 1, recording day 2).
Discussion
We extended a spatiotemporal analysis method developed for EEG data to characterize the dynamics of mesoscale signals (LFP) in multiple brain regions, and to search for states of activation and their interactions, during novelty seeking in awake behaving rats. We identified microstates analogous to EEG microstates, which were physiologically interpretable and modulated by the behavioral state. This establishes a novel link between network-level activations and dynamics and concurrent behavioral states, and may facilitate bridging human EEG and rodent invasive research.
“States” and scales in brain activity
The concept of a “brain state” appears in multiple domains of the neuroscience literature, and spans multiple spatial scales and levels of analysis. For example, the “resting state” in fMRI (Biswal et al., 1995; Raichle et al., 2001) refers to coordination across widespread brain regions during rest. Large-spatial scale states have also been defined based on temporally precise measurements such as EEG and MEG. Examples include defining entropy as the basis of states (Kannathal et al., 2005), using the embedded dimension in a state-space interaction to define cortical synchrony (Carmeli et al., 2005); hidden Markov models (Baker et al., 2014; Hunyadi et al., 2019); and the clustering-based approach of EEG data segmentation into quasi-stable states (Lehmann et al., 1987; Lehmann, 1989; Michel et al., 1999; Michel and Koenig, 2018).
At a finer spatial scale, neural states have been defined based on population spiking activity. For example, “up states” and “down states” in the context of depolarization and hyperpolarization (Steriade et al., 1993), the concept of “vocabulary” in population spike patterns (Luczak et al., 2009), propagating wave patterns (Townsend and Gong, 2018), and several others (Greenberg et al., 2008; Poulet and Petersen, 2008; Curto et al., 2009; Okun et al., 2010).
Our approach was to adapt methods developed in the EEG literature to multichannel LFP recordings. Although it is difficult for LFP microstates to directly connect with human EEG microstates, many features of the identified LFP microstates appear to be consistent with EEG microstate findings (Table 1; Michel and Koenig, 2018), including the number of microstates, their temporal characteristics, and their relation to ongoing behavior. Interestingly, state durations in the range of 50–100 ms were observed in our data, in the EEG microstates literature (Michel et al., 1999), and in population spike recordings (Luczak et al., 2009). One of the findings of this study is the presence of zero-and near-zero lag synchronization across the brain areas via microstate activations in combination with statistically robust nonzero correlation magnitudes (Fig. 6). An open question in EEG microstate research is the amount of true synchronization because of volume conduction (Michel and Koenig, 2018), as the volume conduction can elevate cross-region coherence (Srinivasan et al., 2007). The cross-correlation analysis in the current study presents the novel finding that zero-lag neuronal synchronization within and between brain areas could contribute to EEG microstates, as the multiregion LFP recordings remain largely unaffected by volume conduction. Furthermore, this hints at a possible common link underlying observable brain states across different spatial scales.
Microstates and dimensionality reduction approaches
LFP microstates are identified by a data-driven dimensionality reduction approach that is based on long, continuous recordings and yet can be linked to transient behaviors. We found that four microstates can explain 60–80% of the LFP signals. Overall, this is consistent with other dimensionality reduction approaches, for example independent component analysis of LFP in the hippocampal CA1 (Makarov et al., 2010), or principal component analysis in rat barrel cortex (Einevoll et al., 2007), monkey PFC (Machens et al., 2010), and rat orbitofrontal cortex (Kobak et al., 2016). It is remarkable that the conclusion of three to six large components persists over many different approaches. It seems that multichannel LFP data occupy low-dimensional subspaces, and that different dimensionality reduction methods may identify different basis vectors to span a limited set of the most relevant subspaces. The EEG microstate method uses clustering of the empirical data, which facilitates physiological and behavioral interpretations (Khanna et al., 2015; Michel and Koenig, 2018).
LFP microstates exhibit both stability and variability over time
The reproducibility of the microstate maps across 5–9 d suggests that the microstates are stable over time within animals, perhaps reflecting the activation of circuits that are anatomically patterned. On the other hand, some microstate maps were not in the top four components in all recording sessions. These might have a lower signal-to-noise ratio and thus were relegated to lower-ranked microstates, or they might reflect state-specific functional networks that were less constrained by stable anatomic connectivity.
EEG microstate studies have shown that four “standard” microstate maps are observed consistently over individuals and across many studies (Michel and Koenig, 2018). In our LFP data, only some of the maps were topographically reproducible across animals. However, this apparent discrepancy must be interpreted with caution: the LFP electrodes are much smaller than EEG electrodes, and therefore even small (submillimeter) differences in the placement of probes can cause topographical discrepancies across animals. Furthermore, EEG microstates reflect a macroscopic phenomenon (Betzel and Bassett, 2017), and the relationship between macroscopic and mesoscopic dynamics like LFP is nontrivial (Nunez and Srinivasan, 2006; Alishbayli et al., 2019). Generally, LFPs can be more sensitive to local effects compared with EEG signals, which are spatiotemporally filtered manifestations of a multitude of LFP signals (Buzsáki et al., 2012).
LFP microstates reflect mesoscale rhythmic coordination across brain regions
We identified a novel manifestation of rhythmic synchronization between mesoscopic-scale networks in spatially distinct brain regions (Fig. 9). Remarkably, the interactions exhibited two prominent peaks in the lower-frequency spectrum that mapped onto different aspects of behavior. The ∼2 Hz peak was more strongly modulated by interacting with the novel object, whereas the ∼8 Hz peak was more strongly modulated by movement. This suggests that partly overlapping networks use spectral multiplexing to coordinate inter-regional interaction during novelty using two distinct routes of mesoscopic information exchange (Akam and Kullmann, 2014).
The PFC is crucial for higher cognition (Uylings et al., 2003; Laubach et al., 2018) and is involved in executive functions, including learning and attention (Passetti et al., 2002), decision-making (De Bruin et al., 2000; Sul et al., 2010), working memory (Bekinschtein and Weisstaub, 2014), and inhibitory response control (Brown and Bowman, 2002; Hardung et al., 2017). The PFC also plays a crucial role in communicating with other cortical and subcortical regions to implement executive control (Pasupathy and Miller, 2005; Antzoulatos and Miller, 2014; Feingold et al., 2015; Zhang et al., 2016; Gao et al., 2017). Interactions among the PFC, VTA, and striatum are also crucial, particularly for learning and working memory, through dopaminergic pathways (Karreman and Moghaddam, 1996; Zahrt et al., 1997; Peters et al., 2004; Bentivoglio and Morelli, 2005; Björklund and Dunnett, 2007; Gao et al., 2007; Fujisawa and Buzsáki, 2011; Ott et al., 2014).
Although traditional analyses often involve synchronization between a selected pair of electrodes, our multichannel recordings allowed us to examine synchronization at a larger spatial scale. Because the microstates reflect a weighted combination of data from all channels, these can be considered as functional patterns that span multiple millimeters and reflect coordinated activity across possibly thousands of cells. Because the microstates are defined independently in each region, the synchronization results were not biased by selecting electrode pairs that exhibit a particular pattern of correlation. One major advantage of this method is the ability to decompose multiple partially spatially overlapping networks. Indeed, it is possible that traditional electrode-based synchronization methods mix neurally and cognitively distinct interactions because a single electrode measures activity from multiple networks (Cunningham and Yu, 2014; Cohen, 2021).
Implications for understanding novelty seeking
The experimental paradigm used here is a variant of the commonly used object recognition test (Aggleton, 1985). Novelty seeking has been shown to evoke exploratory behavior (Berlyne, 1950). The reduction of exploratory behavior during the repetition of a previously novel object is taken as a measure of single-trial memory formation. The role and interplay of PFC, STR, and VTA is well known for memory and motivation (Karreman and Moghaddam, 1996; Palmiter, 2008; Liljeholm and O'Doherty, 2012; Antzoulatos and Miller, 2014; Bekinschtein and Weisstaub, 2014; Morales and Margolis, 2017). The access to these brain regions makes this experimental paradigm suitable for exploring the neural states underlying the novelty-seeking behavior. Moreover, the prolonged stationary behavioral state, as the trials progressed, allowed us to investigate and compare the neural activity during novelty seeking against the “background activity.” In this study, we propose a role for mesoscale microstates in novelty-seeking behavior, which can be attributed to the activation of different networks during novelty processing.
Limitations and future directions
First, although LFP microstates have a clear methodological basis, their precise involvement in cognition remains to be determined. EEG microstates have been theorized to reflect building blocks of brain-wide events (Lehmann, 1990). LFP microstates may also be a mesoscopic equivalent of default-mode networks (Biswal et al., 1995). Second, the true number of microstates and their overall shape (the spatial structure) during different sensory stimulation and spontaneous conditions is difficult to determine even with multielectrode arrays. Future studies with larger and denser arrays will help to map the spatial reach of microstates. Third, the relationship between intermicrostate synchronization, as we observed it, with other manifestations of neural synchronization—and their relation to various behaviors—remain to be explored. Third, it remains unclear whether the rhythmic coordination that we observed reflects long-range excitation or inhibition, although these patterns do seem to play a meaningful role in relation to behavior. Finally, the exact relationship between LFP microstates and microstates at smaller (e.g., spiking) and larger (EEG) spatial scales remains to be investigated.
Footnotes
M.X.C. is funded by European Research Council Timeframe Starting Grant 638589 awarded to M.X.C. A.M. is funded by a joint doctoral grant awarded to M.X.C. and B.E. by the Donders Center of Neuroscience. B.E. was supported by European Commission Marie Curie Grant 660328, NWO Vidi Grant 016.189.052, and NWO ALW Open Grant ALWOP.346. We thank Dr. Paul Anderson for contributions to the surgeries and data collection process.
The authors declare no competing financial interests.
- Correspondence should be addressed to Ashutosh Mishra at a.mishra{at}donders.ru.nl
