Abstract
Understanding the effects of anesthesia on cortical neuronal spiking and information transfer could help illuminate the neuronal basis of the conscious state. Recent investigations suggest that the brain state identified by local field potential spectrum is not stationary but changes spontaneously at a fixed level of anesthetic concentration. How cortical unit activity changes with dynamically transitioning brain states under anesthesia is unclear. Extracellular unit activity was measured with 64-channel silicon microelectrode arrays in cortical layers 5/6 of the primary visual cortex of chronically instrumented, freely moving male rats (n = 7) during stepwise reduction of the anesthetic desflurane (6%, 4%, 2%, and 0%). Unsupervised machine learning applied to multiunit spike patterns revealed five distinct brain states. A novel desynchronized brain state with increased spike rate variability, sample entropy, and EMG activity occurred in 6% desflurane with 40.0% frequency. The other four brain states reflected graded levels of anesthesia. As anesthesia deepened the spike rate of neurons decreased regardless of their spike rate profile at baseline conscious state. Actively firing neurons with wide-spiking pattern showed increased bursting activity along with increased spike timing variability, unit-to-population correlation, and unit-to-unit transfer entropy, despite the overall decrease in transfer entropy. The narrow-spiking neurons showed similar changes but to a lesser degree. These results suggest that (1) anesthetic effect on spike rate is distinct from sleep, (2) synchronously fragmented spiking pattern is a signature of anesthetic-induced unconsciousness, and (3) the paradoxical, desynchronized brain state in deep anesthesia contends the generally presumed monotonic, dose-dependent anesthetic effect on the brain.
SIGNIFICANCE STATEMENT Recent studies suggest that spontaneous changes in brain state occur under anesthesia. However, the spiking behavior of cortical neurons associated with such state changes has not been investigated. We found that local brain states defined by multiunit activity had a nonunitary relationship with the current anesthetic level. A paradoxical brain state displaying asynchronous firing pattern and high EMG activity was found unexpectedly in deep anesthesia. In contrast, the synchronous fragmentation of neuronal spiking appeared to be a robust signature of the state of anesthesia. The findings challenge the assumption of monotonic, anesthetic dose-dependent behavior of cortical neuron populations. They enhance the interpretation of neuroscientific data obtained under anesthesia and the understanding of the neuronal basis of anesthetic-induced state of unconsciousness.
Introduction
Recent studies of large-scale brain activity found that multiple brain states appear at a constant anesthetic concentration; and conversely, one brain state can be observed at different anesthetic levels (Hudson et al., 2014; Hudson, 2017; D. Li et al., 2019). The degeneracy in the relationship between brain state and anesthetic concentration suggests that the neuronal network spontaneously shifts between two or more transient attractors indicating metastability (Hudson et al., 2014; Hudson, 2017; D. Li et al., 2019); that is, a system appeared to be in an equilibrium state at short timescales but switches to other states at longer timescales (Bovier, 2006).
Despite these observations, most studies of unit activity assume a one-to-one relationship between brain state and anesthetic concentration and investigate dose-dependent changes of neuronal activity (Vizuete et al., 2012, 2014; Sellers et al., 2013). This would be surprising and would suggest that the spiking dynamics of individual unit activities in different brain states under anesthesia has been poorly explored. Detailed information about the spiking dynamics during shifting brain states is arguably important for interpreting neuroscientific data obtained under anesthetized conditions and for understanding the neuronal mechanisms of altered states of consciousness.
In an attempt to fill this gap of knowledge, we measured single-unit spiking patterns of neuronal populations in chronically instrumented rodents subjected to multiple levels of anesthesia and applied machine learning to identify brain states independent of the actual anesthetic concentration. We hypothesized that brain states identified by specific features of population activity will show degeneracy in the relationship with anesthetic concentration and that these states will be characterized by distinct spike activity patterns.
Materials and Methods
Electrode implantation
The study was approved by the Institutional Animal Care and Use Committee in accordance with the Guide for the Care and Use of Laboratory Animals of the Governing Board of the National Research Council.
Eight adult male Long-Evans rats of 300-350 g weight were housed in a reverse light-dark cycle room for 5-7 d before surgical implantation. Ad libitum access to food and water was provided. A multielectrode array consisting of 8-shank, 64-contact silicon probes (shank length 2 mm, width 28–60 µm, thickness 15 µm, shank spacing 200 µm, row separation 100 µm, contact size 413 µm2; on-shank reference; custom design 8 × 8_edge_2mm100_200_413 (NeuroNexus Technologies) was chronically implanted in the deep layers 5/6 of the visual cortex of each rat. A craniotomy of rectangular shape of ∼2 × 4 mm was prepared over the right primary visual monocular area (center coordinates AP: 1.92 mm, ML: 3.0 mm, from λ), the exposed dura mater was resected, and the electrode array was inserted using a micromanipulator until the electrode tip reached the final position of 1.6 mm below the pial surface. The perimeter was covered with silicone gel (Kwik-Sil, World Precision Instruments). A sterilized stainless-steel screw in the cranium was used for ground, and additional screws were inserted for mechanical stability. The assembly was embedded with Cerebond (MyNeurolab). A pair of insulated wires (PlasticsOne), exposed at the tips, were positioned bilaterally into the nuchal muscles for EMG recording. Carprofren (5 mg/kg s.c., once daily) was administered during and 1 d after surgery. The animals were observed for 7-10 d for any infection or other complications.
Experimental design
One to eight days after surgery, the animals were placed in a closed, ventilated anesthesia chamber for continuous recording of extracellular potentials in a dark environment. Desflurane was administered at stepwise decreasing concentrations of 6%, 4%, 2%, and 0% and continuously monitored (POET IQ2 monitor; Criticare Systems). Core body temperature was maintained at 37°C by subfloor heating. After each step change, a 15 min equilibrium period was allowed for the anesthetic concentration to reach steady state (Fig. 1A). At each level, 20 min of neural recording was done in the resting state followed by visual stimulation. Visual inspection suggested that 20 min of recording at each anesthetic concentration was sufficient to reveal spontaneous fluctuations of brain state. In one experiment that was performed in the beginning of the study, only 10 min per anesthetic concentration was recorded. Because all measurements of neuronal activity (spike rate, burst ratio, etc.) were quantified from 10 s epochs, the shorter data length in this one experiment was considered inconsequential. Analysis of the neural response to stimulation was beyond the scope of this study and will be reported elsewhere.
Electrophysiological recording and identification of single units
Extracellular potentials were recorded using SmartBox (NeuroNexus Technologies) at 30 kHz sampling rate. The data were filtered for detecting unit activities (>300 Hz) and local field potentials (LFPs, <100 Hz). For artifact rejection, ±1 s of data around every time stamp with signal amplitude > 10 SDs was removed. The recordings were also visually inspected for noticeable noise episodes that were manually excluded from the analysis. One experiment was excluded from the analysis because of severe noise contamination (remaining experiments, n = 7). Single-unit activity (SUA) was identified using the clustering software Spiking Circus (Yger et al., 2018). Data were low-pass filtered at 7500 Hz and median-referenced. Units were then detected by thresholding and sorted by template matching. On average, 36 ± 14 (mean ± SD) single units were obtained per animal. The SUAs were further classified into wide-spiking (WS) and narrow-spiking (NS) units (Watson et al., 2016) based on the spike waveform, autocorrelogram, and cross-correlogram (CCG). Units with short half-amplitude width, short trough-to-peak time, and fast-spiking pattern (a prominent peak near 10-30 ms of autocorrelogram) were classified as a NS units; the rest of the units were classified into WS units (Fig. 1B,C). NS and WS presumably correspond to putative inhibitory and excitatory units, respectively (Csicsvari et al., 1998; Sirota et al., 2008). The common approach to identify putative excitatory and inhibitory monosynaptic connections is by calculating the CCG between neuron pairs (Fujisawa et al., 2008; Vizuete et al., 2012), but the chance of finding such connections among distant recording sites (100-200 µm) is relatively small. Therefore, we used an alternative approach by calculating the CCG between individual units (reference unit) and multiunit activity (MUA; the summation of SUAs) and then comparing the level of MUA before and after spike events of individual units (Fig. 1D). This was based on the conjecture that inhibitory (excitatory) units, on average, inhibit (excite) other units, resulting in a negative (positive) asymmetry in CCG. The asymmetry of CCG was defined as (X – Y)/(X + Y), where X (Y) is the number of spike events of all the other units 1-5 s after (before) the spike of the reference unit. All properties of LFP and spikes were calculated for nonoverlapping 10 s epochs by assuming stationarity over the timescale of anesthetic-induced slow oscillations and burst-suppression pattern.
Spectral analysis of LFP
For the spectral analysis, high-quality channels were manually chosen; the channels were chosen to have minimal cable movement, muscle, or cardiac artifacts. Furthermore, for a direct comparison with unit activity, only channels in which at least one SUA was detected were selected. On average, 11.0 ± 5.4 number of LFP channels (n = 7) were used for the analyses. For every time stamp with signal amplitude > 5 SDs, the periods ±1 s of those time stamps were removed; total 259.8 ± 184.0 s of LFP signals (5.5 ± 3.7%) were excluded. The fourth-order Butterworth filter (high-pass at 0.5 Hz) was applied to LFP signal. Power spectral density (PSD) of LFP in each epoch was obtained using the multitaper method in Python spectrum analysis software (pyspectrum.readthedocs.io) (Thomson, 1982). We used a time-half bandwidth product = 7.5, number of tapers = 14, yielding a frequency resolution of 1.5 Hz. The calculated PSDs from each epoch were concatenated to visualize the time-varying pattern of PSD (spectrogram). For the comparison of PSDs between different brain states, PSDs from epochs in each brain state were averaged.
To assess synchronization between individual spikes and LFP, spike-field coherence (SFC) was calculated. First, spike-triggered averages from ±2 s segments of LFP around each spike were obtained. Then, the SFC was estimated by the ratio of PSD of the spike-triggered average to the average of the PSD obtained from each individual LFP segment (Fries et al., 2001). Units with SR < 1 Hz were excluded from the SFC calculation.
EMG activity
The EMG signal was recorded with 1–500 Hz analog bandpass filter and 30 kHz sampling rate. This was used as a surrogate measure of the vigilance level. EMG signal was first downsampled to 3 kHz, and PSD was calculated using the same parameters with the PSD calculation of LFP signal. The PSD values with frequencies < 250 Hz were discarded because of cardiac artifact contamination. Next, overall EMG activity level in each consecutive epoch was estimated by the sum of the log-transformed PSD values. For a comparison across animals, EMG activity was transformed to a range between 0 and 1.
Single-unit spike properties
Spike rate (SR)
Because SR is known to follow lognormal distribution, linear-scale SR values were log-transformed (Buzsáki and Mizuseki, 2014) and averaged for nonoverlapping consecutive 10 s epochs. Zero SR values were substituted by SR = 10−2 Hz before the log transformation.
Gini coefficient
The Gini coefficient was used to estimate the dispersion of the SR distribution. It was originally intended to represent the income or wealth disparity, and is commonly being used in measurement of inequality. For non-negative values, the Gini coefficient can theoretically range from 0 to 1, 0 being complete equality and 1 being complete inequality. Gini coefficient was calculated from raw (not log-transformed) SR data, by plotting the neuronal population sorted by SR on the x axis and cumulative SR on the y axis (Lorenz curve) (Lorenz, 1905). The area below the Lorenz curve of the empirical SR data (area A) was then compared with the area below the Lorenz curve of an ideal SR data (area B), in which all neurons have an equal SR value. The Gini coefficient value was finally defined as a ratio, (B – A)/B.
Burst ratio (BR)
Two spikes with short interspike intervals (ISIs) (<10 ms) were considered as an indication of bursting spike. BR was defined as the number of ISIs < 10 ms (bursting spikes) divided by the total number of ISIs in each epoch. Units with SR < 1 Hz in each epoch were considered as inactive units and excluded from the BR calculation. BR values were log-transformed and averaged for consecutive 10 s epochs.
Local variation (LV)
Spike timing variability, or spike irregularity, was estimated by LV (Shinomoto et al., 2003) from each spike train of SUA. LV was defined as follows:
where 
is the duration of ith ISI and 
is the number of ISIs. LV is 0 for constant 
and approaches 1 for a sufficiently long Poisson ISI sequence. LV is thought to differentiate the degree of intrinsic spiking randomness of individual neurons more effectively than the other measures, such as coefficient of variation of ISI (Shinomoto et al., 2003). LV values were not log-transformed because they did not show lognormal distribution.
Transfer entropy (TE)
TE was used to estimate directional functional connectivity among individual units (Schreiber, 2000; Ito et al., 2011). For two spike trains of units x and y, TE can be estimated as follows:
where m denotes embedding dimension (pattern size), and p(·) implies probability. 
denotes m size spike pattern. For example, for m = 3 cases, there are 8 (=23) possible spike patterns ([0,0,0], [0,0,1], …, [1,1,1]). 
(
) measures the statistical influence of unit y (10) on unit x (y). 
is the reduced amount of uncertainty in the future of x by knowing the past of y given the past of x. TE can also be thought of as mutual information (I) between 
and 
given the past of 
, as follows:
We used m = 3, and all spike trains were downsampled to 125 Hz before calculation of TE; that is, the individual value in each bin of the spike train is 1 if there is one or more spikes within the 8 ms bin and 0 otherwise.
Multiunit spike properties
Three measures, total number of spikes (TNS), longest period below mean (LPBM), and sample entropy (SpEn), were calculated with MUA signal (i.e., sum of all SUAs). TNS represents the amount of total spike events that occur in a sampled neural network. For the calculation of LPBM and SpEn, the MUA signal was convolved with Gaussian kernel with SD of 25 ms (Vyazovskiy et al., 2011), and continuous spike signal was obtained. Details of LPBM and SpEn estimation are described below.
LPBM
LPBM measures the time length of the longest inactive periods (or active periods depending on the time series characteristics) in a given time series (Lubba et al., 2019). First, the time lengths of all consecutive values below the mean of time series are calculated; then the maximum of the time lengths is obtained as an LPBM value. LPBM is known to be one of the important temporal statistics in time series analysis (Lubba et al., 2019), and was used in this study to be a measure of persistent inactiveness of spike activity. A high LPBM value in MUA signal implies an existence of a long inactive period suggesting synchronous fragmentation of spike activities, whereas a low LPBM indicates more continuous activity suggesting an irregular and asynchronous spiking pattern. Therefore, LPBM is expected to yield a high value when spikes are synchronously fragmented in time (e.g., slow oscillation and burst suppression). In addition, LPBM will further increase as burst-suppression ratio increases with deepening of anesthesia.
SpEn
SpEn was used to estimate the statistical irregularity of MUA as a time series. SpEn is an approximation of Kolmogorov entropy that measures the predictability of consecutive time series values based on their past values (Richman and Moorman, 2000). A high SpEn value implies random or unpredictable dynamics, whereas a low SpEn value indicates regular or deterministic dynamics. SpEn has been used to quantify depth of anesthesia and the level of consciousness in electroencephalogram (EEG) studies (Liang et al., 2015; Liu et al., 2018), and it generally decreases as anesthetic deepens. To calculate the SpEn, first an embedded time series is obtained as follows:
where 
is the time series value (convolved MUA signal in this study) at time t, and m is the embedding dimension (pattern size). Second, the correlation sum is calculated from the embedded time series as follows:
where 
denotes a Heaviside step function, 
implies Euclidean distance between two vectors, and r represents the distance criteria. We used m = 3, and r = 0.2 SD of amplitudes within each epoch following previous literature (Liang et al., 2015; Liu et al., 2018). Finally, the SpEn is defined as follows:
Before the SpEn calculation, the MUA signal was convolved as in the case of LPBM calculation and downsampled to 125 Hz.
Classification of brain states
The primary focus of the study was to examine brain state-dependent changes in spike activity patterns during and after anesthesia. To this end, spike train data were first segmented into 10 s nonoverlapping epochs. Then five features of spike activity of neuronal population were measured in each epoch. The five features were chosen with the intent to account for well-known effects of anesthesia on cerebral cortex.
First, with an exception of some dissociative anesthetic agents (e.g., ketamine), anesthesia generally reduces spike activity of neurons. This was measured by mean SR (SRm) of individual units and TNS of the neuron population. Because the effect of anesthetic on the SR of individual units was variable, SRm and TNS showed distinctively different traces (Fig. 2A); thus, we included both quantities. Second, spikes of individual neurons are temporally fragmented (Lewis et al., 2012; Vizuete et al., 2014). LV, a measure of spike timing variability, was suited to quantify this change (Shinomoto et al., 2003). Third, population activity (EEG, LFP, and MUA) slows down and an oscillatory pattern emerges. These changes result in an increase of predictability and suppressed complexity, and can be well captured by SpEn (Richman and Moorman, 2000). Last, deep anesthesia is characterized by burst suppression, an alternation of brief EEG activity and electrocortical silence (Clark and Rosner, 1973). The duration of silent periods monotonically increases during the transition from slow oscillation to isoelectric state. To quantify this pattern in terms of population spiking activity, we applied the LPBM method (Lubba et al., 2019) to MUA. SR and LV were calculated from SUA, then averaged over the neuronal population. TNS, SpEn, and LPBM were measured from the MUA signal.
Although the five features reflect the key effect of anesthesia on neural activity, it is possible that one would obtain different clustering results with a choice of different features. Moreover, there could be other anesthetic effects that were omitted from our feature selection. With the use of additional clustering features, the brain states defined could be further split or combined. In a preliminary study (unpublished data), we performed another clustering analysis with adding either permutation entropy (Bandt and Pompe, 2002) or coefficient of variation of MUA signal. These analyses revealed five clusters similar to those presented here.
The unnormalized values of the five features have different ranges with each other, and the single feature values have different ranges across the animals. Therefore, to normalize the feature values and to mitigate the effect of outliers, the feature data were divided into sextiles for each animal, and transformed by linearly scaling to a given range (0-1), so that the median of the data in the first sextile was considered 0 and the median of data in the last sextile was considered 1 in each experiment. This procedure is based on an assumption that the range of each of the five features is similar across different experiments with the same anesthetic protocol.
The five features were then used for unsupervised clustering to delineate distinct brain states. Hierarchical agglomerative algorithm with Ward's linkage method were applied for the clustering of brain states, using Python package Scikit-Learn (www.scikit-learn.org). Each data point of a 10 s epoch was first treated as a single cluster in feature space; then the points were successively merged until all clusters merged into a single cluster. The method does not require a specific number of clusters (K) at the beginning step, and the clusters can be easily identified from the hierarchy tree (dendrogram) that is built from the algorithm. We determined the optimal number of clusters based on the dendrogram and so-called elbow method. A within-cluster distance was plotted against the number of clusters, and the point where the curve sharply bends was chosen as an “elbow” point. We used the maximum of the second-order difference of the distance-K curve to find the elbow point. We neglected the K = 2 cases, in which the brain state simply represents anesthetized (6%, 4%, and 2% desflurane) and waking state (0% desflurane).
Prior studies attempted to define brain states based on LFP or EEG signals (Chander et al., 2014; Hudson et al., 2014; Watson et al., 2016; Hudson, 2017; D. Li et al., 2019). We additionally applied LFP spectrogram-based clustering and compared it with the spike-based clustering result. The LFP spectrogram was log-transformed, and principal component analysis was performed. As in spike-based clustering, the agglomerative algorithm was applied to the first five principal components, which contained 70.2 ± 4.1% of total variance (n = 7) similar to previous findings (Hudson et al., 2014). We chose K = 5 to be consistent with spike-based clustering. Because of the interindividual variability of LFP spectrograms that might have originated from inadvertent differences in electrode placement, impedance, or anesthetic sensitivity, the clustering analysis was performed in each animal separately. For the comparison of spike-based and LEP-based clustering analyses, the resultant clusters of the two methods were matched using Hungarian algorithm (Munkres, 1957).
Statistical analysis
All statistical analyses were conducted using StatsModels library (www.statsmodels.org) in Python 3.7. For all measures, to test the difference across brain states, statistical comparisons were first performed using linear mixed models based on restricted maximum likelihood estimation. For all linear mixed models, the brain states (categorical independent variable) were used as a fixed effect. For the properties of population activity (i.e., PSD of LFP, the five input features, and EMG), the random effect included all 7 animals. For the individual unit properties, such as SR and LV, the random effect included the animals and units. Statistical significance was tested with F tests. For post hoc pairwise comparisons between brain states, the Bonferroni-corrected p value 0.05 (number of hypotheses = 10) was used, unless stated otherwise.
Results
CCG between SUA and MUA
WS and NS units were classified based on their spike waveform and autocorrelogram; 36 of 251 units were classified as NS units (14.4%). In order to test whether NS (WS) corresponds to putative inhibitory (excitatory) neurons, we examined the asymmetry of CCG between SUA and MUA. As predicted, the asymmetry of NS (WS) units showed negative (positive) values on average (Fig. 1E). Statistical significance was found both in WS and NS units (t test without Bonferroni correction after linear mixed model, t(214) = 2.0, p = 0.045 for WS, and t(35) = −5.17, p < 10−6 for NS). This suggests that NS (WS) units tend to inhibit (promote) population activity, reassuring that they might be classified as putative inhibitory (excitatory) neurons.
Anesthetic procedure and identification of WS and NS units. A, Schematic representation of the anesthesia protocol. The volatile anesthetic desflurane was administered at steady-state concentrations and decreased in a stepwise manner. This study used spontaneous activity data only. B, WS and NS units were classified based on spike waveform features and the autocorrelogram (data not shown). C, Average spike waveform of WS (n = 215) and NS (n = 36) units. D, Average CCG between SUA and MUA. The CCG was normalized by dividing the number of units in MUA. E, Asymmetry of CCG confirms the identification of WS and NS units. Error bar indicates 95% CI across units. *p = 0.045. ***p < 10−6.
Brain state shifts during anesthesia
In order to identify local brain states from the electrophysiological recording independent of the nominal anesthetic concentration, we visualized how LFP spectrogram and MUA characteristics changed over time during the experiment. In all animals, the LFP spectrogram, TNS, SRm, mean LV (LVm), LPBM, and SpEn profoundly changed during (6%, 4%, and 2% desflurane) and after (0% desflurane) anesthesia. Figure 2A illustrates the time course of these variables in 1 animal, for example. Importantly, both the spectrogram and the MUA features changed not only between but also within each recording period at constant anesthetic concentration. For instance, in the middle of the recording at 6% desflurane, low-frequency (<4 Hz) power in the LFP spectrogram and LPBM abruptly increased for no evident reason. Additional abrupt transitions were seen at 2% desflurane. Other animals also showed similar transitions (Fig. 3E) at various anesthetic levels. This example demonstrates that a simple one-to-one relationship between the chosen LFP/MUA variables and the anesthetic concentration does not exist, suggesting the need for a more nuanced identification of brain states from these variables.
Classification of brain states. A, LFP spectrogram and traces of five features from a representative animal. Colors in the horizontal bar above the spectrogram represent different brain states (S1-S5) as classified by clustering. B, Scatter and histogram plots of the five features used for state clustering from pooled data. C, Dendrogram generated from hierarchical agglomerative clustering. D, Elbow method suggests the optimal number of clusters as K = 5. E, Brain state-dependent changes of the five features. Error bar indicates 95% CI across animals. Inset (in each panel), Statistically significant difference between pairs of brain states. F, Relative frequency of four anesthetic concentrations supporting each brain state, pooled data from 7 animals. G, Relative frequency of five brain states at each anesthetic concentration, pooled data from 7 animals.
LFP and EMG in five brain states. A, Example LFP traces from 1 animal. B, PSD versus frequency plot shows power-law relationship in the five brain states. Data were averaged from the 7 animals. C, Relative PSD versus frequency plot in a linear scale. D, L/H power ratio across five brain states. The L/H ratio is defined as log10{(PSD at 0.25-4 Hz)/(PSD at 30-59 Hz)}. Error bar indicates 95% CI across animals. Inset (in each panel), Statistically significant difference between pairs of brain states. E, EMG activity (black trace) during 6% to 4% desflurane in 4 animals. Horizontal bars with different colors represent different brain states: magenta, yellow, cyan, and green for S1, S2, S3, and S4, respectively. F, EMG in S3 is significantly higher than that of S1 and lower than that of S5. Error bar indicates 95% CI across animals. *p = 0.023; ***p < 10−6; compared with S3 (Bonferroni correction). G, H, Brain states were obtained using LFP spectrogram-based clustering. Two representative animals are shown. Left, Horizonal color bars represent brain states defined by spike-based and LFP-based clustering methods. The traces of two representative quantities, TNS and the first principal component (PC1), are shown. Right, Contingency table is represented as a 2d histogram of spike-based and LFP-based clustering results.
To achieve this goal, we used agglomerative clustering on five MUA variables as input features from data pooled from all animals to identify distinct, unitary brain states. The scatter plots in Figure 2B illustrate pairwise relationships of the five MUA features in 5 clusters. The choice of 5 clusters was based on the dendrogram (Fig. 2C), which illustrates that between-cluster distances were large and within-cluster distances were small at K = 5. We also calculated within-cluster distance as a function of K (Fig. 2D). The second-order difference of the distance curve was maximized at K = 5, suggesting that it was an optimal choice consistent with the dendrogram distances.
The five clusters identified by unsupervised clustering were designated as brain states S1 to S5, and the mean values of MUA variables among these states were statistically compared (Fig. 2E; Table 1). S1 was characterized by the lowest TNS, SRm, and SpEn and the highest LPBM, indicating that S1 corresponded to burst suppression (Fig. 3A), which is typical to deep anesthesia. Indeed, S1 was mostly observed at 6% desflurane (Fig. 2F,G). S5, on the other hand, was mostly observed at 0% desflurane (Fig. 2F). This was characterized by high spike activity (high TNS and SRm) and asynchronous firing patterns (high SpEn and low LVm). S2 and S4 had intermediate feature values between those of S1 and S5. S2 was mostly observed at 4% desflurane and S4 was mostly seen at 2% desflurane (Fig. 2F).
The five features (TNS, SRm, LVm, LPBM, and SpEn) in the five brain statesa
Interestingly, S3 was mostly found at 6% desflurane (Fig. 2F) similar to S1. However, S3 showed a distinct pattern from S1. It was characterized by high SpEn, relatively low SRm, and very low LBPM. TNS was not reduced as much as SRm (TNS indicates total number of spikes in the neuronal population, and SRm is the mean of log-transformed individual SRs). The discrepancy suggests that, in S3, some neurons were inactive, but a few neurons emitted a large number of spikes, whereas these outliers were absent in S1. The details of SR distribution of neuronal population are analyzed further in the next section. The low LPBM indicates no clear distinction of active and inactive period in S3. The low LPBM, together with the high SpEn, implies that spiking pattern in the S3 was asynchronous.
Although brain state and anesthetic concentration were not uniquely related, a general trend of the occurrence probability of brain states with anesthetic level was evident; the S1, S2, S4, and S5 in order occurred mostly at correspondingly decreasing desflurane concentration. Therefore, it is reasonable to surmise that the occurrence probability of brain states, with an exception of S3, generally reflected the depth of anesthesia. Interestingly, however, they were also observed in other concentrations (Fig. 2F,G). For instance, at 6% desflurane, S1 was present 61% of the time, whereas at 4% desflurane, S1 was present 39% of the time, with the balance occupied by other brain states. It was seen that several different brain states occurred at each constant anesthetic concentration. For example, at 6% desflurane, S1, S2, S3, and S4 appeared in non-negligible proportion (Fig. 2G). The many-to-many relationship between brain state and anesthetic concentration suggests a general need for brain state-dependent investigation of unit activity.
LFP properties of the five brain states
Because LFP generally reflects the state in anesthesia, we examined LFP patterns and PSD in each brain state. Typical LFP traces in five brain states are shown in Figure 3A from the same animal as in Figure 2A. The LFP in S1 exhibited burst suppression. S2 and S4 revealed relatively high amplitude, slow activity as generally expected in anesthesia. In contrast, S3 showed relatively low amplitude and desynchronized LFP pattern. The PSD averaged over all animals showed a power law relationship with frequency (Fig. 3B). The slightly higher slope in S2 and S4 was associated with the increased low frequency (<4 Hz) and decreased high frequency (>30 Hz) PSD (Fig. 3C) consistent with the anesthetic-induced suppression of high-frequency γ power and enhancement of δ and slow oscillation in EEG/LFP. S5 was characterized by increased theta (5-9 Hz) and high-frequency (>20 Hz) power, the typical signatures of EEG/LFP in wakefulness. For a quantitative comparison, we calculated the L/H ratio as log10{(PSD at 0.25-4 Hz)/(PSD at 30–59 Hz)} (C. Y. Li et al., 2009). S2 and S4 showed significantly higher L/H ratio than S5 (F test: F(4,30) = 3.94, p = 0.011; post hoc test: t = −3.46, p = 0.005 and t = −3.12, p = 0.018 for S2 and S4, respectively, with Bonferroni correction; Fig. 3D). In sum, the brain states S1, S2, and S4 were consistent with known LFP features of deep, moderate, and light anesthesia, respectively. However, the LFP in S3 was unexpected and contrary to the generally presumed dose-dependent effect of anesthesia.
High EMG activity in paradoxical desynchronized state
The asynchronous firing pattern and relatively high LFP γ power found in S3 raise the question whether the systemic arousal level may also be elevated in S3 as it is in S5. Generally, the EMG follows the level of arousal; therefore, the vigilance state of animals was estimated by EMG activity. Although both S1 and S3 occurred mostly in 6% desflurane, EMG of S3 was substantially higher than that of S1. The rescaled EMG traces from each animal exhibited higher muscle activity in S3 than in S1 and sometimes even higher than in S2 (Fig. 3E). Statistically significant differences in the rescaled EMG were found for S3 versus S1 and S3 versus S5 (F test: F(4,30) = 21.8, p < 10−7; t test of S3 vs others: t = −2.76, p < 0.023 and t = 5.48, p < 10−6 for S1 and S5, respectively, with Bonferroni correction [number of hypotheses = 4]; Fig. 3F).
Comparison of spike-based and LFP-based clustering
The LFP spectrograms (Fig. 2A) and PSD (Fig. 3B–D) suggest that brain states might be defined using the LFP or EEG signal as they were in prior studies (Chander et al., 2014; Hudson et al., 2014; Hudson, 2017). Therefore, we also applied LFP spectrogram-based clustering and compared it with the spike-based clustering result. First, consistent with the findings from spike-based clustering (Fig. 2), there was many-to-many correspondence between LFP-based brain state and anesthetic concentration (Fig. 3G,H, left panels). One state was seen at multiple desflurane concentrations, and multiple brain states were found at the same concentration. Second, LFP-based brain states and spike-based brain states displayed qualitatively similar changes during the whole experiment (Fig. 3G,H, left panels). The contingency matrix (Fig. 3G,H, right panels) from 2 representative animals illustrates the similarity of results from the two clustering methods. The correlations between were higher than the expected frequencies, based on the marginal sums of the contingency tables; χ2 test: χ2(16) = 346.0, 623.8, 420.6, 748.6, 635.9, 379.4, and 230.5; p < 10−63, 10−121, 10−78, 10−148, 10−124, 10−70, and 10−40 for each animal. The match ratio (trace of the contingency matrix divided by the number of epochs) was 0.484 ± 0.105 (n = 7). However, there were also occasional discrepancies between LFP-based brain states and spike-based brain states. For example, S2 and S4 of the two methods were often interchanged (Fig. 3H, right panels). The findings from LFP-based clustering confirmed the degeneracy of the relationship between brain state and anesthetic concentration.
Spike rate distribution across brain states and neuron types
We compared the five brain states in terms of both the number of emitted spikes and the average SR of individual units. Generally, desflurane suppressed spike activity (Fig. 4A,B). Figure 4A illustrates the time course of SRm (average of log-transformed spike rate) and total spike number TNS from the same animal as in Figures 2A and 3A. The traces of SRm and TNS deviated from each other, especially at 6% desflurane. The TNS showed a pronounced decrease when the brain state transitioned from S3 to S1, whereas the SRm remained the same. Figure 4B depicts SR distribution of the neuronal population in the five brain states. When brain state changes from S5 through S4 and S2 to S1, there was a progressive, overall shift in the SR distribution, indicating that SR of all units, whether they are high-firing or low-firing, was reduced by desflurane. In S3, many units were inactive, even more than in S1 and S2, but there were a few units with very high SR (Fig. 4B). Accordingly, the variation of SR across individual units was the highest in S3. The variation in SR distribution was quantified by the Gini coefficient, and the value of S3 was significantly larger than all others (Fig. 4C; Table 2). Specifically, Figure 4D implies that the highly active units in S3 are WS units. The SR of active units in S3 was comparable to that in S5 (Fig. 4E, left), for WS units; however, the SR of NS active units in S3 was lower than that of NS units in S5 (Fig. 4D, right). For the mean SR of WS units, there were significant differences among the states, except S1 versus S3; SRm increased from S1 or S3 through S2 and S4 to S5 (Fig. 4D; Table 2). For NS units, SRm was significantly higher in S5 than in all other states. Thus, in general, desflurane profoundly suppressed SR of both WS and NS units, but a few WS units in S3 remained highly active, resulting in very high SR variation in this brain state.
SUA features in the five brain statesa
SR properties of five brain states. A, Example trace of SRm (black) and TNS (blue) from 1 animal. Black shaded area represents 25th and 75th percentiles of SR distribution. Horizontal bars with different colors represent different brain states. B, Frequency distribution of log-transformed SR in five states. C, Gini coefficient of untransformed SR distribution. Error bar indicates 95% CI across units. Inset, Pairwise statistical significance. D, Average SR as a function of SR threshold for WS (left) and NS (right) units. E, Comparison of SRm across five states. Inset, Statistically significant difference for pairs of brain state. F, Pairwise SR correlation among all epochs from the same animal in A. Two color bars above and right of the matrix represent the five brain states (same color code as before). G, Orthogonal linear regression of SRs between brain states. For the intrastate comparisons (diagonal elements), the SR of the first half of the data were compared with the SR of the second half of the data. Top right, Inset, Pearson correlation values for within- and between-brain states.
The changes in SR distribution (Fig. 4B) suggest that the decrease in average firing rate could be generalized across all units. Therefore, we further tested whether units had a tendency to preserve their firing rate rank across brain states. SR similarity between any two epochs was estimated by calculating the Pearson correlation, and is presented in Figure 4F. The correlation matrix for individual units showed to be high within-state similarity and relatively low between-state similarity. The results from correlation analysis of all units from all animals (n = 251) were consistent with the result from the representative animal. Within-state comparisons of SR between the first half of the state and the second half of the same state are shown in the diagonal panels of Figure 4G. Between-state comparisons between different states are shown in the off-diagonal panels of Figure 4G. Orthogonal linear regression indicated that within-state similarity of SR (R > 0.93 for all the five states) is generally higher although still significant (p < 0.001) than between-state similarity, except S1 versus S2 (Fig. 4G, top right inset). These findings indicate that SR profiles of individual units are preserved both within and between brain states.
Temporal dynamics of spike activity
Anesthesia not only suppresses the average SR as reported in the previous section (Fig. 4) but also modulates the temporal dynamics of spike activity (Vizuete et al., 2014). Raster plots in Figure 5A illustrate the changing temporal dynamics of spike activity. In S1, S2, and S4, but not in S3, spike activity is more synchronized and temporally fragmented compared with S5. Figure 5B displays raster and distribution of ISIs in seven representative units (four WS and three NS) from the same animal at different desflurane concentrations. The shape of ISI distribution was profoundly altered by the anesthetic. In wakefulness (S5), the ISI distribution was unimodal, whereas in the other four states it was bimodal or multimodal. This was partially because of the silent periods in spike activity. The large ISI values in the ISI raster plot (ISI ∼ 103 ms) in Figure 5B (especially, S1 and S2) correspond to silent periods that contribute to a second (upper) peak in ISI histogram (Fig. 5B). In addition, some WS units in S1 and S2 tended to fire in brief bursts that were associated with short ISI (ISI < 101 ms; Fig. 5B). Burst activity had also contributed to a peak near ISI ∼ 10 ms in the ISI histograms (Fig. 5B). In S3, two WS units exhibited very high SR (represented by dense points, asterisk in Fig. 5B) that was comparable to the SR in S5, consistent with the findings in the previous section (Fig. 4B–D). Both units showed unimodal ISI distribution as in S5.
Temporal dynamics of individual units. A, Example LFP trace and corresponding raster plot of single units in 1 animal in brain states S1-S5. B, Scatter plot and frequency distribution of ISI of seven units from the same animal in A. Red histograms represent WS units. Blue histograms represent NS units. Two WS units in S3 (*) displays extraordinarily high SR with unimodal ISI distribution. C, Average normalized autocorrelogram on log-scale calculated from all active units (SR ≥ 1 Hz). Inset, The same autocorrelogram on linear scale for short time lags (−50 to 50 ms). D, Distribution of BR of active units (SR > 1 Hz). E, Comparison of BR across five states. Error bar indicates 95% CI across units. Inset, Statistically significant pairs of brain state. F, An illustration of the change of BR and in(activeness) with brain state in individual WS units shows that WS units become either bursty or inactive under anesthesia; each line connects data from the same unit. A dark (light) colored line indicates low (high) BR value in S5. Inactive units (SR < 1 Hz) are shown at the bottom of each panel. For better visualization of the gradual change from S1 through S2 and S4 to S5, data in S3 were separated to the left. Right, Emphasis of two extreme cases (S1 and S5). G, Distribution of LV of active units (SR > 1 Hz) only. H, Comparison of LV across five states. Error bar indicates 95% CI across units. Inset, Statistically significant difference for pairs of brain state.
To determine whether burst activity and long silence periods generally occurred in all units and all animals, an autocorrelogram was calculated from all active units (SR ≥1 Hz). The averaged autocorrelogram showed a gradual increase of burst activity (ISI < 10 ms) in WS units under anesthesia from S5 to S1 (Fig. 5C, left). A prominent peak was observed near 6 ms that progressively decreased from S1 to S5 (Fig. 5C, left). As expected, the autocorrelogram showed little or no evidence of bursting of NS units. Another measure of bursting of WS units, the BR, generally decreased from S1 to S5 (Fig. 5D,E). Inactive units (SR < 1 Hz) were excluded from the autocorrelogram and BR calculation. Similar to the SR distribution, BR did not follow normal distribution but skewed to the right; thus, it was log-transformed. A statistically significant difference in BR of WS units was found for all pairwise comparisons of brain states, except S4 versus S2, and S3 (Fig. 5E; Table 2). The increases in BR of NS units were less pronounced (Fig. 5C; right; Fig. 5E). For NS, BR of S3 was significantly lower than that of S1 and S5. The suppression in SR (Fig. 4), together with the changing temporal pattern of spiking, indicates that neurons were either inactive or bursty at deeper levels of anesthesia. As the brain state changed from S1 to S5, more units became active and burst activity of active units in anesthesia was reduced in S5 (Fig. 5F). Again, S3 was an exception; BR of WS units in S3 was comparable to BR in S4.
The state-dependent changes of ISI distribution were also characterized by LV, a measure of spike timing variability, a measure that is sensitive to changes in both burst activity (small ISI) and to the presence of long silent periods (large ISI). LV showed a similar trend to BR across brain states but more effectively distinguished the five brain states, especially for NS units (Fig. 5E, inset; Fig. 5H). LV of both WS and NS units decreased from S1 through S2 and S4 to S5. All pairwise comparisons, except S1 versus S2, were statistically significant for WS units; and all pairwise comparisons, except S1 versus S2, and S3, were significant for NS units (Fig. 5H; Table 2).
In summary, spiking activity of most units was profoundly inhibited by desflurane, and the remaining active units showed an enhanced burst activity (for WS) and prolonged silence period (both WS and NS). In the paradoxical state S3, these effects were marginal such that units exhibited an irregular spiking pattern similar to that seen in wakefulness or S5.
Individual neurons conform to population activity
It is well known that desflurane as well as other anesthetics enhance correlation between spiking activity and LFP (Vizuete et al., 2014). We reexamined this effect as a function of brain state and found that desflurane increases spike-LFP coherence in all anesthetized states (Fig. 6A). Specifically, the SFC at low frequency (<5 Hz) was high in S1, S2, and S4. Although S3 was mostly found in 6% desflurane, the SFC of S3 was lower than that of S4. We also examined spike-triggered MUA (Fig. 6C). The MUA in S5 was high and relatively flat across time lags, with a small oscillatory pattern in theta frequency range (5–9 Hz). From S4 through S2 to S1, the overall MUA level gradually decreased, indicating a suppression of overall spike activity, while the MUA peak near 0 ms remained almost the same, indicating synchronous firing. In addition, the “dip” in MUA observed before and after spike events became deeper and wider as brain state moved from S4 through S2 to S1. In S1, the number of spikes near ±500 ms to spike events is close to zero, consistent with the near-silent periods of LFP burst suppression. Again, distinct from the other three anesthetic states, S3 did not have a large trough on either side of spike events, whereas the MUA was substantially lower than in S5.
Functional connectivity in five brain states. A, SFC as a function of frequency. B, Mean SFC of slow-δ oscillation (<4 Hz). Error bar indicates 95% CI across units. C, Spike-triggered average MUA. D, Pearson correlation between SUA and MUA. Error bar indicates 95% CI across units. Inset, Statistically significant difference for pairs of brain state. E, All pairwise TE of individual units from all animals on the axis of SR of individual units. Color of each dot in each panel represents TE value between two units: red represents high TE; blue represents low TE. The x axis (y axis) corresponds to SR of information-receiving (-sending) unit. F, Comparison of TE of active units (SR ≥ 1 Hz) across five brain states. Error bar indicates 95% CI across animals. Inset, Statistically significant difference between brain states. G, TE as a function of SR threshold; that is, TE values of active units are averaged while criteria for active and inactive units vary. States correspond to those as in the panel above. H, Average log-transformed SR of WS and NS units as a function of SR threshold. I, Lorenz curve of pairwise TE (including inactive units) from pooled data. J, Gini coefficient of the TE distribution. Error bar indicates 95% CI across animals. Inset, Statistically significant pairs of brain state. K, Information Venn diagram to illustrate the relationship between two neurons (X and Y) in wakefulness and anesthesia. Yellow area represents the amount of information transfer from Y to X (
). High spiking activity in wakefulness suggests high TE (left). In anesthesia, neurons are either inactive or intermittent/burst firing with enhanced synchronization. TE between inactive neurons is extremely small, whereas TE between active neurons is larger than the TE in wakefulness. This results in a large variation of TE values as seen in I, J.
To evaluate the extent to which individual neurons conform to population activity, we quantified SUA-MUA correlation. Both spike train of SUA and MUA signal were convolved with Gaussian kernel (SD = 25 ms); then the Pearson correlation between the two convolved signals was calculated. A substantial change in correlation was observed in both WS and NS units (Fig. 6D); all pairwise comparisons were statistically significant for WS units, and all pairwise comparisons except S3 versus S5 were significant for NS units (Table 2).
Information transfer depends on spike rate
The correlation results described so far indicate a nondirectional relationship between convolved SUA and MUA signals. In order to estimate directional functional connectivity of neuronal interaction, TE between individual binary spike trains (SUAs) was calculated. Because spike activity itself is in part a result of neural communication, TE is presumed to depend on the degree of overall spike activity. Indeed, TE was high for units with high SR and low for units with low SR (Fig. 6E). However, for a same range of SR, TE in anesthesia was higher than that in wakefulness. For example, for units having SR in a range of 100 to 101 Hz (the colored points inside the black squared box in each panel in Fig. 6E), TE values in S1 were higher than TE values in S5.
Synchronous firing correlates with enhanced information transfer
The sum of TE values (TEs) for active units was highest in S5 and lowest in S1 and S3 (Fig. 6F, left; Table 2), which is consistent with the general reduction of SR in anesthesia. However, the mean TE (TEm; the mean of TE values for active units) showed an opposite trend as it was the highest in S1 (Fig. 6F, second panel from the left; Table 2). Inactive units (SR < 1 Hz) were excluded from all results in Figure 6F, and the averaging was done across animals (n = 7). To see whether the SR exclusion threshold had any effect, TE results with different SR thresholds were compared in Figure 6G. When all units were included (SR threshold = 0 Hz), the TEm was the highest in S5 (Fig. 6G, second panel from the left) as it should be qualitatively the same as in TEs with SR threshold = 0 Hz (Fig. 6G, first panel from the left). However, as the SR threshold increased, TEm of S5 became the lowest and TEm of S1 the highest (Fig. 6G, second panel from the left). Figure 6H illustrates that the differences among brain states do not merely reflect the SR changes. The mean SR in S5 was the highest of all states across all SR thresholds, distinct from the TEm for both WS and NS units. For WS units, SR of S3 was comparable to that of S5 when SR threshold > 0 Hz.
The variation of all pairwise TE values (including inactive neurons) estimated by Gini coefficient was higher in anesthesia (S1-S4) than wakefulness (S5) (Fig. 6I). This is because, in anesthesia, many of the units were silent or inactive; thereby, these units had very low TE, while the remaining active units had relatively high TE. Statistically significant differences in the Gini coefficient were found for S5 versus S1, S2, S3, and S4 and S4 versus S1, and S3 (Fig. 6J; Table 2).
The reason for the difference in change across brain states between the TEs and TEm can be attributed to (1) the number of active neurons and (2) synchronous activity (Fig. 6C,D), as explained by Venn diagram of information in Figure 6K. In wakefulness, there are many active neurons; therefore, the sum of transfer entropies of active neurons (TEs) is high (left in Fig. 6K). In anesthesia, on the other hand, there are far less number of active neurons (Fig. 4B); therefore, the sum of transfer entropies of active neurons (TEs) is relatively small. For the small number of active neurons (Fig. 6K, top right), the enhanced synchronization in anesthesia produces a large value of TE; this contributes to the high value of the TEm between active neurons. For inactive neurons, however, there is few information to be transferred, so that TE is extremely small (Fig. 6K, bottom right). This also explains why the variation in TE is high in anesthesia (Fig. 6I,J).
Information transfer along different connection types is state-dependent
We further investigated whether desflurane differentially affects different connection types, by examining TE of neurons pairs, WS-to-WS, NS-to-WS, WS-to-NS, and NS-to-NS (Fig. 6F,G). The most pronounced change with state was observed in WS-to-WS connectivity indicated by a gradual decrease of TE from S1 to S5 (Fig. 6F). Statistical significance was seen in S1 versus all the other states and S2 versus S5 (Table 2). The NS-to-WS connectivity was also higher in S1 versus all the other states. The increase of TE in WS-to-NS and NS-to-WS connections was not as pronounced as in WS-to-WS and NS-to-WS cases. S3 showed relatively low TE such that WS-to-NS of S3 was lower than that of S1 (Fig. 6F, fifth panel from the left; Table 2) and NS-to-NS of S3 was lower than that of S4 (Fig. 6F, sixth panel from the left; Table 2). The findings suggest that desflurane exerts a more substantial effect on WS-to-WS and NS-to-WS connections than WS-to-NS and NS-to-NS connections.
Discussion
The main goal of this work was to determine how cortical unit activity changes with dynamically transitioning brain states under anesthesia. Using an unsupervised machine learning method, we identified five brain states with distinct neuronal spiking behavior. Multiple brain states were observed at a constant anesthetic concentration; and conversely, the same brain state occurred at different anesthetic concentrations. The spontaneous shift of brain states at fixed anesthetic levels suggested that the neuronal network underwent metastability (Bovier, 2006). Recent studies of large-scale brain activity during anesthesia argued that neuronal dynamics may be at equilibrium on short timescales (seconds) but undergo state switching at longer timescales (minutes) (Hudson et al., 2014; Hudson, 2017). In our study, spiking behavior of individual neurons at multiple anesthesia levels also showed metastability consistent with the previous EEG studies.
Spike-based and LFP-based brain states
Brain states obtained by LFP-based clustering method also showed many-to-many correspondence with the level of anesthesia, and resulted brain states similar to those obtained by spike-based clustering, although a complete agreement between the two methods was not established. The difference between the two methods might be explained by the difference in spatial extent that the spikes and LFP signals originate from. The LFP can spread to distant sites reaching >1 cm from its origin (Kajikawa and Schroeder, 2011; Herreras, 2016). Furthermore, the LFP spread may also be state-dependent. Computational modeling suggested that the spatial extent of the region generating the LFP signal is large when distant populations are synchronized (Lindén et al., 2011) (e.g., slow oscillation or burst suppression). In contrast, the spike-based brain states reflected local cortical state around the electrode recording sites, allowing a more direct comparison of the spiking behavior of individual neurons. Last, quantities, such as SR and TNS, cannot be directly captured by LFP.
Synchronous fragmentation and information capacity
The intermittent firing pattern indicated by the increased LV and bimodal ISI was synchronous among the neurons. The intermittency was also reflected in the MUA by an increase in LPBM. This synchronized fragmentation of spike activity estimated by the increase of SUA-MUA correlation was more effective in distinguishing the four brain states (except S3) than all the other examined properties of spike dynamics, suggesting that the synchronously fragmented spike activity is the most pronounced effect of anesthesia (Fig. 6D). Moreover, high SUA-MUA coupling implies that the spike activity of most neurons is constrained by the population activity. Because in this case the entire population acts like a single neuron, information capacity of the population is very low (Tononi, 2004; Izhikevich, 2006). Although high SUA-MUA coupling suggests an increase of shared information, the information content of each neuron is extremely limited. Thus, one could surmise that, in such condition, the animal would be unconscious. This also explains why surgery is preferred in anesthetic states when the EEG displays slow oscillations. Synchronized fragmentation of spike activity has also been reported with other anesthetics, such as ketamine/xylazine (Steriade et al., 1993), urethane (Steriade et al., 1993; Kasanetz et al., 2002; Clement et al., 2008), propofol (Lewis et al., 2012), in addition to desflurane (Vizuete et al., 2014), despite the agents' diverse molecular structure and pharmacological targets. Thus, our study suggests that the synchronously fragmented spike pattern seen with most anesthetics is a common signature of impaired information processing associated with loss of consciousness.
Unlike sleep, anesthesia affects all neurons
In natural sleep, the SR of high-firing neurons substantially decreased while the SR of low-firing neurons was enhanced (Watson et al., 2016). It has been suggested that high-firing neurons appear to be comprised of so-called choristers, which conform to the mean SR of the neuronal population, while the low-firing neurons called soloists respond to stimulation with firing rate changes distinct from that of the population (Bachatene et al., 2015). Specifically, the preferential augmentation of SR of low-firing, stimulus-selective neurons during rapid eye movement sleep has been thought to contribute to an increase of the signal-to-noise ratio of sensory processing (Clawson et al., 2018). The fact that desflurane decreased the SR of virtually all neurons, including those with low baseline firing rate, may be one of the reasons why sensory functions are more completely blocked in anesthesia than in sleep.
Anesthesia facilitates bursting of WS neurons
A recent modeling study of spiking neuronal network demonstrated that burst spikes of individual neurons are more influenced by their presynaptic environment than by their cell type (Tomov et al., 2016). For example, regular spiking neurons could exhibit burst firing because of a network effect. Burst spikes of individual neurons can also shape global network dynamics (C. Y. Li et al., 2009). Nevertheless, the causal relationship between the intrinsic spiking pattern of individual neurons and network synchronization in vivo has yet to be elucidated. While some neurons showed bursting, approximately three-fourths of the neurons were essentially inactive (fired at <1 Hz) in deep anesthesia. Several modeling studies postulated the suppression of metabolic rate in the brain as a key mechanism of anesthetic-induced low-frequency oscillations and burst suppression (Cunningham et al., 2006; Ching et al., 2010, 2012). Another study reported that the occurrence of burst spikes is highly associated with suppression of spike activity, such that hyperexcitable state at the end of suppression period enables an emission of burst spikes (Kroeger and Amzica, 2007). However, it remains unclear how anesthesia almost completely suppresses the majority of neurons while causing burst and synchronized activity in the remaining active neurons. When synaptic inhibition is enhanced by anesthesia, neural network exhibits either slow synchronous firings or isoelectric activity (Lukatch and MacIver, 1996; Lukatch et al., 2005). Synchronous firing may allow neurons to receive enough number of spikes from connected neurons within a short period of time, thereby preventing the rapid decay of spike activity. Neurons with many and strong synaptic connections would generate a large number of synchronized action potentials, and those with weak connections would be inactive. Neurons with strong population coupling (choristers) receive many synaptic inputs from their neighbors (Okun et al., 2015), and show high firing rate both in wakefulness and anesthesia (Bachatene et al., 2015; Clawson et al., 2018). However, the existence of a paradoxical desynchronized state in deep anesthesia suggests the possibility of an alternative scenario in which neurons fire asynchronously while the mean firing rate is profoundly suppressed to a level comparable to the average firing rate during burst suppression. A future modeling study should offer the mechanism of how paradoxical desynchronized state as well as the typical synchronized state can emerge under anesthetized condition.
The paradoxical desynchronized state and consciousness
In the paradoxical desynchronized state S3 mostly found in deep anesthesia (6% desflurane), the mean SR was as low as in S1 (burst-suppression period); however, a small portion of neurons showed distinctly high SR (Fig. 4B,C). Interestingly, the firing pattern of neurons was asynchronous, similar to wakefulness (Fig. 5). Together with the relatively high EMG activity (Fig. 3E,F), this suggested that S3 could be considered a paradoxically aroused state. One may surmise that there was at least a theoretical possibility of transient awareness in this state similar to that proposed for the UP states in slow-wave sleep (Destexhe et al., 2007). Given the unexpected nature of this paradoxical state, replication of this finding together with a more systematic behavioral assessment, such as pupillometry (Shimaoka et al., 2018) or purposeful responses to vibrissal, olfactory, visual/corneal stimuli (Jugovac et al., 2006), should lead to clinically significant implications and advance the general understanding of neuronal network dynamics in different states of consciousness.
Spiking behavior changes monotonically with brain state
It has been widely reported that most properties of large-scale brain activity (EEG, LFP) exhibit a biphasic pattern with anesthetic depth; that is, an initial increase (decrease) of EEG/LFP variable is followed by decrease (increase) as anesthesia deepens. (Kuizenga et al., 1998, 2001; Lee et al., 2017). Likewise, in our study, the low-frequency dominance of LFP (L/H ratio) first increased and then decreased (ignoring S3) as the anesthetic was stepwise withdrawn. Unexpectedly, the spiking behavior did not follow this biphasic pattern; the changes were always monotonic for all spike properties from burst suppression (S1) to full wakefulness (S5) (again, ignoring S3). The reason for this discrepancy is not known, but it may be because of a limitation of measurements at a large-scale level. For example, SR decreases monotonically as the anesthesia deepens, but this cannot be measured directly by EEG or LFP because the latter measurements mostly reflect synchronous population activity. In addition, the theory of complex systems predicts that, for a system consisting of many interacting elements, such as the neuronal network, an incremental change in local variables can lead to abrupt, qualitative change in macroscopic variables, in this case, leading to biphasic macroscopic behavior (Kuizenga et al., 1998, 2001).
In conclusion, the identified brain states displayed degeneracy in their relationship with the anesthetic concentration, suggesting the presence of metastable dynamics. A previously unidentified, transient, paradoxically desynchronized state was found in deep anesthesia. The overall shift in SR distribution by anesthesia is distinct from sleep. The synchronously fragmented spiking in anesthesia appears to be a signature of the state of unconsciousness.
Footnotes
This work was supported by National Institute of General Medical Sciences, National Institutes of Health Award R01-GM056398 and the Center for Consciousness Science, Department of Anesthesiology, University of Michigan Medical School. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. We thank members of the Center for Consciousness Science, University of Michigan Medical School, for valuable comments; and Kathy Zelenock for assistance in laboratory operations and manuscript editing.
The authors declare no competing financial interests.
- Correspondence should be addressed to Anthony G. Hudetz at ahudetz{at}umich.edu
