Dynamic reconfiguration of frequency-specific cortical coactivation patterns during psychedelic and anesthetized states induced by ketamine.

Recent neuroimaging studies have demonstrated that spontaneous brain activity exhibits rich spatiotemporal structure that can be characterized as the exploration of a repertoire of spatially distributed patterns that recur over time. The repertoire of brain states may reflect the capacity for consciousness, since general anesthetics suppress and psychedelic drugs enhance such dynamics. However, the modulation of brain activity repertoire across varying states of consciousness has not yet been studied in a systematic and unified framework. As a unique drug that has both psychedelic and anesthetic properties depending on the dose, ketamine offers an opportunity to examine brain reconfiguration dynamics along a continuum of consciousness. Here we investigated the dynamic organization of cortical activity during wakefulness and during altered states of consciousness induced by different doses of ketamine. Through k-means clustering analysis of the envelope data of source-localized electroencephalographic (EEG) signals, we identified a set of recurring states that represent frequency-specific spatial coactivation patterns. We quantified the effect of ketamine on individual brain states in terms of fractional occupancy and transition probabilities and found that ketamine anesthesia tends to shift the configuration toward brain states with low spatial variability. Furthermore, by assessing the temporal dynamics of the occurrence and transitions of brain states, we showed that subanesthetic ketamine is associated with a richer repertoire, while anesthetic ketamine induces dynamic changes in brain state organization, with the repertoire richness evolving from a reduced level to one comparable to that of normal wakefulness before recovery of consciousness. These results provide a novel description of ketamine's modulation of the dynamic configuration of cortical activity and advance understanding of the neurophysiological mechanism of ketamine in terms of the spatial, temporal, and spectral structures of underlying whole-brain dynamics.

Introduction

Unlike canonical anesthetics and psychedelics, ketamine is thought to work primarily by antagonizing N-methyl-D-aspartate receptors and HCN1 channels (Zhou et al., 2013), with psychedelic and anesthetic properties in a dose-dependent manner (Corssen and Domino, 1966; Domino et al., 1965). A large body of work has investigated the neurophysiologic effects of ketamine at anesthetic (Akeju et al., 2016; Blain-Moraes et al., 2014; Bonhomme et al., 2016; Colombo et al., 2019; Lee et al., 2013; Purdon et al., 2015; Vlisides et al., 2017) and subanesthetic (de la Salle et al., 2016; Driesen et al., 2013; Höflich et al., 2015; Liao et al., 2012; Muthukumaraswamy et al., 2015; Niesters et al., 2012; Pallavicini et al., 2019; Rivolta et al., 2015; Scheidegger et al., 2012; Vlisides et al., 2018) doses. Specifically, ketamine has been demonstrated to induce frequency-specific changes in electrophysiological oscillatory activity and network-specific changes in hemodynamic responses. However, these findings are largely derived from the analysis of time-averaged electrophysiological or functional magnetic resonance imaging (fMRI) data that span at least a few minutes, thus the resultant measure disregards the vast amount of time-varying information that is present in the data. Several studies have incorporated the temporal dynamics. Using Lempel-Ziv complexity (LZC) (Lempel and Ziv, 1976; Welch, 1984; Ziv and Lempel, 1978) as a measure of neural signal diversity, ketamine anesthesia appears to maintain the spatiotemporal LZC of cortical activation evoked by transcranial magnetic stimulation as observed during normal wakefulness, while a psychedelic dose of ketamine is associated with increased single-channel temporal and spatiotemporal LZC in spontaneous magnetoencephalographic (MEG) signals (Schartner et al., 2017). We recently demonstrated the dose-dependent effect of ketamine on spatiotemporal complexity in spontaneous electroencephalographic (EEG) signals (Li and Mashour, 2019) in a single cohort of healthy volunteers. LZC assesses the spatial diversity in a short time period of a few seconds and reduces the spatiotemporal pattern to a unidimensional measure, thus providing a simple and powerful way of indexing ketamine-induced altered states of consciousness. However, how the different dosing of ketamine modulates the complex spatiotemporal dynamics over a longer time scale has not yet been explored.

The resting brain is an inherently dynamic system. Recent neuroimaging studies have demonstrated that spontaneous brain activity exhibits rich spatiotemporal structure that can be characterized as the exploration of a repertoire of spatially distributed patterns, known as dynamic brain states, that recur over time (Baker et al., 2014; Gutierrez-Barragan et al., 2019; Vidaurre et al., 2018). These dynamics may be a manifestation of covert cognitive information processing (Cabral et al., 2017; Gonzalez-Castillo et al., 2019), and reflect the information capacity of the brain that is thought to be central to consciousness (Cavanna et al., 2017; Tononi et al., 2016). Such explorative dynamics have also been observed in both humans and animals during deep sleep (Brodbeck et al., 2012; Houldin et al., 2018; Stevner et al., 2019), anesthesia (Hudson et al., 2014; Li et al., 2019; Varley et al., 2021; Vlisides et al., 2019; Wenzel et al., 2019) or disorders of consciousness (Demertzi et al., 2019; Khanmohammadi et al., 2018). Specifically, despite differing molecular mechanisms, general anesthetics such as propofol, ketamine, sevoflurane, and isoflurane have been reported to suppress the flexible brain configurations during wakefulness to a rigid brain configuration that is predominantly shaped by brain anatomy (Barttfeld et al., 2015; Ma et al., 2017; Uhrig et al., 2018; Zhang et al., 2019). Conversely, psychedelic drugs typically induce the sustained occurrence of rich phenomenology and have been associated with an expanded repertoire of brain states (Atasoy et al., 2018; Tagliazucchi et al., 2014). Despite these encouraging findings, the modulation of the dynamic repertoire across varying states of consciousness has not been systematically assessed in a unified manner. As a unique drug with both psychedelic and anesthetic properties depending on dose, ketamine provides a pharmacological tool to characterize brain dynamics from the restricted brain configurations associated with general anesthesia to the expanded repertoire of the psychedelic state.

Ketamine is a fast acting anesthetic, with a distribution half-life of only 10–15 min (Li and Vlisides, 2016). We and others have shown that an anesthetic dose of ketamine induces fast transient changes in EEG signals that is characterized by alternating slow delta oscillations and gamma oscillations (Akeju et al., 2016; Garwood et al., 2021; Li and Mashour, 2019; Schwartz et al., 1974). The objective of the current study was to use high-density EEG signals that directly measure neural activity at high temporal resolution to investigate the dynamic brain states and their reconfigurations during altered states of consciousness induced by ketamine. Given the frequency-dependent effect of ketamine on cortical activity (Brito et al., 2020; Li and Mashour, 2019), we began by extracting the band-limited envelope of source-localized signals within the canonically defined frequency bands (delta, theta, alpha, beta, and gamma bands). With the envelope data aggregated over subjects, we employed an unsupervised clustering technique to identify recurrent brain states of frequency-specific spatial coactivation patterns and analyze the temporal dynamics of the occupancy and transitions of these brain states during wakefulness as well as during different doses of ketamine. We hypothesized that the repertoire, or diversity, of the brain states would increase with a subanesthetic dose of ketamine (associated with psychedelic experiences) and decrease with an anesthetic dose of ketamine. Furthermore, based on the pharmacokinetics and phenomenology of ketamine, as well as our previous study (Li and Mashour, 2019), we hypothesized that the repertoire of brain states would not be consistently depressed but rather fluctuate during ketamine anesthesia, with the possibility of returning to baseline level before recovery of responsiveness.

Materials and methods

In this study, we re-analyzed EEG data in 15 healthy volunteers who received subanesthetic and anesthetic dosing of ketamine administration in prior studies, in which we examined the dose-dependent effect of ketamine on sensor-level EEG oscillatory/connectivity patterns (Vlisides et al., 2017) and spatiotemporal complexity (Li and Mashour, 2019), as well as the time-averaged changes of source-level spectral power during subanesthetic ketamine, at a dose that induces psychedelic experiences (Vlisides et al., 2018). Unlike previous studies, the current study investigated the dynamic repertoire of cortical coactivation patterns at the source level and its modulation by different dosing of ketamine.

Experimental protocol

Details of the study population and experimental protocol can be found in prior publication (Vlisides et al., 2017, 2018). In this study, we analyzed the source-localized EEG during the following four periods. The first period (baseline wakefulness) was 5 min of rest with eyes closed. The second period was the last 5 min eyes-closed epoch of a continuous infusion of subanesthetic (0.5 mg·kg−1 total) ketamine over 40 min. Only the last 5 min were analyzed because, based on a pilot study with the same subanesthetic dosing regimen, we assume the ketamine infusion reached pharmacological steady-state conditions by that epoch. The third period was the anesthetic period from loss of consciousness to the recovery of consciousness (defined by loss or recovery of responsiveness, with acknowledgment that disconnected subjective experience could occur during ketamine anesthesia), with the duration of 10.1 ± 3.4 [mean ± standard deviation (SD)] min across subjects. General anesthesia was initiated with a 1.5 mg·kg−1 bolus dose of ketamine. The fourth period was 5 min of eyes-closed rest after return of consciousness. We assume the baseline, subanesthesia, and recovery periods to be steady state, except for the anesthetic period because the bolus injection may produce non-stationary effect-site drug concentrations.

EEG acquisition and preprocessing

The EEG data were acquired with 128-channel EGI Hydrocel Nets (Eugene, OR, USA) digitized continuously at 500 Hz with a vertex reference; channel impedances were kept below 50 kΩ as recommended by the manufacturer. The raw EEG signals were exported into MAT-LAB (version 2017a; MathWorks, Inc., Natick, MA), and downsampled to 250 Hz. Electrodes on the lowest parts of the face and head were removed, and the remaining 109-channel recordings on the scalp were preprocessed as follows. First, the signals were detrended using a local linear regression method with a 10-s window at a step size of 5-s in Chronux analysis toolbox (Mitra and Bokil, 2007), which were then low-pass filtered at 50 Hz via a 5-order Butterworth filter using a zero-phase forward and reverse algorithm. Second, bad channels were detected by visual inspection and interpolated by using the spherical spline interpolation method in EEGLAB toolbox (Delorme and Makeig, 2004), and then the EEG signals were re-referenced to the average reference. Third, noisy data segments were automatically detected and rejected in a step-wise manner (Li et al., 2019). The signals were bandpass filtered at 0.5–45 Hz via a 4-order Butterworth filter (only for noisy segment detection) and then segmented into 1-s epochs. The 1-s data were rejected if (1) the average amplitude was > 3 ·average amplitude (or its SD > 2 ·SD value) of the whole recording, and (2) this was present in at least 4 channels. This step was performed for the baseline, subanesthesia, anesthesia, and recovery periods separately, and 14.1 ± 12.0%, 20.9 ± 12.3%, 3.4 ± 4.1% and 14.4 ± 9.4% (mean ± SD) of the data were rejected, respectively. The remaining non-bandpassed data were used for further analysis, with the data length of 4.5 ± 0.6, 4.0 ± 0.6, 9.7 ± 3.4 and 4.5 ± 0.5 min for the four studied periods.

For the preprocessed EEG signals, the analysis pipeline is illustrated in Fig.1, with each processing step described as follows.

Fig. 1.

Schematic overview of the analysis pipeline. The 128-channel EEG data of 15 subjects during baseline, subanesthetic and anesthetic ketamine, as well as recovery period, were preprocessed in sensor space and concatenated for the analysis. Cortical sources were estimated using the weighted minimum norm estimation (wMNE) method, followed by the extraction of regional time series by averaging the cortical activity across voxels within each of the 100 regions based on the Yeo atlas. Each region can be matched to one of the seven resting-state networks (RSNs): FPN, frontoparietal network, DMN, default mode network, DAT, dorsal attention network, LIM, limbic network, VAT, ventral attention network, SOM sensorimotor network, and VIS, visual network. Band-limited envelope data were derived from Hilbert transformation and aggregated across all studied periods and subjects. K-means clustering algorithm was applied to generate a discrete set of brain states and the state time series for each studied period and subject. The spectral and spatial properties of the brain states, the temporal dynamics of state time series and their alterations with different dosing of ketamine were subsequently characterized.

Source reconstruction

The source reconstruction was performed with weighted minimum norm estimate method (wMNE) (Baillet et al., 2001) in Brainstorm (Tadel et al., 2011). This method imposes l2-norm constraints on the source distribution and solves the inverse problem by choosing the configuration with the minimum energy, with the incorporation of a depth-weighting factor compensating for the tendency of classical MNE to favor weak and surface sources. The wMNE method has established characteristics in statistics and mathematics, which is less sensitive to in-accuracy of the forward model and has been a widely used method for EEG source imaging, especially when no individual structural MRI anatomy and digitized electrode positions were available (Hassan and Wendling, 2018; He et al., 2018).

First, a realistic head model was derived from the FSaverage template, which is the average of 40 subjects (Fischl et al., 1999), by using a symmetric boundary element method in OpenEEG (Gramfort et al., 2010; Kybic et al., 2005). The head model includes three shells representing scalp, skull, and brain, by assuming the electrical conductivity within each layer being homogeneous, which is 0.33, 0.0042, 0.33 S/m, respectively. The noise covariance matrix was set as an identity matrix by assuming equal and unit variance of noise over all sensors. The source solution space was constrained to the cerebral cortex, with the orientations constrained to be perpendicular to the cortical surface. The current density distribution was estimated on a three-dimensional grid of 15,002 vertices, which was then segmented into 100 parcels or regions-of-interest based on the Yeo atlas (Schaefer et al., 2017; Thomas Yeo et al., 2011). Activity of each region was obtained by averaging the current source density of all voxels within that region. This parcellation resulted in regions with high homogeneity, where each region is matched to one of the seven resting-state networks (RSNs): frontoparietal network (FPN), default mode network (DMN), dorsal attention network (DAT), ventral attention network (VAT), sensorimotor network (SOM), visual network (VIS), and limbic network (LIM).

Supplementary analyses were performed using alternative approaches in terms of template anatomy (ICBM152 head (Fonov et al., 2011)), inverse estimation (standardized low-resolution electromagnetic tomography [sLORETA] (Pascual-Marqui, 2002)), and parcellation method (Desikan-Killiany atlas (Desikan et al., 2006)).

Identification of brain states

To reveal the frequency-specific spatial coactivating patterns, the reconstructed source data were band-pass filtered into five frequency bands: delta (1–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), beta (13–25 Hz), and gamma (25–45 Hz). For each frequency band, the envelope time series was derived by computing the magnitude of the Hilbert transform of the band-limited signal for each cortical region, which was then demeaned and down-sampled to 40 Hz by temporally averaging within sliding windows of 100 ms and 75% overlapping (Baker et al., 2014), equivalent to low-pass filtering the envelope time series below 10 Hz. We corroborated the window selection strategy by performing supplementary analysis with window length ranged from 50 to 1000 ms (i.e., low pass filtering below 20 to 1 Hz).

To assess the temporal progression of spatial patterns of cortical activation and its correlates with ketamine, the envelope data were normalized by the standard deviation of each region (to ensure that the sought dynamical effects are not confounded with inter-state and inter-individual variability, and that each state and individual is given a similar weight), and then concatenated temporally across all studied periods (baseline, subanesthesia, anesthesia, and recovery) and subjects, yielding a matrix X with m features and N observations. Here, m equals 500, which is the number of frequency bands by the number of spatial regions, and N is the number of time points of all 4 periods from 15 subjects (in this case, 838,084). The concatenated matrix was subject to k-means algorithm as implemented in Matlab, which partitioned the N observations into N clusters that have minimized within-cluster distance (as defined by one minus Pearson’s correlation coefficient between two observations). We used the correlation metric by following previous fMRI studies on coactivation patterns (Chen et al., 2015; Cornblath et al., 2020; Gutierrez-Barragan et al., 2019; Huang et al., 2020; Liu et al., 2013; Liu and Duyn, 2013), and corroborated the selection by performing supplementary analysis with an alternative metric (one minus cosine similarity). The clustering procedure was repeated 100 times, each with a new set of initial centroids to avoid local minima problems. We varied the number of clusters N from 5 to 15 and selected the solution with 10 clusters that demonstrated most consistency across subjects (see Supplementary materials, where we also performed a series of tests to evaluate the clustering solutions).

With clustering analysis, the envelope matrix X was transformed into (1) a time sequence of labels with each element indicating the brain state (cluster) corresponding to each time point of N observations across studied periods and subjects, (2) a discrete set of brain states with distinct spectral and spatial properties. Each brain state can be represented by a vector of m features, which was generated by normalizing the mean activation pattern by the standard deviation of the observations that were assigned to that cluster. The resultant activation patterns were then mapped to the cortical surface that allowed for the comparisons between brain states and between frequency bands (within brain state). For each state, we assumed the spatial activation pattern is spectrally heterogeneous and defined the dominant frequency band in which the highest activation (positive envelope values) or deactivation (negative envelope values) occurs. In the dominant frequency band, we assessed the spatial variability of the activation pattern by (1) computing the standard deviation of the activation values across regions, and (2) applying Shannon entropy to the distribution of the activation values (using a histogram), with a higher entropy value indicating a broader distribution and thus a greater level of spatial diversity. To further compare the spatial patterns with resting state networks, we first generated a binary vector for each of the 7 RSNs (as described in Source reconstruction) with each element corresponding to each of the 100 cortical regions, which was set to 1 if the region belonged to the RSN or 0 otherwise. We then computed the cosine similarity (Cornblath et al., 2020) between the binary vectors and the spatial pattern in the dominant frequency band for each brain state.

Analysis of brain state dynamics

For each subject, the state time series represents the temporal progression of brain states during baseline, subanesthesia, anesthesia, and recovery. Across all studied periods, we first characterized the temporal properties relating to the occurrence of brain states, by quantifying the fractional occupancy, mean dwell time, and mean interval time for each brain state in each subject. The fractional occupancy is defined as the fraction of time spent in each brain state. The mean dwell time is defined as the average amount of time spent in each state before transitioning out of that state, while the mean interval time is the average amount of time spent between consecutive visits to a certain state (Baker et al., 2014). We then examined how the occupancy of brain states was modulated by ketamine and computed the fractional occupancy of the brain states in each of the studied periods, respectively. Furthermore, we applied Shannon entropy to the distribution of fractional occupancy across brain states (Demertzi et al., 2019), which ranged from 0 to 1, with the highest value indicating that it is equally probable to visit any of all possible states, while lower values indicating that it is more probable to visit some states than others.

We further investigated the effect of ketamine on state transitioning. We assumed the brain state time series to be a Markov chain (i.e., the state transition depends only on the current state; with the Markov property tested by following (von Wegner et al., 2017) in Supplementary materials) and computed the state transition probability matrix for each studied period in each subject. Transition probability measures the likelihood of the current state transiting to another state at a future time. Specifically, for the transition from state i to state j(i, j = 1, 2, …, N) the probability was estimated by counting the times of this transition divided by the total times of all possible state transitions (which is equal to the length of the state time series minus 1). In the resultant matrix, the elements on the diagonal line indicate the probability of staying in a certain state (referred to as persistence probabilities), while the off-diagonal elements indicate the probability of switching between two distinct states (referred to as transition probabilities), with all persistence and transition probabilities summing up to 1 (Li et al., 2019). Moreover, we focused on between-state transitions, and computed the state transition matrix from the retained state time series after removing state stays; the resultant matrix represents the transition probability for each pair of distinct states on off-diagonal elements only.

We used the following summary measures to quantify the state transition dynamics. First, the transition rate is defined as the number of state transitions between distinct states divided by the total time a subject spent in each of the studied periods (Li et al., 2019; Stevner et al., 2019), with lower (higher) values suggesting that the cortical activity is less (more) likely to transition to a distinct state, or equivalently, more (less) likely to be persistent or “sticky ” in the same state. Second, the entropy rate (Markov entropy) that measures the unpredictability of state transitions (Demertzi et al., 2019) is defined as where p (i, j) corresponds to the transition probability from state i to state j, and NST is the number of all state transitions. We computed the entropy associated with the complete state transition matrix (including both persistence and transition probabilities) and the reduced matrix with only state transition probabilities, respectively. A lower (higher) entropy value indicates a lesser (greater) amount of information encoded in the state transition dynamics, or equivalently, the brain state at time t + 1 is more (or less) predictable.

Statistical analysis

Linear mixed model (LMM) analysis was performed with IBM SPSS Statistics version 24.0 for Windows (IBM Corp. Armonk, NY) for the following comparisons. (1) To compare the brain states in terms of temporal properties, a LMM was fitted with the brain state as the fixed effect for fractional occupancy, mean dwell time and mean interval time. (2) To compare the occupancy of brain states across baseline, subanesthesia, anesthesia, and recovery, a LMM was fitted with the studied period as the fixed effect for each brain state individually. (3) To compare the summary measures for brain state dynamics across the studied periods, a LMM was fitted with the studied period as the fixed effect for the entropy associated with state occupancy, transition rate, entropy associated with both persistence and transition probabilities, and the entropy associated with only transition probabilities. For all models, we used restricted maximum likelihood estimation. We did not model repeated effects but included a random intercept specific for each subject. Student’s two-tailed paired t tests were used for post-hoc pairwise tests, and Bonferroni method was used to correct for multiple comparisons.

To test whether the state transitions we observed could occur by chance, we performed the following surrogate data analysis (Li et al., 2019). First, we generated N = 1000 surrogate time series by randomly shuffling the state time series for each subject, which permutated the temporal order of the state occurrence while keeping the occupancy of the states. With each surrogate time series, a state transition matrix was calculated, and the statistical significance for a state transition was obtained through the cumulative distribution function where α denotes the original transition probability and p (ℎ) denotes the estimated normal distribution by assuming normality of surrogate data. To further test the randomness of between-state transitions, we removed the state stays and computed the state transition matrix from the retained state time series. We assumed that it is equally probable for the transitions to occur between any pair of distinct states, and used right-tailed, Wilcoxon signed rank test to obtain the statistical significance for a state transition with probability higher than the averaged transition probability (i.e., 1/90, where 90 is the number of all possible state transitions) across subjects. For the statistically significant state transitions as detected above, we used two-tailed, Wilcoxon signed rank test to compare the state persistence and transition probabilities in baseline and each of the studied periods (subanesthesia, anesthesia, and recovery). The false discovery rate (FDR) approach was applied to correct for multiple comparisons. The statistical analyses described above were performed using MATLAB. All data sets were tested for normality of distribution by Lilliefors corrected Kolmogorov-Smirnov tests. A p value of < 0.05 was considered statistically significant.

Results

Frequency-specific cortical coactivation patterns identified by unsupervised clustering analysis

Spontaneous cortical activity exhibited temporally recurrent cofluctuations in the band-limited envelope time courses across distinct brain regions (Fig. S1). To extract the inherent recurring patterns, we used a k-means algorithm to partition the envelope data into ten clusters based on their spectral and spatial similarity. Each cluster represents a brain state that is characterized by distinct spectral (delta, theta, alpha, beta, and gamma) and spatial (mapped to 100 cortical regions) properties. The brain states were derived from all studied periods (baseline, subanesthesia, anesthesia, and recovery) in all subjects. Each subject had their own state time course, representing the brain state being active at each instant in each studied period. For each brain state, we computed the specificity defined as the likelihood of the state occurring in each studied period (Stevner et al., 2019) (Fig. 2 A) and, based on the difference in specificity between baseline and anesthesia period, the ten brain states were ordered, with their characteristic spectral and spatial patterns as shown in Fig. 2 B.

Fig. 2.

Brain states represent distinct spectral and spatial coactivation patterns in spontaneous cortical activity. (A) Ten brain states, identified by k-means clustering of the band-limited envelope data from 15 subjects during four periods (baseline, ketamine subanesthesia, ketamine anesthesia, and recovery), were sorted according to their specificity for baseline and anesthesia (i.e., the difference in probability of a brain state occurring within baseline and anesthesia periods). (B) Spectral and spatial distribution for each brain state, defined as the centroid of each cluster normalized by the standard deviation of all samples that were assigned to the same cluster. (C) Cosine similarity between each brain state in the dominant frequency band with seven canonical resting-state networks defined in the Yeo atlas, with positive (red) and negative (blue) values corresponding to high and low (above and below average) amplitude activity, respectively.

For each brain state, the cortical coactivation pattern is heterogeneous, both spectrally and spatially. States 1 and 2 had high specificity for the baseline period, showed strong negative spatial similarity (Pearson’s correlation coefficient r = −0.880) (Fig. S2A), and were associated with activation and deactivation of alpha activity in posterior regions including visual cortex and posterior cingulate cortex (Fig. 2 C). States 3 and 4 were characterized by high-amplitude activity in beta and alpha band, respectively, with the activation mapped to sensorimotor cortex but also spreading across regions involving frontoparietal network, dorsal attention network, and ventral attention network. States 5 and 7 demonstrated similar spatial distribution, with the greatest activation mapped to visual cortex, but primarily in beta and gamma bands, respectively. State 6 was dominated by activation in delta band, which was spatially mapped to cortical regions including prefrontal cortex and temporal cortex in the DMN. The remaining states exhibited high specificity for ketamine anesthesia, which showed simultaneous activation and deactivation in distinct frequency bands but spread globally across the whole cortex. State 8 showed global activation (deactivation) in theta (delta) band (Fig. S2B). States 9 and 10 were anticorrelated with each other (r = −0.779), with State 9 exhibiting global activation in delta and deactivation in gamma band, and vice versa for State 10. As compared to other states, the anesthesia-related states demonstrated lower spatial variability, as quantified by the standard deviation and entropy of the activation values across cortical regions (Fig. S2C).

Each brain state was visited by each subject, with the temporal characteristics of the state time series assessed by the fractional occupancy, mean dwell time, and mean interval time for each brain state in each subject (Fig. 3 A–C, with the descriptive statistics summarized in Table S1). Across the subjects, all brain states occupied 5.6–14.1% of the time on average (Fig. 3 A). The states that exhibited higher specificity for baseline (States 1 and 2) or anesthesia (States 8–10) occupied a larger proportion of time than other states (p < 0.05). The mean dwell time was between 84.4 and 202.2 ms on average across states (Fig. 3 B). The dwell times for the states with dominant frequency in beta and gamma bands (State 3,5,7) were much shorter in duration than other states (p < 0.001), but there was no significant frequency-related difference for the other states (p > 0.05). The state that most resembled the anterior component of the DMN (State 6) had a relatively long dwell time of 196.8 ± 43.1 [mean ± SD] ms, but also the maximum interval time between state visits, which was 3.1 ± 0.6 s and significantly greater than all other states (p < 0.001) (Fig. 3 C).

Fig. 3.

Temporal characteristics of the brain states and the effect of ketamine on state occurrence dynamics. (A) Fractional occupancy, defined as the fraction of time spent in each brain state. (B) The mean dwell time, defined as the average amount of time spent in each state before transitioning out of that state. (C) The mean interval time, defined as the average amount of time spent between consecutive visits to a certain state. (D) Fractional occupancy in each brain state across the studied period of baseline, subanesthesia, anesthesia, and recovery. The height of the colored bar and errorbar denote the mean and SD of the values across subjects. (E) Changes in the entropy that measures the distribution of fractional occupancy of brain states across the studied periods. For each subject, the anesthesia period was equally divided into four segments, and the state occupancy and entropy values were calculated for each segment individually. In A–C and E, the central line and edges on each box indicate the median and the interquartile range (IQR) of the values across the subjects, the whiskers extend to the most extreme values, and the outliers are marked as red crosses. In D,E, * indicates statistically significant difference relative to baseline (Bonferroni corrected p < 0.05, linear mixed model analysis).

It is worthwhile to note that the state characteristics (especially the numerical temporal features) must be interpreted in the context of the number of brain states partitioned and the window length for smoothing the envelope data prior to clustering analysis. We assessed the effect of varying the number of brain states (5 to 15 states) and the smoothing-window length (50 to 1000 ms). The increases in the number of brain states resulted in shorter between-state distance, smaller occupancy, and shorter dwell time, indicating that the addition of more states may split existing states and thus offer the characterization of brain dynamics at finer temporal scales (Fig. S3). As for the optimal number of brain states, different evaluation criteria yielded inconsistent recommendations (Fig. S4), in line with previous studies (Huang et al., 2020; Liu et al., 2013). The size of the smoothing window had little effect on the between-state distance and the fractional occupancy, while the mean dwell time became longer for the larger window size (Fig. S5). The employed solution with 10 states and 100 ms smoothing window was selected because it provided a tradeoff between a detailed representation of brain dynamics and a high degree of inter-subject consistency with an acceptable computational complexity.

Effects of ketamine on state occurrence dynamics

To assess the effects of ketamine on brain state dynamics, we first compared the fractional occupancy of each brain state across the studied periods of baseline, subanesthesia, anesthesia, and recovery (Fig. 3 D, with the detailed statistical analysis results summarized in Table S2). We observed dose-dependent decreases in the occupancy of posterior patterns (State 1 and 2), which were reversed at recovery but did not fully return to baseline values (p < 0.05). Similar changes were seen in State 3, while a small decrease was only detected in anesthesia for State 4 (p = 0.006), even though both states were associated with similar spatial patterns of activation. The two states that correlated with beta (State 5) or gamma (State 7) activation in visual cortex also behaved differently. The occupancy of State 5 was increased in subanesthesia and recovery (p < 0.01) but not in anesthesia, while that of State 7 was increased in a dose-dependent manner (p < 0.001). No statistically significant changes were detected in the occupancy of the anterior pattern (State 6). Consistent with high specificity for anesthesia, the global patterns (State 8–10) occupied the largest proportion of time during ketamine anesthesia (p < 0.001), which was followed by recovery (p < 0.05) and the subanesthesia period.

We further examined the entropy of the distribution of the occupancy across states. As demonstrated in Fig. S6, the temporal variability of the entropy values in anesthesia was more than fourfold higher than those in other periods, as might be expected due to the non-stationary effect-site drug concentrations associated with a bolus injection of ketamine. To reveal dynamic changes, we divided the state time series during anesthesia into four equal-length segments for each subject, and the entropy values of the state occupancy distribution were calculated and shown in Fig. 3 E, with the detailed statistical analysis results summarized in Table S3. We found increased entropy at subanesthesia as compared to baseline (p = 0.019, 0.065(0.023 – 0.108) for mean and 95% CI). As expected, the entropy values were dynamic during anesthesia; the second segment demonstrated the lowest entropy value, which was significantly lower than those in baseline (p < 0.001, −0.096(−0.139– −0.053)), and which gradually increased, then returned to the baseline level in the last segment. Like subanesthesia, the entropy values at recovery were significantly higher than those in baseline (p = 0.006, 0.075(0.032–0.119)), which was not unexpected considering the shift from anesthetic to subanesthetic effect-site drug concentrations.

Effects of ketamine on state transition dynamics

We investigated how cortical activity transitions among brain states during baseline wakefulness and how the state transitions were modulated by different doses of ketamine. The Markovianity test results suggest the validity of the first-order Markov assumption in the state time series (Table S4). Fig. 4 (A) shows the mean state transition matrix across subjects; this matrix represents the probability of state transitions for each pair of brain states in each studied period. When compared to random transitions by permutating the temporal order while keeping the occupancy of the states, the persistent probabilities were significantly higher for all states (p < 0.05), suggesting that cortical activity is more likely to be sticky in the same state than expected by chance. Next, we aggregated over all states and found the persistence probability in the same state was significantly higher than that of transitioning to a different state (Wilcoxon signed rank test p < 0.001, baseline: 84.0 ± 2.8% vs. 16.0 ± 2.8%, subanesthesia: 80.3 ± 3.1% vs. 19.7 ± 3.1%, anesthesia: 83.0 ± 3.4% vs. 17.0 ± 3.4%, recovery: 79.3 ± 3.0% vs. 20.7 ± 3.0% [mean ± SD]). Furthermore, we focused on between-state transitions from the state time series after removing state stays; the mean transition matrices are shown in Fig. 4 (B). The transitions between states were not evenly distributed and some state transitions occurred more frequently than others (p < 0.05). Collectively, these results demonstrate that, irrespective of the presence of ketamine or not, it is more probable to stay in a certain state than switching to another state, but when a switch occurs, brain states transition in a structured (i.e., not random) way.

Fig. 4.

The effects of ketamine on state transition dynamics. (A) Group average state transition matrix for each of the studied periods, with each off-diagonal element indicating the probability of transitioning from any state in each row to another state in the given column, while the element on the diagonal line indicates the probability of staying in a certain state. The elements with + indicate the state transitions with the probability statistically higher than those of random transitioning by permutating the temporal order while keeping the occupancy of the states (FDR-adjusted p < 0.05). (B) Group average state transition matrix from the retained state time series after removing the state stays. The elements with + indicate the state transitions with the probability statistically higher than average transition rate (1/90 in this study) (FDR-adjusted p < 0.05, Wilcoxon signed rank test). (C) Changes in state persistence and transition probabilities in subanesthesia, anesthesia, and recovery relative to baseline. Each node indicates a brain state, with different colors denoting statistically higher (red), lower (blue), or no changes (black) in the probability of staying in that state (FDR-adjusted p < 0.05, Wilcoxon signed rank test). The directed arrows in red (blue) denoted statistically higher (lower) probability for that transition (FDR-adjusted p < 0.05, Wilcoxon signed rank test). (D) Changes of transition rate. (E,F) Changes in entropy values associated with state persistence and transition probabilities (E) or state transition probabilities only (F). In D–F, * indicates statistically significant difference relative to baseline (Bonferroni corrected p < 0.05, linear mixed model analysis).

To assess the effect of ketamine on state transitions, Fig. 4 (C) shows statistically significant changes in state persistence and transition probabilities in subanesthesia, anesthesia, and recovery as compared to baseline period. In subanesthesia, we found reduced probabilities of staying in States 1 and 3 as well as transitioning between State 3 and State 1 or State 2, and increased probabilities of staying in States 5, 7 and 10 as well as transitioning between State 5 and State 7 and from State 7 to State 10 (p < 0.05). In accordance with the state specificity, an anesthetic dose of ketamine resulted in an overall decrease in the likelihood of staying in States 1–4, as well as transitioning into/out of these states; conversely, it became more probable to stay in States 7–10, as well as transition into/out of these states (p < 0.05). At recovery, the state transitions did not return to baseline level, in that they remained low for the persistence and transition probabilities involving States 1–3, and high for the persistence and switch probabilities involving States 5, 7, 9 and 10 (p < 0.05).

Apart from individual persistence and transition probabilities, we evaluated the effect of ketamine on state transition dynamics (Fig. 4 D–F, with the detailed statistical analysis results summarized in Table S3). The transition rate, which measures the percentage of time transitioning to a different state, increased at subanesthesia (p = 0.006, 3.7 (1.5– 5.9)% for mean and 95% CI) and recovery (p < 0.001, 4.6 (2.4–6.8)%) as compared to baseline (Fig. 4 D). The transition rate was not static during anesthesia, which increased gradually from the first to the fourth segment, with the values in the fourth segment significantly higher relative to baseline (p = 0.04, 3.1 (0.9–5.2)%). Second, the Markov entropy, which measures the unpredictability of the entire state transition matrix (including both state persistence and transitions), exhibited similar changes as observed in the entropy associated with state occupancy (Fig. 4 E). Relative to baseline, the entropy values increased at subanesthesia (p = 0.005, 0.064 (0.027–0.100)) and recovery (p < 0.001, 0.076 (0.039–0.113)), but were dynamic during anesthesia, with the values in the second segment significantly lower (p = 0.01, −0.057(−0.093—0.021)) and those in the other segments comparable to baseline. Third, the entropy associated with only transition probabilities showed a small but not statistically significant increase at subanesthesia and recovery (Fig. 4 F). The intra-anesthetic changes were also dynamic, with the entropy values significantly lower than those in baseline for all four segments (p < 0.05). In addition, as compared to the entropy for the entire state transition matrix (Fig. 4 E), this measure exhibited an overall increase, with the median entropy value from 0.531 to 0.722 to 0.891–0.940, which is expected because the state transitioning is less predictable after excluding state stays. Taken together, these results demonstrated that an anesthetic dose of ketamine induced dynamic changes in state transitioning, while subanesthetic ketamine was associated with increases in the amount of state transitions and the degree of unpredictability.

Additionally, we obtained confirmatory evidence that demonstrated similar dose-dependent changes in state occurrence and transition dynamics despite different strategies in the selection of smoothing-window length for envelope data and the number of brain states in the clustering analysis (Fig. S7), as well as alternative approaches of template anatomy and parcellation method (Fig. S8), inverse estimation method (Fig. S9), distance metric (one minus cosine similarity) (Fig. S10), and clustering algorithm (k-medoids) (Fig. S11).

Discussion

In this study, we investigated the dynamic organization of spontaneous cortical activity during wakefulness, subanesthetic ketamine associated with psychedelic effects (Vlisides et al., 2018), and ketamine anesthesia. Based on source-localized EEG signals, we identified a set of recurring states that represent frequency-specific spatial coactivation patterns. We assessed the effect of ketamine on the fractional occupancy and transition probabilities of individual brain states and found ketamine anesthesia tends to shift the configuration toward states with low spatial variability. Importantly, by quantifying the temporal dynamics of the occurrence and transitions of brain states, we showed a dose-dependent effect of ketamine on the richness of the dynamic repertoire. Subanesthetic ketamine is associated with a richer repertoire, while ketamine anesthesia induces dynamic changes in brain state organization, with the repertoire richness evolving from a reduced level to one comparable to that during normal wakefulness before recovery of consciousness. Therefore, ketamine induces brain activity configurations that range across the spectrum of conscious experience, making it a unique pharmacological tool to study consciousness. Furthermore, this study connects prior work on brain state repertoire using fMRI (slower dynamics) with studies of neurophysiologic complexity using EEG and MEG (faster dynamics).

The studies of whole-brain functional network organization have been primarily based on fMRI, which has provided consistent evidence for a set of resting state networks consisting of groups of brain regions that display statistically similar temporal fluctuations. An important limitation of fMRI is the slow temporal dynamics of the BOLD signals, with a time scale of seconds precluding a direct investigation of the rich temporal dynamics of the underlying neuronal activity. A recent whole-brain computational modeling analysis found that the most relevant time scale for brain processing is around 200 ms, suggesting that electrophysiological data may be more relevant for discovering spatiotemporal dynamics across the whole brain (Deco et al., 2019). Recent studies using complementary modalities have linked the hemodynamic measures of brain activity to resting-state electrophysiological recordings, and have shown spatial similarity in the covariation patterns between fMRI and envelopes of band-limited electrophysiological signals, in a frequency-dependent manner (Deco et al., 2017; Hipp et al., 2012; Kucyi et al., 2018; Siems et al., 2016; Wens et al., 2019). Motivated by this, as well as the frequency-dependent effect of ketamine on cortical activity (Li and Mashour, 2019), we extracted the band-limited envelope of source-localized EEG signals within the canonically defined frequency bands, and characterized the dynamic organization of cortical activity during wakefulness and during altered states of consciousness induced by ketamine on a time scale of a few hundreds of milliseconds (Fig. 3 B). Supplementary analysis demonstrated that the observed brain reconfiguration dynamics could not be fully revealed at a temporal resolution comparable to that provided by fMRI (Fig. S7A).

Through unsupervised clustering analysis of aggregated envelope data from wakefulness and different doses of ketamine, we derived a set of dynamic states that represent frequency-specific spatial coactivation patterns. During wakefulness, we identified a pair of anticorrelated states (States 1 and 2) characterized by activation and deactivation of alpha activity in posterior region including posterior cingulate cortex in the DMN, and another state (State 6) charactered by delta activation in anterior region including prefrontal cortex and temporal cortex in the DMN (Fig. 2 B). This subdivision of the DMN with characteristic spectral properties is consistent with the higher-order cognitive networks inferred from hidden Markov modeling analysis of spontaneous MEG signals at rest (Vidaurre et al., 2018). However, beside the posterior cingulate cortex, States 1 and 2 also exhibited activation and deactivation in visual cortex, which is not separable even when we increased the number of derived states (Fig. S3A). This may be attributed to differences in several methodological aspects including modality, source localization method, brain state segmentation technique, number of brain states, and the studied behavioral states. It is worth highlighting that Vidaurre et al. (2018) characterized the resting-state MEG (normal wakefulness with eyes open) with 12 brain states, while the current study used 10 states to investigate brain configurations during normal wakefulness (with eyes closed) as well as during subanesthetic and anesthetic doing of ketamine. In addition, the two states associated with sensorimotor cortex demonstrated activation in both alpha and beta frequency bands (States 3 and 4), which is in alignment with the spectral properties of the sensorimotor network identified in Vidaurre et al. (2018).

Ketamine had different effects on individual brain states (Figs. 2 D, 3 A–C). The posterior patterns (States 1 and 2) demonstrated decreases in the occupancy and the probabilities of persistence in and transitioning in/out of these states in a dose-dependent manner, which is in line with the finding that posterior alpha power decreases sequentially during subanesthetic and anesthetic dosing of ketamine (Li and Mashour, 2019; Vlisides et al., 2017). The anterior division of DMN (State 6) did not show significant changes in the occupancy and persistence probability, while transitions involving State 1 were suppressed with ketamine. These results suggest that ketamine disconnects the prefrontal cortex and more posterior regions, which might be associated with the disruption of connected consciousness while preserving vivid disconnected consciousness during ketamine anesthesia. The reductions were also seen in the patterns associated with sensorimotor cortex. While the beta activation state (State 3) showed gradual decreases from wakefulness to subanesthesia and then to anesthesia, the alpha activation state (State 4) only showed decreases during ketamine anesthesia. Prior studies have reported reduced beta activity within sensorimotor cortex using measures based on EEG power envelope (Forsyth et al., 2020), and the discrepancy between the two states might indicate the highly correlated but distinct roles of alpha- and beta-band rhythms in the modulation of the excitability of the sensorimotor cortex (Stolk et al., 2019). In contrast, ketamine increased the occupancy of the states associated with beta and gamma activation in visual cortex (States 5 and 7), which is in line with prior task-related findings that ketamine amplified beta and gamma activity in human visual cortex and that has previously been linked to cortical pyramidal cell inhibition in animal studies (Shaw et al., 2015). The strongest elevation was seen in the global patterns with the prevalence of delta, theta, and gamma activity across the cortex (State 8–10). Overall, we found that ketamine anesthesia tends to shift the configuration toward brain states with low spatial variability (Fig. S2C). This is consistent with the findings in a recent fMRI study, where the dominance of global activation and deactivation states has been speculated to be associated with psychoactive effects of ketamine (Huang et al., 2020).

A key finding of the current study is the modulation of the dynamic configurations of brain states by subanesthetic and anesthetic doses of ketamine (Figs. 3 E and 4 D–F). During subanesthetic ketamine, we observed elevated entropy in terms of state occupancy and state transitions and increased transition rate, indicating that (1) the cortical activity is more likely to visit all possible brain states, (2) at a certain time, networks are more likely to transition to a distinct state than persist in the current state, and (3) the state that the network transitions into is more unpredictable. These observations suggest that subanesthetic ketamine is associated with a more diverse set of brain configurations or a richer repertoire as compared to normal wakefulness, which is consistent with the entropic brain hypothesis that relates the richness of conscious experience with the entropy or diversity of a repertoire of brain configurations (Carhart-Harris, 2018). We also demonstrated an anesthetic dose of ketamine induced dynamic changes as evidenced by entropy that was reduced or comparable to that during normal wakefulness, as well as the increased rate of state transitions comparable to that during subanesthesia. These results, corroborated by our previous findings on local signal diversity (Li and Mashour, 2019), supports the hypothesis that ketamine-induced brain reconfiguration has features of general anesthesia (reduced entropy or diversity), normal consciousness (entropy or diversity comparable to normal wakefulness), and altered states of consciousness (elevated entropy or diversity), irrespective of the different temporal scales.

There were methodological limitations in this study that should be considered. First, the accuracy of EEG source localization is affected by using an anatomy template and default electrode positions in this study, which may be improved by integrating subject-specific MRI anatomy, exact electrode positions, and a more precise forward model that distinguishes between gray and white matter, etc. Alternatively, despite the discrepancies in some brain states, ketamine-induced changes in brain state dynamics are largely reproducible in scalp EEG signals when the volume conduction effect is appropriately controlled (Figs. S12 and S13). Second, we investigated the brain reconfiguration dynamics in terms of both the spectral and spatial properties and derived the brain states that represent frequency-specific spatial coactivation patterns. The spatial resolution or the separability of cortical regions into different states may be constrained as compared to the clustering solutions without differentiation of spectral information. Third, the dynamic brain states were inferred from the envelope data of the source-localized EEG signals. As the fluctuations of envelope and phase have been demonstrated to provide different information about the states of consciousness (Duclos et al., 2021), the brain states inferred from phase information or original signal may reveal distinct reconfiguration dynamics during ketamine anesthesia (e.g., the brain state time series exhibit Markovianity (Table S4) and non-periodicity (Fig. S14), different from the EEG microstate sequences derived from original signal (von Wegner et al., 2017)). However, envelope-based measures have been demonstrated, compared to phase-based measures, to be more consistent and reproducible across source-reconstruction pipelines (see Figs. S8 and S9 that demonstrate the robustness of our results across alternative reconstruction pipelines), subjects, and experiments (Colclough et al., 2016; Mahjoory et al., 2017). Fourth, we characterized the spectral properties of the brain states within the canonically defined frequency bands because of the frequency-dependent effect of ketamine on cortical activity. The EEG power spectrum consists of oscillatory components and a broadband, aperiodic background, which have been shown to associate with distinct resting-state networks (Jacob et al., 2021). We cannot exclude the possibility that the aperiodic 1/f-like activity could be a confounding factor in this study. Lastly, we employed k-means clustering approach to segment brain states, which assumes that all samples can be categorized into a set of mutually exclusive clusters and that the transition from one state into another is discrete. The same assumptions have been adopted in several brain state segmentation methods including EEG microstate analysis that segments EEG scalp maps into a series of microstates (Michel and Koenig, 2018). A recent study investigated these assumptions in the framework of microstate analysis and suggested that EEG microstates may be better conceptualized as spatially and temporally continuous, rather than discrete activations of neural populations (Mishra et al., 2020). It is worth noting that this finding is not in conflict with the existence of inferred brain states but rather suggests that incorporating the insight of a continuous trajectory may improve the performance of brain state segmentation, which is an area for future studies.

In summary, we characterized the dynamic organization of spontaneous electrophysiological activity with a set of frequency-specific spatial coactivation patterns and demonstrated that ketamine anesthesia tends to shift the configuration toward patterns with low spatial variability and that the effect of ketamine on the richness of the dynamical repertoire is dose-dependent. These results advance understanding of the neurophysiological mechanism of ketamine in terms of the spatial, temporal, and spectral structures of the underlying whole-brain dynamics.