Intrinsic Functional Connectivity Resembles Cortical Architecture at Various Levels of Isoflurane Anesthesia

Abstract

Cortical single neuron activity and local field potential patterns change at different depths of general anesthesia. Here, we investigate the associated network level changes of functional connectivity. We recorded ongoing electrocorticographic (ECoG) activity from temporo-parieto-occipital cortex of 6 ferrets at various levels of isoflurane/nitrous oxide anesthesia and determined functional connectivity by computing amplitude envelope correlations. Through hierarchical clustering, we derived typical connectivity patterns corresponding to light, intermediate and deep anesthesia. Generally, amplitude correlation strength increased strongly with depth of anesthesia across all cortical areas and frequency bands. This was accompanied, at the deepest level, by the emergence of burst-suppression activity in the ECoG signal and a change of the spectrum of the amplitude envelope. Normalization of functional connectivity to the distribution of correlation coefficients showed that the topographical patterns remained similar across depths of anesthesia, reflecting the functional association of the underlying cortical areas. Thus, while strength and temporal properties of amplitude co-modulation vary depending on the activity of local neural circuits, their network-level interaction pattern is presumably most strongly determined by the underlying structural connectivity.

Introduction

In the absence of specific tasks or stimuli the brain’s spontaneous, ongoing activity is structured on different spatiotemporal scales. Various measures of functional connectivity reflecting different modes of intrinsically generated coupling between different brain areas have been developed (Deco et al. 2011; Siegel et al. 2012; Engel et al. 2013). Among these intrinsic coupling modes (ICMs), amplitude envelope correlations (in the following termed “envelope ICMs”) represent the level of simultaneity of activity fluctuations between neural populations in 2 cortical areas at a particular carrier frequency (Engel et al. 2013). Amplitude envelopes are typically modulated at frequencies below 0.1 Hz and show low dependence on behavioral state (Leopold et al. 2003). In the awake state, high envelope ICMs are found between bilateral homologous areas (Nir et al. 2008; Hipp et al. 2012) as well as within functional areas of one hemisphere (He et al. 2008). The spatial specificity of the observed patterns can be increased by excluding effects of volume conduction through orthogonalization (see Methods; Hipp et al. 2012). Convolution of the amplitude envelope with the hemodynamic response function yields a close estimate of the functional magnetic resonance imaging (fMRI) blood-oxygen level dependent (BOLD) signal (Liu et al. 2011; Martin 2014). Accordingly, envelope ICMs across 2 brain areas observed with electrophysiological methods are closely linked to the correlations of their BOLD signals (He et al. 2008; Schölvinck et al. 2010; Keller et al. 2013). Thus, for instance, networks defined by envelope ICMs recorded with magnetoencephalography (MEG) show a high similarity to fMRI-resting networks (de Pasquale et al. 2010; Hipp et al. 2012).

Ongoing activity is severely altered during general anesthesia, that is, the reversible loss of consciousness induced through the administration of hypnotic agents. Patterns of cortical activity under general anesthesia change with the level of the hypnotic agent. After loss of consciousness, the EEG first shows continuous patterns with increased delta- and alpha- and decreased beta power (Brown et al. 2010). In deeper anesthesia, brain activity changes to burst suppression, which is characterized by an isoelectric line interrupted by bursts of activity (Llinás and Steriade 2006; Brown et al. 2010). The frequency of bursts decreases with deepening anesthesia, until complete suppression is reached (Brown et al. 2010). Network-level changes caused by anesthesia have been examined using various recording techniques and analysis approaches. Studies focusing on the awake-anesthetized transition have suggested that loss of consciousness is associated with global changes in functional connectivity (Cimenser et al. 2011) and distinct alterations of phase-amplitude coupling (Lewis et al. 2012; Mukamel et al. 2014). Large-scale envelope ICMs under different anesthesia depths have been studied using fMRI. While some of these studies observe an increase of connectivity with deepening anesthesia (Liu et al. 2013), others report a breakdown of coupling (Bettinardi et al. 2015). Micro-scale networks on the cellular level have been characterized using extra- and intracellular microelectrode recordings under different concentrations of anesthetic agents (Llinás and Steriade 2006). Nevertheless, little is known about the corresponding connectivity patterns in mesoscale cortical networks consisting of spatially and functionally related cortical areas.

Here, we characterize envelope ICMs in visual, auditory and parietal areas of the ferret at different levels of isoflurane anesthesia and a constant nitrous oxide administration. We recorded local field potentials from the cortical surface using a custom-built 64-contact electrocorticographic (ECoG) array. Unlike EEG, these recordings are not affected by the distortion of the electrical field in space and frequency by the surrounding tissues (Buzsáki et al. 2012). The isoflurane concentration was varied between 0.4% and 1.6% in steps of 0.2%. For each anesthesia level, envelope ICMs were quantified in different frequency bands using orthogonalized amplitude correlation (Hipp et al. 2012). Similar coupling patterns were visible at different isoflurane concentrations across animals, reflecting inter-individual differences in drug sensitivity. We therefore applied hierarchical clustering to all connectivity matrices of each frequency band and identified similar coupling-defined network states across animals. The resulting 3 clusters corresponded to light, medium, and deep anesthesia. Overall, the strength of envelope ICMs increased across anesthesia depths. Normalized connectivity patterns, however, remained similar across anesthesia depths, resembling in their topology known patterns of structural connectivity in the cortex.

Methods

Data presented in this study were collected from 6 adult female ferrets (Mustela putorius). All experiments were approved by the independent Hamburg state authority for animal welfare (BUG-Hamburg) and were performed in accordance with the guidelines of the German Animal Protection Law.

ECoG-Array

We employed a micro-ECoG-array co-developed with the University of Freiburg (ECoG-array; IMTEK, Freiburg) covering a large portion of the posterior, parietal, and temporal surface of the left ferret brain hemisphere (Rubehn et al. 2009). The probe was designed to record from posterior early and higher visual areas, auditory, somatosensory, and parietal areas, located in the medial, suprasylvian, and lateral gyrus, respectively (Manger et al. 2002; Bizley et al. 2007). Sixty-four platinum contacts (ø: 250 μm; 0.05 mm², 7–25 kΩ @ 1 kHz) were arranged equidistantly (1.5 mm) in a hexagonal manner on a 10 μm thin polyimide foil. The approximately 12 × 10.5 mm sized foil was subdivided in 3 “fingers”, allowing for a smooth adaptation to the slightly curved surface of the ferret cortex (Fig. 1A).

Surgery

Animals were initially anesthetized with an injection of ketamine (15 mg/kg) and medetomidine (0.08 mg/kg). After tracheotomy, a Y-shaped glass tube was placed in the trachea to allow artificial ventilation and the administration of isoflurane (0.4–1.6%, 1:2 N₂O–O₂ mix; see in more detail below). N₂O was added for its anesthetic and analgetic effect. Before fixing the animal’s head in the sterotaxic frame, a local anesthetic (Lidocain, 10%) was applied to the ear cannula. To prevent dehydration throughout the experiment, a cannula was inserted into the femoral vein to deliver a continuous infusion of 0.9% NaCl, 0.5% NaHCO₃, and pancuronium bromide (6 μg/kg/h). Body temperature was monitored rectally and automatically maintained at 38 °C with a homeothermic blanket. After the initial sedation, 1.0% isoflurane was added to the ventilated air to maintain anesthesia throughout the following surgical procedures. Stereotaxic ear- and eye-bars were removed after surgery before the level of anesthesia was decreased. An anesthetic monitoring device (Datex-Ohmeda S/5, Datex-Ohmeda Medical Instrumentation Inc., Tewksbury, USA) was used to continuously measure inspiratory and expiratory concentrations of isoflurane (Supplementary Table 2), MAC, N₂O, O₂, and CO₂ throughout the course of the experiment. All data were digitized and logged to disk via the electrophysiological recording systems.

To get access to the left parieto-occipital cortex, first, the temporalis muscle was reflected laterally. A saline-cooled and soft-tissue-safe micro-saw (Mectron, Piezosurgery) was used to create a craniotomy (~13 × 11 mm) over the left posterior cortex. After careful removal of the dura, the ECoG-array was placed on the cortex. The ECoG-array was covered with artificial dura (Lyoplant/Braun) and the previously extracted bone flap was returned into place and fixed with KWIK-CAST (World Precision Instruments, WPI), a tissue-safe silicon-based two-component glue. Finally, the temporal muscle and the skin were flapped back in anatomical position and provisionally fixed with a suture. This ensured that the cortex would reside in physiological cortico-spinal fluid and avoided cooling of the cortex to nonphysiological temperatures throughout the course of the experiment (Kalmbach and Waters 2012).

After the end of the whole experimental session (~24–36 h), the animal was deeply anesthetized with 5% isoflurane and sacrificed.

Data Acquisition

The continuous stream of neuronal data (~4.5 h) was processed with an AlphaLab “SnR” electrophysiological recording system. Signals were amplified 200-fold, hardware bandpass filtered (0.5 Hz; 2-pole—8 kHz; 3-pole) and digitized at 44.6 kHz sampling rate. It was then immediately software lowpass filtered (347 Hz), downsampled by a factor of 32 (1.4 kHz) and saved to disk. Concurrently, vital signs and anesthesia parameters as described above were sampled with the same system (2.8 kHz) and also stored to disk.

Data acquisition started between 5 h and 9.5 h after the initial ketamine anesthesia following at least 4 h of isoflurane ventilation. Before the recordings started, the isoflurane level was kept at 0.4% for at least 1 h. During the recording period the isoflurane concentration was changed every 20 min in steps of 0.2%, first increasing up to a level of 1.6% and then decreasing again down to 0.4%. This resulted in 13 different conditions of “increasing” and “decreasing” levels of isoflurane. After each dosage change, an equilibration of isoflurane concentration between the inspiratory and expiratory airflow was reached after about 2–4 min. The administered isoflurane concentrations differed on average less than 0.001% from the preset values (STD: 0.018%; see Supplementary Table 2).

Data Analysis

All data analysis was performed with Matlab (The MathWorks Inc.), partially using the Fieldtrip (Oostenveld et al. 2011) and the Chronux toolbox (http://chronux.org/; Mitra and Bokil 2008). Offline data were bandpass filtered in a phase preserving manner between 2 Hz and 300 Hz and further downsampled to 700 Hz. Notch filters were applied at 50 Hz, 100 Hz, and 150 Hz (1 Hz width, 4-pole butterworth). Broken contacts, as indicated by low signal amplitudes or high-amplitude noise, were excluded from further analysis. All data were visually inspected before further analysis to exclude ECG- or other electrical artifacts.

We computed the power spectra with a multitaper (4/1) analysis of all channels and conditions between 2 Hz and 256 Hz at 6 frequency bands per octave. Power in the classical EEG frequency bands—theta (4-8 Hz), alpha (8–13 Hz), beta (13–30 Hz), low gamma (30–64 Hz), and high gamma (64–128 Hz)—was calculated as the average of the respective frequency ranges.

The underlying time–frequency representation was calculated based on Morlet wavelets. The temporal resolution was set to 1/2 of the maximum frequency’s cycle duration (7.8 ms), while the initial frequency resolution was at 8 logarithmically spaced bands per octave, resulting in 56 Morlet wavelets of 7 cycles. For further analysis, we calculated the mean correlation coefficients of each channel pair in the theta, alpha, beta, low-gamma, and high-gamma band as specified above.

We analyzed the spectral properties of the amplitude modulation by calculating a multitaper spectrum of the amplitude envelope of each channel. To estimate the amplitude correlation in a given modulation frequency band, we orthogonalized one signal with respect to the other as described above and applied a bandpass filter to both amplitude time courses before calculating the correlation coefficient.

Assignment of Contacts to Cortical Areas

To account for inter-animal differences in the position of the ECoG-array and thus make cross-animal comparisons possible, area-level correlation matrices were calculated for cross-animal analysis. Using the sulci as landmarks, a schematic representation of the array position for each experiment was generated relative to a functional map of the ferret cortex (Bizley et al. 2007). We then assigned each contact to a cortical area by a “winner-gets-all” method and calculated the average connectivity values between areas. See Supplementary Fig. 1 and Supplementary Table 1 for more precise description of the method.

Hierarchical Clustering

For each frequency band, we applied agglomerative hierarchical clustering as implemented in Matlab to the area-level coupling matrices of all animals and conditions to identify common network states. To ensure the robustness of the result, we compared the emerging clusters for 2 different distance measures (“pdist”: “euclidian” and “cityblock”) and using 2 linkages (“ward” and “average”; Ward 1963). An example of a clustering result is shown in Supplementary Fig. 2.

We chose to partition the dataset of 78 conditions (6 subjects × 13 isoflurane levels) into 3 clusters of 44, 19, and 15 members, respectively. At this level, the consistency of the defined clusters across algorithms was high, and the dendrograms showed a clear separation of the clusters from the next level. Further branching would have led to stronger, inconsistent parcellation with groups of very few members (<6), such that the comparability would suffer from statistically noisy results. The results presented here are based on “average” clustering by “euclidean distance”.

Normalization of Correlation Values

To examine correlation differences which could not be explained by the distance between contacts, each correlation coefficient in the contact-level connectivity matrix was re-expressed relative to the correlations of all contact pairs with the same inter-contact distance. Contact distances were rounded to the first decimal place. We only considered distances which were represented by at least 25 contact pairs in all animals. All distances above the first one not fulfilling this criterion were excluded. This resulted in 18 distances between 1.5 mm and 9.8 mm (Supplementary Fig. 3). The z-score z = (x − μ)/σ was calculated for each of the correlation values.

Burst Detection and Avalanche Definition

We performed an avalanche-style analysis of epochs of burst suppression (Beggs and Plenz 2003; Shew et al. 2009). The time-domain data of the individual channels were segmented in periods of “burst” and “suppression” using a modified version of the method proposed by Lewis et al. (2013). We analyzed data between 1.0% and 1.6% isoflurane, since under these concentrations all animals displayed burst-suppression patterns. All data were z-transformed. For each animal, periods of unambiguous suppression were labeled in the first 120 s of the recording from one contact under 1.0% isoflurane. Data were filtered above 8 Hz and the instantaneous amplitude was calculated via a Hilbert transform and smoothed with a moving average filter of a width of 200 ms. The threshold for bursts epochs was set 4 standard deviations above the mean of manually labeled suppression epochs. Suprathreshold epochs of least 200 ms were defined as bursts. Suppression phases were defined as epochs in which the value fell below the threshold for at least 200 ms. The result was cross-checked by visual inspection. Avalanches of neural activity were defined as epochs of burst activity propagating across the grid. They were by definition separated by periods of global suppression, in which none of the channels was in “burst mode”. Channels which displayed bursting activity for more than 95% of the recording time were discarded before avalanche definition.

Results

We recorded ongoing cortical activity from the left parieto-temporo-occipital cortex of 6 anesthetized female ferrets using a custom-made 64-channel subdural electrode array (Fig. 1A). Animals were ventilated with a mixture of isoflurane and N₂O. We determined the changes of the brain network connectivity at different isoflurane levels.

Isoflurane Increases Low-Frequency and Decreases High-Frequency Power

We examined the effect of deepening anesthesia on power change relative to the first data segment in the typical frequency bands. For this purpose, we calculated the modulation index (MI) (A − B)/(A + B) within each contact—with “A” being the power of the segment of interest, “B” of the first data segment—scaling the values to a range from −1 to 1. With increasing isoflurane concentration all animals showed a relative increase in theta, alpha, and beta power and a relative decrease in gamma power (Fig. 1D).

We next assessed whether there were any topographical differences in the power change occurring from light to deep anesthesia. We pooled the power values across 5 functional groups as follows (cf. Supplementary Fig. 1; Bizley et al. 2007): early visual (areas 17, 18, 19), higher visual (area 21, SSY), auditory (areas A1, AAF, ADF, PPF, PSF), posterior parietal (areas PPr, PPc), and somatosensory (areas SII, SIII). For these areal groups, we compared the MI at 1.6% isoflurane. For all frequency bands, the power change relative to 0.4% isoflurane was significant for all functional groups (Komolgorov–Smirnov test (KS-test), Bonferroni–Holm corrected, P < 0.05). However, there were only a few significant differences between the functional groups (Fig. 2). Overall, the observed relative change in power occurred relatively homogeneously across the whole temporo-parieto-occipital cortex. This result is in accordance with other studies showing an occipito-frontal power shift under anesthesia rather than local differences between the areas of the occipital cortex (Cimenser et al. 2011; Supp et al. 2011; Sellers et al. 2013).

Cluster Analysis of Envelope ICMs Reveals 3 Anesthesia-Dependent Network States

For each 4-min data segment, we calculated orthogonalized amplitude correlations between all channels. This measure excludes the zero-phase-lag component of the 2 signals and is thus not susceptible to potential effects of volume conduction (Hipp et al. 2012). Correlation values were calculated on a wavelet-transformed time–frequency representation and then pooled to the following frequency bands: theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), low gamma (30–64 Hz), and high gamma (64–128 Hz).

For these frequency bands, we then computed coupling matrices reflecting the envelope ICMs at different anesthesia depths. It has long been known from clinical settings that different anesthetic concentrations are needed to elicit the same anesthetic effect in different individuals. We hypothesized that due to inter-individual differences in responsiveness to the anesthetic agent, similar coupling patterns might emerge in the individual animals under different concentrations of isoflurane. Therefore, we used hierarchical clustering to identify similar coupling patterns for each frequency band across animals and anesthesia depths.

For each frequency band, we calculated the envelope ICM coupling matrix for each data segment and animal. Amplitude correlation values were averaged within 14 cortical areas covered by the ECoG-array in all animals. This approach eliminated inter-individual differences of the contact positions and the number of recording sites per cortical area, allowing cross-animal comparison. For each frequency band, the matrices were vectorized and subjected to hierarchical clustering (Supplementary Fig. 2). Based on the dendrograms, we chose to partition the dataset into 3 clusters corresponding to network states defined by functional connectivity under light, medium, and deep anesthesia (Fig. 3A). The deep-anesthesia connectivity state was consistent across all frequency bands but appeared in only 4 of the 6 animals. In contrast, the light and medium anesthesia connectivity pattern were observed in all animals. They were nearly identical for the theta, alpha, and beta bands. In the low- and high-gamma band, the medium anesthesia connectivity state occurred at lower isoflurane concentrations than in the lower frequency bands. This indicates a frequency band related difference in network-level effects of isoflurane.

For each anesthesia network state, we calculated the area-level coupling matrices averaged across animals (Fig. 3B). Within each frequency band, global correlation strength increased from the light to the deep-anesthesia network state for all fields of the coupling matrix (light: 0.08, medium: 0.23, deep: 0.41; HedgesG: light–medium: −1.57, medium–deep: −1.66, light–deep: −3.25). A KS-test between the clusters confirmed that this difference was significant for all connections in the matrix (Bonferroni–Holm corrected, P < 0.05).

The structure of the coupling matrices suggested that particularly high correlations occurred among functionally closely related areas. To test this hypothesis, we pooled the single-contact correlation values into 5 functional groups (early visual, higher visual, posterior parietal, auditory, and somatosensory) and compared correlation values within these functional groups with those across functional groups (Fig. 3C). Within functional group correlations were consistently higher than across functional group correlations for all network states and frequency bands (KS-test, Bonferroni–Holm corrected, P < 0.05).

Network States are Associated with Changes of Burst-Suppression Characteristics

Inspection of the raw ECoG signal for the 3 connectivity states revealed patterns known from single-cell and EEG recordings (Steriade et al. 1994). The data segments in the light anesthesia state mostly revealed continuous activity patterns, which in some animals were overlaid by bursts of higher-frequency activity. For the medium anesthesia state frequent bursts with short periods of suppression were observed and for deep-anesthesia groups of bursts interrupted by long suppression intervals (Supplementary Fig. 4). At the lowest isoflurane dose of 0.4%, 3 animals already showed burst-suppression patterns. At the highest concentration of 1.6%, all animals had entered burst suppression, 4 displaying the deepest burst-suppression state. Spectral analysis of the anesthesia states showed an increase of low-frequency power and a decrease of high-frequency power in the medium and deep-anesthesia state (Fig. 1D, see also Supplementary Fig. 5).

We further characterized the burst-suppression segments of the dataset using the approach of Lewis et al. (2013). Data segments recorded under isoflurane concentrations of 1.0% and higher were included since under these conditions all animals displayed burst-suppression patterns.

We calculated the average burst frequency, mean burst length, and mean suppression length for each data segment. Burst frequency and suppression length showed a distribution similar to the medium and deep-anesthesia state as obtained through our connectivity-based analysis (Supplementary Fig. 6A, B). Burst lengths showed no clear dependence on isoflurane concentration (Supplementary Fig. 6C).

We quantified the high visual similarity of connectivity-based anesthesia states, burst frequency, and suppression length through linear correlation. For all frequency bands, the matrices were significantly correlated at P < 0.001 (Bonferroni–Holm corrected, Supplementary Table 3). Thus, the anesthesia-dependent changes of burst-suppression characteristics of local activity were closely linked to the states observed on the network level. Deep-anesthesia states were characterized by a decrease of burst frequency and an increase of the length of suppression episodes.

Overall Strength of Envelope ICMs Shows Dependence on Distance and Frequency Range

To capture general properties of envelope ICMs within the 3 network states, we examined their dependence on frequency and spatial distance of the contacts. For this purpose, we calculated the mean correlation values across all electrode combinations and within 6 contact distance categories for the 3 clusters (Fig. 4). For all 3 anesthesia-dependent network states, correlation showed a drop-off as a function of distance and frequency. In the light anesthesia network state, envelope ICMs showed a peak in the alpha-frequency range and continuously diminished with increasing frequency. During medium depth anesthesia, the highest correlation values were observed in the beta band. During deep anesthesia, mean strength of envelope ICMs was similar in the theta, alpha, and beta bands, and only dropped within the gamma and high-gamma band.

Spatial Structure of Envelope ICMs Remains Robust to Anesthesia Depth and Distance Effects

As described above, the most obvious effect of anesthesia was an overall increase of envelope ICMs between and within all areas. The relative spatial pattern within the connectivity matrix, however, appeared similar within each frequency band across anesthesia-dependent network states. To further examine this aspect, we re-expressed each value within each area-level coupling matrix as a z-score relative to all correlation values of the particular coupling matrix. The resulting mean z-scored connectivity matrix (Fig. 5A), representing only the spatial structure of the envelope ICMs separated from the absolute coupling strength, remained similar across anesthesia depths. The topography of the envelope ICMs was similar for all frequency bands.

Since the functional groups comprised neighboring cortical areas, the general decrease of amplitude correlation with distance would induce a bias to higher intra- compared with cross-group correlations. Therefore, we tested whether the spatial patterns visible in the coupling matrices could be explained exclusively by the lower distance of the contacts comprised in areas and functional groups. For this purpose, we re-expressed each correlation coefficient of the contact-level coupling matrix as a z-score relative to all correlations between contact pairs with the same distance (Supplementary Fig. 3). On the resulting values, we recomputed the coupling matrices for the 3 network states defined above. As shown in Fig. 5C, the topographical patterns observed previously were still largely preserved. Furthermore, within-group z-scores were still higher than across-group z-scores (KS-test, Bonferroni–Holm corrected, P < 0.05). This suggests that higher envelope ICMs between functionally related areas cannot exclusively be attributed to their spatial proximity.

Envelope Modulation Frequency Changes with Anesthesia Depth

To further investigate the temporal dynamics of amplitude envelopes in the 3 network states, we calculated the mean power spectrum of the amplitude envelope for each frequency band and each state (Fig. 6A). The spectrum of the signal amplitude was relatively flat for all frequency bands in the light anesthesia network state. During medium anesthesia, a peak around 0.7 Hz was visible, which was most pronounced in the gamma band. For the deep-anesthesia state, the peak frequencies of the amplitude modulation were below 0.1 Hz.

This change in modulation frequency raised the question whether amplitude correlations differed for specific modulation frequencies, and whether this dependency changed with anesthesia. We therefore bandpass filtered the time course of the amplitude envelope in eleven frequency bands between 0.04 Hz and 2 Hz and calculated the mean correlation within each of them for each frequency and each cluster (Fig. 6B). Correlation at all modulation frequencies increased with anesthesia. During light anesthesia, the resulting spectrum was flat. For the medium anesthesia state, it showed a peak around 0.5 Hz, which was more pronounced in the beta- and gamma-band than in the theta- and alpha-band. During the deep-anesthesia network state, correlation showed a constant decay with modulation frequency.

Deep-Anesthesia Network State is Characterized by Global Spread of Avalanches

To link our analysis of burst-suppression epochs to the global network dynamics we used an avalanche-based approach, established in the analysis of activity of neuronal and slice cultures (Beggs and Plenz 2003; Shew et al. 2009). Avalanches are cascades of activation propagating through a network. Once one node of the network is activated (in our case, when the area of cortical tissue under one ECoG contact switched into “burst mode”), it can activate other nodes of the network, thus initiating a cascade of activity. This cascade can reach from the activation of only one neighboring node to the activation of all nodes of the network. By definition, individual avalanches are separated by epochs of global “suppression” across all contacts of the ECoG.

Data were segmented into burst- and suppression epochs as described above. We included data recorded at isoflurane concentrations of 1.0% and above, under which all the time-series data showed burst-suppression patterns. Avalanches were defined as sequences of noninterrupted, overlapping burst activity in at least one of the channels of the grid. They were by definition separated by epochs of global suppression, in which none of the channels was in “burst-mode”. The number of avalanches per data segment of 4 min was 87.2 ± 37 std. As a parameter relatively independent on the small number of avalanches in each data segment, we chose the ratio of avalanche activity time to suppression time to further classify the data. This ratio showed a pattern similar to the connectivity states (Fig. 7A). Hierarchical clustering of the values showed 2 states almost identical to the medium and deep-anesthesia connectivity state (Fig. 7B). Linear correlation between the state matrices was significant for all frequency bands at P < 0.001 (Supplementary Table 3).

Avalanches differed in their spatial extent, ranging from activity in at least one channel without further spread of activity to activation of all channels during one avalanche. We analyzed the distribution of the number of channels active during avalanches for the 2 anesthesia states. For both medium- and deep-anesthesia state, the shape of the distribution was similar. It had a small peak at bursts including only few channels and a larger peak at global bursts including almost all channels of the grid. In deep anesthesia, the peak for global bursts was steeper and higher than in medium anesthesia (Fig. 7C). Significance was tested using a bootstrap statistic. We randomly shuffled the data segments into clusters of similar size to the original ones 10 000 times. For each randomization, we calculated the area between the channel number distributions of the 2 random clusters. We then compared the original area between the curves (see Fig. 7C) with the resulting distribution of areas from the randomizations; they differed significantly (P = 0.0018; Supplementary Fig. 7C). Distribution of avalanche length showed a decay from smaller to larger values in both medium and deep anesthesia (Supplementary Fig. 7A) without a significant difference between the states (P = 0.063; Supplementary Fig. 7B).

Overall, avalanche-based analysis of the burst-suppression data segments showed the same anesthesia states as the connectivity-based approach. Under deep-anesthesia avalanches showed a higher tendency to spread globally across the whole cortex, which coincides with the high global connectivity observed in this state.

Discussion

This study shows, using the ferret cortex as a model, that quantification of envelope ICMs allows to characterize network states emerging across different concentrations of isoflurane anesthesia under nitrous oxide/oxygen-ventilation. Based on hierarchical clustering of the coupling matrices, we could differentiate 3 network states under light, medium, and deep anesthesia. The particular anesthetic dosage at which a state was reached varied across animals, reflecting individual differences in the response to the anesthetic. Overall, our results suggest that envelope ICMs in ongoing activity reflect the functional organization of cortical networks in their topographical structure. Although the average coupling increased with anesthesia depth, the topography of functional networks was conserved across different network states. The dominant modulation frequencies of the amplitude envelopes shifted across anesthesia depths, and so did the frequency bands at which the strongest envelope ICM was observed. Anesthesia states classified by envelope ICMs were associated with different dynamics of local neural population activity visible in the electrocorticogram. The light anesthesia state with low global correlations was characterized by continuous activity patterns, the high correlation levels of the deep-anesthesia state were related to low-frequency burst suppression.

Methodological Considerations

It has to be emphasized that activity we observe is not exclusively a result of isoflurane anesthesia, but as well of the co-administered nitrous oxide. The network we examined consists of early as well as higher auditory, visual, and somatosensory areas, thus comprising a relatively large fraction of the cortex. However, since frontal and occipital areas show different reactions to anesthetic agents in humans (Cimenser et al. 2011; Supp et al. 2011) as well as in the ferret (Sellers et al. 2013), our results cannot be generalized to the entire neocortex without further experimental evidence.

We partition the signal into frequency bands that can be differentiated in ferrets during the awake state and sleep (Stitt et al. 2017) and are as well observed in primates, cats, and rodents. This provides a common framework for comparison across states, but at the cost of neglecting possible state-specific boundaries. The approach for quantification of the envelope ICMs used here includes the orthogonalization of the phases at each frequency which eliminates effects of volume conduction and common input, but comes at the cost of also removing additional physiological in-phase activity, which plays an important role in neural processing (Singer 1999; Engel et al. 2001). As shown before, this also implies that our approach to most underestimates functional connectivity (Hipp et al. 2012). To evaluate the effects of orthogonalization on connectivity, we also calculated the envelope ICMs using non-orthogonalized amplitude correlations. The clustering result was similar and the ICMs displayed similar coupling patterns with a slight overall increase of connectivity (Supplementary Fig. 8). This finding points against the possibility that orthogonalization eliminates crucial physiological characteristics of network activity under anesthesia. Furthermore, it shows that the ECoG signal is only minimally influenced by volume conduction. Our state classification by hierarchical clustering interprets deepening anesthesia as switches between clearly distinguishable states. Due to the relatively large logical distance between the 3 clusters visible in the dendrograms, we believe this approach to the data is justified, even though it ignores minor within-state differences of the coupling patterns. Furthermore, our approach does not take into account spontaneous fluctuations in network connectivity strength over time (de Pasquale et al. 2010).

Changes of Local Dynamics and Spectral Characteristics Under Anesthesia

The ECoG signal can predominantly be attributed to the synchronous population activity of cortical pyramidal neurons (Buzsáki et al. 2012). We observed a change of the ECoG signal from continuous activity under light anesthesia to burst-suppression patterns under higher isoflurane concentrations, an effect long known from surface EEG recordings in clinical settings (Brown et al. 2010). Spectral characterization of local activity revealed changes known from previous studies that have employed isoflurane anesthesia (Sellers et al. 2013). With increasing isoflurane concentration we observed a relative increase in theta, alpha, and beta power and a relative decrease in gamma power. These changes in the local dynamics relate to potentiation of GABAergic inhibition in the thalamo-cortical system (Lukatch et al. 2005). At the single-cell level, bursts correspond to epochs of phasic depolarization, while periods of suppression correspond to electrical silence resulting from anesthetic-mediated suppression of glutamatergic transmission (Steriade et al. 1994; Lukatch et al. 2005).

The spectrum of the amplitude envelopes changes across anesthesia depth, likely representing the drastic alteration of activity pattern visible in the surface ECoG signal. The continuous spontaneous activity of the low anesthesia state was reflected in “envelopes”, not raw or band-passed activity, predominantly modulated at frequencies below 0.1 Hz, similar to those described for the awake brain (Leopold et al. 2003; Nir et al. 2008). In medium anesthesia, the “envelope spectrum” showed an overall power increase with a peak around 0.6 Hz especially for the alpha-, beta-, and gamma-band, likely a spectral fingerprint of the emerging burst-suppression patterns. The further increase of modulation frequency, especially below 0.1 Hz in deep anesthesia, might be the correlate of the sparser occurrence and longer duration of bursts, however, a more systematic examination of burst properties would be necessary to verify this assumption. The correlation of the envelopes is highest at the dominant modulation frequencies of the particular state.

Envelope ICMs and Cortical Network State

Complementary to the states of local neuronal population activity visible in the patterns of the raw ECoG signal, we defined network states through amplitude envelope correlations. With deepening anesthesia we observed a global increase of coupling strength with constant relative spatial patterns. The spatial specificity of envelope ICMs we observed across anesthesia depths has been described before in the awake state and sleep. Homologous areas display high interhemispheric correlations in envelope ICMs calculated from MEG (Hipp et al. 2012). Within the same hemisphere, amplitude correlations between ECoG contacts within the sensorimotor cortex have been observed to be higher than those to areas of other modalities in the awake state, slow-wave-, and REM-sleep (He et al. 2008). On a finer spatial scale, the resemblance of tonotopy by envelope ICMs has been shown in the macaque (Fukushima et al. 2012). To our knowledge, our study is the first to examine high-resolution envelope ICMs at different depths of isoflurane anesthesia.

In the medium- and deep-anesthesia network states, which were characterized by burst-suppression patterns, cortical activity was not equally coupled across the whole set of recorded areas, but amplitude correlation was still higher for neighboring cortical sites than for distant ones. Correspondingly, synchrony and similarity of bursts have been shown to depend on cortical distance. Bursts occur nearly simultaneously in adjacent cortical areas (Contreras et al. 1997). However, large-scale ECoG recordings in humans have shown that burst suppression is not completely homogeneous across the whole cortex (Lewis et al. 2013). First, not all areas enter and exit burst suppression simultaneously. Second, under burst suppression, not all bursts are global, and the global ones are less synchronized for more distant cortical areas. These results are well compatible with our observation that even under deep burst suppression, the cortex is not homogeneously synchronized, but the spatial and functional structure of the network is preserved.

In the awake brain, amplitude envelopes and their correlations have been suggested to remain stable across behavioral states (Leopold et al. 2003). Our results show that anesthesia modulates some of their basic characteristics, while others remain constant: While mean correlation increases with anesthesia depth, it shows a negative dependence on distance between recording sites in all 3 states. Similarly, a drop-off of envelope ICMs for all frequency bands as a function of distance has been described in the awake rhesus monkey (Leopold et al. 2003). Furthermore, we observed a drop-off of correlation as a function of frequency, which suggests a more local amplitude co-modulation of higher frequencies. In contrast, in the awake brain particularly high coherence values have been reported for amplitude correlations of the gamma band (Leopold et al. 2003). Whether this is due to the different scale of our recordings, which cover a larger area of much smaller ferret brain, caused by the orthogonalization, or is a property of the anesthetized brain remains open to further investigation.

Relation to fMRI Studies on Functional Connectivity Under Anesthesia

To our current understanding, connectivity revealed by amplitude correlation of electrophysiological signals is comparable to that defined by BOLD correlation in resting-state fMRI (Leopold and Maier 2012; Engel et al. 2013). The low-frequency modulation of neural activity represented in the amplitude envelope is tightly linked to the temporal dynamics of local cerebral blood flow and oxygen extraction mirrored in the BOLD signal via the hemodynamic response function (Martin 2014). ECoG amplitude- and BOLD-derived correlations of ongoing activity are highly correlated (He et al. 2008; Schölvinck et al. 2010; Keller et al. 2013), and resting networks defined by independent component analysis (ICA) of MEG power envelopes resemble those known from fMRI studies (de Pasquale et al. 2010, 2012). The close association of ECoG amplitude envelopes and cerebral blood flow via the hemodynamic response function has been described in isoflurane-anesthetized rats and is preserved in the anesthetized state (Liu et al. 2011).

However, it should also be noted that the comparability of fMRI studies to electrophysiological results is limited by some methodological differences. Most resting-fMRI studies segment activity patterns into distinct networks using spatial ICA or examine networks of few, sometimes distant seed regions. Due to the limited temporal resolution of fMRI, they typically do not consider frequencies above 0.1 Hz, which we show to contribute to amplitude correlations under anesthesia. For the same reason, these studies usually do not examine the temporal dynamics of the signals.

Since different anesthetics potentially differ in their effects on brain connectivity (Brown et al. 2010; Williams et al. 2010; Jonckers et al. 2013; Grandjean et al. 2014), we confine our discussion to studies using isoflurane. It has been shown that the spatial specificity of fMRI-coupling is preserved under isoflurane anesthesia (Vincent et al. 2007; Hutchison et al. 2010). An increase of functional connectivity with deepening anesthesia was observed in a study examining BOLD patterns in the sensorimotor cortex of the rat under 3 isoflurane concentrations (1.0%, 1.5%, 1.8%) (Liu et al. 2013). Distinct sensorimotor subnetworks were distinguishable under the lowest isoflurane concentration and BOLD correlations decreased as a function of distance. With increasing dosage, spatial specificity of correlation patterns decreased and high brain-wide correlations emerged. The same study described a similar effect on power correlations from 3 recording sites, again in accordance with our own findings (Liu et al. 2013).

These strong long-range BOLD correlations observed under deep anesthesia have been shown to be associated with electrophysiological burst suppression and to breakdown under isoflurane concentrations leading to an isoelectric local field potential (Liu et al. 2011). High global correlations have been described under 1.5% isoflurane, which is likely comparable to our medium- or deep-anesthesia state (Kalthoff et al. 2013). While there are no previous electrophysiological data on envelope ICMs under isoflurane anesthesia, increased interhemispheric amplitude correlations have been observed during sleep (Nir et al. 2008). Although these findings are well compatible with our results, other studies report a breakdown of functional connectivity with increasing isoflurane concentration (Wang et al. 2011; Hutchison et al. 2014), and similar observations have been made under other anesthetics (Lu et al. 2007; Bettinardi et al. 2015).

Taken together, it is difficult to derive a consistent picture of resting-fMRI networks under anesthesia from the available literature. These contradictory findings might be due to variations in experimental and analytical methodology. Experimental protocols differ in species, regions between which coupling is calculated, isoflurane concentrations, and co-medication. This makes an assignment of similar brain states across studies very difficult. In addition, neurovascular coupling and thus the BOLD signal are altered by blood CO₂- and O₂-level (Nasrallah et al. 2015), a mechanism by which different ventilation protocols might influence the observed coupling patterns. Additionally, different fMRI preprocessing strategies seriously influence the observed functional connectivity. Since high correlations in deep anesthesia are based on synchronous activity across large brain areas, they disappear if the global brain signal is removed in the preprocessing of the fMRI data, giving the opposed impression of sparsely correlated networks (Kalthoff et al. 2013; Liu et al. 2013).

Avalanche Analysis of Burst-Suppression

The high connectivity levels observed in the deep-anesthesia state were associated with a higher tendency of avalanches to spread across the whole network. Avalanches of multi-unit spikes have been characterized in recordings of spontaneous activity from neuronal cultures and cerebral slices cultures (Beggs and Plenz 2003; Shew et al. 2009). The authors report that avalanches in these networks follow power-law distributions in the number of active electrodes and duration. Their further findings suggest that the network operates in a critical state with a branching coefficient close to 1, which would maximize the processing capacities of the network (Haldeman and Beggs 2005). Notably, reduction of inhibition by a ɣ-aminobutyric acid (GABA)_A-antagonist causes a peak a shift from a power-law-distribution to a bimodal distribution (Beggs and Plenz 2003; Shew et al. 2009). We observed similar effects under burst-suppression caused by the GABA_A-agonist isoflurane. This can likely be explained by the fact that isoflurane also acts on a variety of other ion channels.

Envelope ICMs and Cortical Network Function

At first glance, it might seem counterintuitive that deepening anesthesia leads to increased functional connectivity. However, increased envelope ICMs in deep anesthesia are associated with the rather stereotypic local patterns of burst suppression, which are highly synchronized throughout the network, as already suggested by Liu et al. (2011). Both the more stereotyped local dynamics and the overall increase in functional connectivity are in line with previous hypotheses on the loss of consciousness under deep anesthesia (Alkire et al. 2008; Supp et al. 2011). It has been suggested that anesthesia has 2 effects on the cortex: a breakdown of integration, and the entrainment to stereotypic activity patterns leading to a loss of information (Alkire et al. 2008). Our data agree with this notion and support the following interpretation. While the topographic patterns of the envelope ICMs remain stable, the signals propagating across the network change drastically with the transition into deeper anesthesia. In the awake state and under light anesthesia, dynamics of ongoing activity is complex and external stimuli elicit differentiated responses representing stimulus features. Under burst suppression, specific local responses are replaced by stereotypic bursts similar to those occurring spontaneously (Hudetz and Imas 2007). These evoked bursts spread to areas otherwise unresponsive to stimuli of the particular modality (Hudetz and Imas 2007; Land et al. 2012). Thus, on the one hand, the stereotypic patterns prevailing throughout the network hinder the emergence and propagation of “informative” intrinsic and stimulus-related neural activity. On the other hand, global connectivity increases, possibly leading to a reduction of the functional segregation of areas. Both aspects lead to a reduction in information processing capacity which may account for the loss of consciousness under anesthesia (Alkire et al. 2008).

Our results reinforce the hypothesis that envelope ICMs closely reflect patterns of structural connectivity between the underlying areas (Engel et al. 2013). We show that this effect is robust even to a severe perturbation such as that induced by anesthesia. For resting-fMRI, the high similarity of functional to structural connectivity has been shown by tracer injection (Vincent et al. 2007). Nonetheless, as suggested previously, ongoing activity has a role in shaping processing capacities of functional networks (Engel et al. 2013). Contrary to classical views, ongoing activity is likely to carry information and to be endowed with meaningful spatiotemporal structure, which reflects previous learning and can bias the processing of stimuli (Engel et al. 2001; Deco et al. 2011). While some aspects of this information are carried by envelope ICMs, others are represented in complementary connectivity measures such as phase ICMs (Engel et al. 2013), the s estimator (Carmeli et al. 2005, 2007) or transfer entropy measurements (see Wibral et al. 2013, Wollstadt et al. 2014). These aspects of ongoing connectivity are modulated on faster time scales than those reflected in envelope ICMs and thus may be responsible for adapting stimulus processing to rapidly changing contextual settings. Envelope ICMs seem to represent coherent excitability fluctuations that may lead to coordinated changes in the activation of brain areas. As we have proposed, this might regulate the availability of neuronal populations or regions for participation in an upcoming task and gate the emergence of more specific interactions for processing of stimulus- or response-related signals (Engel et al. 2013). However, these aspects clearly await further investigation.

Supplementary Material

Supplementary material is available at Cerebral Cortex online.

Funding

This work was supported by grants from the Deutsche Forschungsgemeinschaft (DFG): SFB 936/A1/A2, SPP 1665/EN533/13-1, and SPP 2041/EN533/15-1.

Notes

We thank Dorrit Bystron for excellent technical support. Conflict of Interest: None declared.

References

Author notes