Global and local complexity of intracranial EEG decreases during NREM sleep.

Key to understanding the neuronal basis of consciousness is the characterization of the neural signatures of changes in level of consciousness during sleep. Here we analysed three measures of dynamical complexity on spontaneous depth electrode recordings from 10 epilepsy patients during wakeful rest (WR) and different stages of sleep: (i) Lempel-Ziv complexity, which is derived from how compressible the data are; (ii) amplitude coalition entropy, which measures the variability over time of the set of channels active above a threshold; (iii) synchrony coalition entropy, which measures the variability over time of the set of synchronous channels. When computed across sets of channels that are broadly distributed across multiple brain regions, all three measures decreased substantially in all participants during early-night non-rapid eye movement (NREM) sleep. This decrease was partially reversed during late-night NREM sleep, while the measures scored similar to WR during rapid eye movement (REM) sleep. This global pattern was in almost all cases mirrored at the local level by groups of channels located in a single region. In testing for differences between regions, we found elevated signal complexity in the frontal lobe. These differences could not be attributed solely to changes in spectral power between conditions. Our results provide further evidence that the level of consciousness correlates with neural dynamical complexity.

Introduction

There is increasing evidence for a strong correlation between signal complexity of human electrophysiological activity and the level of consciousness. This is largely irrespective of the precise definition of complexity, which is variously given in terms of diversity and/or unpredictability of some signal feature across spatial regions and/or time in either spontaneous signals (Zhang ; Burioka ; Ferenets , 2007; Liu and Sun, 2015) or the response signals to perturbation via transcranial magnetic stimulation (TMS) (Sarà and Pistoia, 2010; Casali ; Sarasso ). The relevant studies consider levels of consciousness that differ either due to pharmacological action (Zhang ; Ferenets , 2007; Casali ; Sarasso ), chronic disorders of consciousness (Sarà and Pistoia, 2010; Casali ) or natural sleep states (Burioka ; Casali ; Abásolo ; Liu and Sun, 2015; Andrillon ). These observations support integrated information and complexity theories of consciousness that emphasize diversity of phenomenology (conscious experience) as a key property of consciousness to be reflected in its neural correlates (Tononi ; Tononi and Edelman, 1998; Seth ; Tononi, 2008).

In a recent study, we found evidence for a robust and spatially uniform decrease in three distinct flavours of spontaneous multi-dimensional electroencephalogram (EEG) complexity during propofol-induced general anaesthesia (Schartner ). These three distinct flavours were captured by: (i) a form of Lempel–Ziv complexity (LZc), which derives from the (lack of) compressibility of the data matrix; (ii) amplitude coalition entropy (ACE), which reflects the entropy over time of the constitution of the set of most active channels; (iii) synchrony coalition entropy (SCE), which reflects the entropy over time of the constitution of the set of synchronous channels.

In this study, we applied these measures to depth electrode recordings taken during non-rapid eye movement (NREM) sleep, rapid eye movement (REM) sleep and wakeful rest (WR) from 10 epilepsy patients undergoing pre-surgical evaluation. Complementing our analysis of spontaneous complexity of EEG under propofol, we here capitalize on this new dataset to (i) investigate complexity changes during sleep phases and (ii) examine regional versus global complexity changes, taking advantage of the high spatial resolution of depth electrodes compared to scalp EEG. Electrode locations varied across participants, and spanned cortical and sub-cortical regions. We tested the extent to which changes in the complexity measures were consistent and detectable irrespective of which regions were covered. This regional analysis of signal complexity is relevant for comparison with other signatures of the different sleep stages, some of which have been reported to differ strongly across cortical regions [e.g. slow-wave and sleep spindle propagation (Andrillon ; Nir )], while other signatures are exhibited more evenly across the cortex [e.g. average power spectra (Cavelli )].

A recent study of avalanche events during wake and sleep states has suggested that the diversity of spatiotemporal scales under which neural dynamics unfold is preserved during NREM sleep (Priesemann ). Given the similarity of our data set with the one studied in (Priesemann ), we additionally replicated their analysis on our data, and discuss the implications in conjunction with the observed behaviour of the complexity measures.

We found that all three complexity measures scored substantially lower in all participants during NREM than during WR, when computed across sets of channels broadly distributed across multiple cortical and sub-cortical regions. This global pattern was in almost all cases mirrored at the local level when analysing the measures across groups of channels located in just a single region. These differences persisted when controlling for changes in spectral power across changes in conscious level. We also found evidence for higher signal complexity in the frontal lobe than the other cortical lobes and, further, an increase in large avalanches during NREM.

Methods

Ethics statement and data protection

In agreement with the HORIZON 2020 requirement, the protocol used to collect the data analysed here has been drawn up in accordance with the EU standards of good clinical practice and with the Declaration of Helsinki (current revision) and is approved by the Ethics Committee of the Niguarda Hospital of Milan (protocol number: ID 939, Niguarda Hospital, Milan, Italy). All data related to the study participation are treated confidentially in compliance with good clinical practice as well as in compliance with Italian specific national laws on the protection of individuals. Participants are informed that personal data are collected and stored electronically, which can be used for purposes of scientific research, and that dissemination of the results can take place only in an anonymous and/or aggregate form. Participants are informed that they have the right to access the stored data, and to update or modify erroneous data.

Participants and data acquisition

The data were derived from a data set collected during the pre-surgical evaluation of 10 neurosurgical patients with a history of drug-resistant, focal epilepsy. All participants were candidates for surgical removal of the epileptogenic zone. The recordings were obtained from stereotactically implanted depth multi-lead electrodes (stereo-EEG, SEEG), inserted for the precise localization of the epileptogenic zone and connected areas (Cossu ). The investigated hemisphere, the duration of implantation, the location and number of recording sites were determined based on non-invasive clinical assessment.

SEEG activity was recorded from platinum–iridium semi-flexible multi-contact intracerebral electrodes, with a diameter of 0.8 mm, a contact length of 1.5 mm, an inter-contact distance of 2 mm and a maximum of 18 contacts per electrode [(Dixi Medical, Besancon France), see Fig. 1]. The individual placement of the electrodes was ascertained by post-implantation computerized tomographic imaging scans (post-CT), and Montreal Neurological Institute (MNI) coordinates obtained for each contact. Full details on the contact localization procedure can be found here (Arnulfo ). Briefly, post-CT was co-registered to pre-implant magnetic resonance imaging (MRI) by an algorithm based on rigid affine transformation and mutual information. Next, the post-CT scan was thresholded and skull-stripped in order to find and remove radiological artefacts. Given the planned entry points on the skull and electrode pin dimensions, the axis direction was estimated and each recording position iteratively computed. A region-of-interest around each contact was defined and a most probable anatomical label assigned, using the Destrieux atlas. Single participant channel positions were projected to MNI space by internal transformation files defined in the Matlab toolbox Freesurfer.

Figure 1

Depth electrode dimensions, location and channel definition. a) Outline of a multi-lead intracerebral electrode. b) Horizontal section of MRI scan coregistered to CT scan of participant 1, showing a multi-lead intracerebral electrode (visible inside the red rectangle) (Arnulfo ). c) In MNI space in millimetres, the location of initial contacts (blue dots) are partially overlaid by red crosses, marking the contacts chosen for bipolar-montage. Numbers index the resulting channels as used for analysis. Channels were obtained by applying bipolar-referencing to neighbouring contact pairs that were not discarded (see text for details).

In addition, scalp EEG activity was recorded from two platinum needle electrodes placed during surgery on the scalp at standard 10–20 positions Fz and Cz. Electroocular activity was recorded from the outer canthi of both eyes, and submental electromyographic activity was also recorded. Both EEG and SEEG signals were recorded using a 192-channel recording system (NIHON-KOHDEN NEUROFAX-110) with a sampling rate of 1000 Hz. Data were recorded and exported in EEG Nihon Kohden format. Recordings were referenced to a contact located entirely in the white matter.

Selection of recording contacts and data preprocessing

In each participant, recordings were made from up to 194 contacts (blue dots in Fig. 1c). For the present analysis, selection of recording sites was based on the following criteria: we excluded from the analysis those contacts that (i) were located in the epileptogenic zone (as confirmed by post-surgical assessment), (ii) were located over regions of documented alterations of the cortical tissue (e.g. Taylor dysplasia), as measured by the radiographic assessment, or (iii) exhibited spontaneous or evoked (Valentin ) epileptiform SEEG activity during wakefulness or NREM. Contacts located in white matter, assessed by MRI, were also excluded from analysis. Data samples were taken from each participant from four different states: WR, non-rapid eye movement sleep early at night (NREMe), non-rapid eye movement sleep late at night (NREMl) and REM. Sleep scoring was obtained as in Silber ), using one scalp EEG derivation together with one bipolar electrooculographic (EOG) and one electromyographic (EMG) derivation. WR recordings were taken at various times of day between 8 am and 6 pm, and the subjects were sitting on a bed with eyes closed. All NREM epochs were collected during stage N3, as defined in Silber ). NREMe corresponds to the first stable NREM (stage N3) episode and NREMl to the last stable NREM (stage N3) episode of the night (Pigorini ). By using only NREM in stage N3, possible fluctuation of the level of consciousness due to subliminal processing in NREM stages N1 and N2 (Andrillon ) were avoided. The data samples were imported from EEG Nihon Kohden format into Matlab and converted using a customized Matlab script.

Bipolar montages were calculated by subtracting the signals from adjacent contacts of the same-depth electrode (see Fig. 1a) to minimize common electrical noise and to maximize spatial resolution (Cash ; Gaillard ). To further minimize volume conduction artefacts, at most every third (bipolar) channel from each electrode was retained for analysis. The number of retained channels per participant varied between 18 and 31 (red crosses in Fig. 1c show an example choice of contact pairs; this is participant 1 in Fig. 2). For most analyses we used 18 channels per participant, selected as follows. A first electrode was chosen at random and the m channels on that electrode were all selected, ordered from the innermost outwards. The process was repeated until 18 channels had been selected (see Fig. 1c for an illustration of this ordering). Figure 2 depicts the channels’ anatomical locations.

Figure 2

Channel locations for the 10 participants. For each participant, coloured dots indicate the positions of channels that were used for analysis. The locations are plotted using MNI coordinates on a standard glass brain [python nilearn (Abraham )]; see text for details. Numbers in black blocks indicate participant number, coloured digits count channels per participant and region.

The data from the selected channels was further preprocessed as follows. No epochs were removed, as visual inspection did not result in the detection of severely artefacted epochs. The data samples were downsampled to 250 Hz and divided into 10 s segments. Linear de-trending, baseline subtraction and normalization by standard deviation was performed for each channel of each segment. After preprocessing, the length of the retained data sample for each participant and state varied between 7 and 16 min. Figure 3 illustrates representative channel activity for the different states.

Figure 3

Sample segments from five channels for each state. Segments are shown after pre-processing, prior to normalization by standard deviation. The recordings for NREMe display visibly stronger slow-waves (low-frequency components) as compared to those for REM or WR.

LZc

We compute LZc following our previous study (Schartner ). Briefly, for a given segment of data, LZc quantifies the number of distinct patterns of activity in the data, so that it is maximal for completely random data. It can be thought of as being proportional to the size of a computer file containing the data, after applying a compression algorithm. Computing the Lempel–Ziv compressibility of data requires a binarization of the multidimensional time series. Our threshold was based on the instantaneous amplitude of the Hilbert transform, i.e. the absolute value of the analytic signal of the channel’s time series, in order to capture the amplitude of the activity. The threshold T for the i channel was chosen as the mean of the absolute value of the analytic signal of the i channel. The data segment is then treated as a binary matrix, with rows corresponding to channels (time series) and columns corresponding to time (observations). A Lempel–Ziv compression algorithm obtains a list of words (binary subsequences that appear at least once) in the data matrix, as sketched in Fig. 4. The LZc is then proportional to the number of binary words. The greater the degree of randomness, the greater the number of different sub-sequences that will be present, and thus the higher the LZc.

Figure 4

Schematic of the computation of LZc and SCE. LZc: a) x is the activity of the i channel and a is the (Hilbert) amplitude of x. b) b is a binarized, using the mean activity of a as binarization threshold. c) The data after binarization of all n channels. d) The multidimensional time series is concatenated observation-by-observation into one binary sequence k, and then e) a list of distinct patterns is obtained and listed as a dictionary of binary words via a Lempel–Ziv algorithm. LZc is proportional to the size of this dictionary. SCE: a) Two time series. b) The analytic signals of these two, which are complex signals, with the real part being the original signal and the imaginary part being the Hilbert transform of the original signal. c) A binary synchrony time series is created for this pair of signals; 1 indicates that the phases of the complex values of the analytic signals are similar (absolute difference of less than 0.8 modulo .) d) Such time series are obtained to represent each channel’s synchrony with seed channel i. e) is the entropy over observations in the resulting data matrix . The overall SCE is then the mean value of across choices of seed channel i.

LZc is obtained by rearrangement of the binarized multidimensional time series, observation-by-observation, into a binary one-dimensional sequence as described in Fig. 4, and then applying a standard open-source Lempel–Ziv compression algorithm (Rosetta Code) to this sequence. LZc tells us how much variety there is in the patterns of activations over time.

We normalize LZc by dividing the raw value by the value obtained for the same binary input sequence (k in Fig. 4) randomly shuffled. LZc’s raw score for a binary sequence of fixed length is maximal if the sequence is entirely random. (The standard deviation of LZc for 50 different shufflings of the same input sequence was less than 0.002% of the mean value, thus the data matrices we analysed were sufficiently large such that there was negligible variance arising from basing the normalization on just a single random shuffling). Thus the normalized LZc values indicate the level of complexity on a scale of 0–1.

ACE

ACE is defined as in Schartner ) as the entropy (over time) of the constitution of the set (coalition) of channels that are ‘active’, given the binarization scheme described above for LZc for defining ‘active’/ ‘inactive’ channels. As for LZc we normalize ACE by its value for a random shuffling of the data (see section ‘LZc’). This measure is a variant of that first introduced in Shanahan (2010) and Wildie and Shanahan (2012); see Schartner ) for further discussion.

SCE

Synchrony coalition entropy (SCE) is also defined as in Schartner as the entropy over time of the constitution of the set of channels that are in synchrony (see Fig. 4 for a schematic). For data consisting of channels we consider two channels to be in synchrony at time t if the absolute value of the difference between their instantaneous Hilbert phases is less than 0.8 radians (≈45°). Then we define coalition time series by taking the value 1 if channels i and j are synchronized at time t and taking the value 0 otherwise. The coalition entropy of with respect to channel i is the entropy of (over time), normalized as a proportion of its maximum possible value N:

The overall SCE is then the mean value of the across channels. The upper bound SCE would arise from completely random coalition time-series in which each entry is 1 with probability 0.5. Such time-series are generated (with the same dimensions as those chosen for the data) to obtain the normalization factor N. Note that SCE does not score exactly 1 for shuffled input data—unlike ACE and LZc—as the probability at a given time point of two shuffled channels being in synchrony is less than 0.5.

Phase-randomized surrogate renormalization

In order to exclude that any changes in complexity merely reflect changes in spectral properties of the data, we also considered the complexity measures renormalized by their values for phase-randomized surrogate data. The procedure for carrying out the renormalization was previously described in Schartner ), and is as follows. From the complete data of a given participant in a given state, a segment was randomly chosen. Each channel was then expressed as a superposition of sinusoids using fast Fourier transform. Then the phase of each sinusoid was independently and randomly changed, before applying an inverse Fourier transform. The complexity measure was computed for 100 such phase-randomized data segments, and then the mean of these 100 scores was used to normalize the measure’s score for the original data segments. Averaging over 100 such data segments sufficed to obtain negligible randomness in the normalization factor. Phase-randomized surrogate normalized measures are denoted with a subscript .

Statistics

Analyses were performed using non-overlapping segments of length 10 s for a total length between 7 and 16 min of SEEG recording per participant and per state. The mean and standard error of the complexity measures’ scores were computed over these segments. At the single participant level, the effect size of differences between states was measured using Cohen’s d (Cohen, 1992). We call an effect size high if d > 0.8. For group level comparisons, a t-test was applied, with correction for false discovery rate (FDR) via the Benjamini–Hochberg procedure.

Results

Global analyses

We first computed the complexity measures across broadly distributed sets of channels. Specifically, for each participant we used the first 18 channels (out of the 18 to 31 available; see section ‘Methods’ for channel ordering), whose precise sets of locations varied across participants, but which spanned at least two cortical regions (see Fig. 2). Figure 5a shows the mean values across 10 s segments of LZc, ACE and SCE for each participant during WR, NREMe, NREMl and REM, normalized to scores for WR. Classification of effect sizes for differences between states at the single participant level are shown in Table 1. For all participants, the three measures ACE, SCE and LZc score higher for WR than NREMe with high effect size (Cohen’s d > 0.8). For all but one measure and participant, scores were also higher for REM than NREMe with high effect size. Values for NREMl typically lie in between those for WR or NREMe, with high effect sizes for most participants (see state pair NREMl/NREMe in Table 1). The exception is LZc for NREMl versus NREMe, for which the effect size is small for most participants. The difference in values between WR and REM is small for most cases (d < 0.8). Overall, signal complexity across all available brain regions is higher for REM and WR—states associated with consciousness—as compared to NREMe and NREMl.

Figure 5

LZc, SCE and ACE computed across broadly distributed sets of 18 channels. a) Individual participant plots. Plotted points show mean over 10 s segments. Scores are normalized by the score for WR (in addition to each measures’ normalization as specified in their definition). The measures score higher for WR as compared to NREMe or NREMl, consistently across participants. The measures score similar values for REM and WR. Error bars indicate standard error across segments. See Table 1 for effect sizes at the single participant level. b) Mean scores across the 10 participants. For each measure, the mean and standard error across participants are displayed, here not normalized relative to scores for WR. Significant differences between state pairs are shown by a solid line if P < 0.01 and a dotted line if P < 0.05 (t-test, FDR corrected for multiple comparisons). c) For comparison, the three complexity measures normalized by their values for phase-randomized surrogate data, demonstrating that complexity changes across conscious states are not due to changes in the spectral power profile alone.

Table 1

Effect size comparison per measure and state pair

	ACE	SCE	LZc
WR/NREMe	10, 0, 0	10, 0, 0	10, 0, 0
WR/NREMl	8, 2, 0	8, 2, 0	8, 2, 0
WR/REM	4, 5, 1	2, 8, 0	1, 9, 0
REM/NREMe	10, 0, 0	10, 0, 0	9, 1, 0
REM/NREMl	9, 1, 0	8, 2, 0	8, 2, 0
NREMl/NREMe	7, 2, 1	8, 2, 0	2, 8, 0

For each measure and state pair, the three numbers correspond to how many participants out of 10 had higher score for the left state with Cohen’s d  >  0.8 (left digit), no substantial difference, d  <  0.8 (middle), and higher score for the right state with Cohen’s d  >  0.8 (right), computed across segments. The results were obtained from applying the measures to 10 s segments from 18 channels as in Fig. 5.

In the Supplementary Material, the behaviour of the complexity measures is compared to that of normalized spectral power in the various frequency bands (δ, θ, α, β, γ), that of a simple correlation measure, sumCov, which equals the mean of the absolute values of the correlation coefficients between each pair of channels, and that of LZsum, the mean LZc of single channels (see Supplementary Table S1). As expected, normalized delta power is consistently substantially higher in NREMe than WR, while beta and gamma power are in almost all cases substantially lower. The high delta band power during NREM sleep reflects the presence of slow-waves (Alain et al., 1999; Massimini ; Nir ; Tom, 2013; Vyazovskiy and Harris, 2013; Buzsáki ).

It was important to verify that the observed changes in complexity during NREM sleep are not merely reflecting well-known changes in spectral properties of the data, in particular the increase in delta power. To address this we performed a surrogate data analysis, following our previous study (Schartner ). Specifically, we computed the measures normalized by their values for phase-randomized surrogate data, denoted by and , see section ‘Phase-randomized surrogate renormalization’. When repeating the above analyses with the phase-randomized surrogate normalization, the consistency of the measures’ behaviour was somewhat reduced, yet importantly at the group level all three measures still scored significantly higher for WR and REM than NREMe, see Fig. 5c. (See Supplementary Fig. S8 for plots showing the individual participants.) Thus, we conclude that the observed changes in complexity partly reflect spectral changes, but in particular the difference in complexity between WR/REM and NREMe goes significantly beyond that which can be accounted for by the difference in power spectrum.

For each state we further investigated the extent to which the complexity measures’ scores across subjects correlated with each other, as well as with normalized spectral power in the various frequency bands. Figure 6 indicates the pairs of measures and states that showed the strongest correlations (Pearson coefficient r of absolute value >0.7) across subjects. Strong positive correlations (r  >  0.7) were found for ACE/SCE for all states, and ACE/LZc for all states except REM. Weaker, yet still significant correlations were observed between SCE and LZc (see Supplementary Table S2). The complexity measures SCE and ACE also showed strong negative correlations () with delta power during NREM sleep states, but only weaker correlations during WR and REM, and r  >  0.7 was observed for SCE versus gamma power during NREMl. Weaker, yet still significant correlations were observed for LZc versus delta power, ACE versus gamma power and LZc versus gamma power (see Supplementary Table S2). The imperfect correlation between all the complexity measures indicates that they are capturing not entirely equivalent properties of the dynamics.

Figure 6

Thresholded correlations between complexity measures and frequency bands for each state. For a given state and participant the Pearson correlation between pairs of measures was computed across 10 s segments. For each state and measure pair (ACE, SCE, LZc and spectral power in the five canonical frequency bands, correlation of two spectral power bands not shown) a point is displayed if the absolute value of the Pearson correlation r, averaged across participants, is > 0.7. (These correlations are all significant at P < 0.05, corrected for false discovery rate.)

From one participant, data from NREM sleep stage N2 (NREM2) were available in addition to N3 (NREM3), obtained from a whole night recording. For each of the states REM, WR, NREM2 and NREM3, we concatenated all corresponding epochs of the whole night recording and obtained mean values for each complexity measure across 10 s segments, see Fig. 7. In line with a recent study demonstrating that LZc indexes different stages of NREM sleep when applied to scalp EEG (Andrillon ), we found for all three measures that WR scored higher than NREM2, which scored higher than NREM3. Perhaps surprisingly REM scored slightly higher than WR, although given that this did not extend across the whole set of 10 participants (Fig. 5), we do not discuss this further.

Figure 7

Complexity changes between NREM2 and NREM3. ACE, SCE and LZc calculated for one participant (number 2) for WR, REM, NREM2, NREM3. Epochs for each state throughout the whole night were concatenated and the measures computed for 10 s segments, all three scoring higher for NREM2 than NREM3. Solid lines indicate substantial (Cohen’s d  >  0.8) differences in scores between states.

Local analyses

Differences between states

Figure 5 shows that, when computed across widely distributed sets of channels, ACE, SCE and LZc take lower values during NREM than WR. In order to test if, or to what extent, there are regional differences in the decrease of complexity with NREM sleep, we carried out the following more local analyses.

First, we applied ACE, SCE and LZc to individual regions, following the classification in Fig. 2. For each participant, each region with four or more channels was analysed; where there were more than four channels, a quartet was picked at random. Depending on the participant-specific distribution of electrodes, there were two, three or four regions analysed per participant. The results (just for WR and NREMe) are shown in Fig. 8, and in general mirror the global result. The scores for LZc and ACE were in almost all cases greater for WR than for NREMe with large effect size. The consistency of SCE was weaker than for the other two measures (see Supplementary Table S3 for score counts). Importantly however, despite several instances of SCE showing increased values during NREMe, there was no region for which this was observed for more than one participant. We repeated this analysis using the phase-randomized surrogate normalization. With this alternative normalization, LZc still scored higher for WR than NREMe in almost all cases, although ACE and SCE no longer exhibited a consistent difference between these two states (see Supplementary Fig. S9). This further mirrors the global result that the observed change in complexity between these states is partly but not entirely reflecting spectral changes.

Figure 8

ACE, SCE and LZc scores computed for channel quartets restricted to single regions. The results for WR (white) and NREMe (red) were here obtained from applying the measures to quartets of channels within single regions, here specified by the first three letters of the label (as listed in Fig. 2). LZc and ACE score for all regions and participants higher for WR than NREMe (with just a single exception for ACE); SCE behaves less consistently here.

Secondly, in order to perform an analysis on the maximum possible number of regions for each participant, we computed the LZc of each single channel (LZs, same computation as LZc with trivial concatenation). Averaged across all channels per region and participant, LZs scored for 48 out of 50 region/participant pairs lower during NREMe than WR (see Supplementary Table S3 and Supplementary Fig. S2). Considering individual channels further, we noted also that normalized delta power is higher for NREMe than WR for 49 out of 50 region/participant pairs (compare Supplementary Fig. S2 with Fig. S3).

Thirdly, we investigated the complexity of the interaction of a channel in one region with a group of channels in another region. To this end we utilized the local SCE, , of a group of target channels in one region with respect to a seed channel i in another region (see section ‘Methods’ and Supplementary Fig. S4). When taking three target channels we observed lower values during NREMe than WR for 107 out of 129 choices of seed and target regions for all participants. Despite this inconsistent decrease of with NREM sleep, there were no exceptional choices of seed and target region that went against the trend for more than one participant (Supplementary Fig. S5). Similar results were obtained when taking two target channels (Supplementary Fig. S6), with 109 out of 143 choices of seed and target regions for all participants showing lower values during NREMe than WR, also without discernible pattern for the 34 exceptional cases with LZs higher for NREMe than WR. We conclude that there is no consistent local deviation from the global behaviour of SCE. Finally, we tested complexity of synchrony (CS) between pairs of channels, defined as the LZc of their synchrony time series (see Supplementary Material and Supplementary Fig. S7 for details). CS scored higher for WR than NREMe for 122 out of 139 pairs of regions.

In summary, all locally applied complexity measures behave predominantly according to the main trend, with just a few anomalies (Supplementary Table S3). There was no consistent pattern to anomalies across regions or participants.

Differences between regions for a given state

Differences in complexity scores between regions were overall less pronounced and less consistent than differences in complexity between states (see Supplementary Fig. 8). Here we investigate whether there are trends of regional complexity differences at the group level.

Mean scores for LZc, ACE and SCE were computed from channel quartets restricted to each of the four cortical lobes. All measures scored higher in the frontal lobe compared to other lobes for REM and WR, although an analysis of variance (ANOVA) test for a significant variation across the lobes gave a significant P value (< 0.05; P values uncorrected in this section) only for LZc when pooling across states (F = 4.1; P = 0.01). Given the available data, this analysis had , respectively for frontal, parietal, temporal and occipital cortex.

We repeated the analysis with LZs (single-channel LZc), computable whenever the participant has just one or more channels located in the given lobe. This led to the slightly larger sample sizes of , respectively for frontal, parietal, temporal and occipital cortex. For each state, LZs scored highest for the frontal lobe, see Fig. 9. A one-way ANOVA test for significant variation between lobes yielded for state WR, REM, NREMl, NREMe, respectively (F is the statistic, P the associated P-value from the F-distribution). A t-test as a posthoc analysis indicated a significant difference at P < 0.05 (uncorrected) between the frontal and temporal lobe in states WR and REM. When pooling values from all states for a given lobe, the ANOVA test did indicate a significant variation of LZs across lobes, with . By contrast, an ANOVA across lobes for the difference LZs(WR)-LZs(NREMe) yielded , confirming that the drop in complexity with NREMe does not differ substantially across lobes. For comparison, an ANOVA between lobes for delta power yielded , with a significant (t-test, P < 0.05, uncorrected) difference between the frontal and temporal lobe (delta higher in temporal) in states WR and REM; pooled across states we found .

Figure 9

Mean LZs scores for the four cortical lobes in each state. For each participant, LZs scores were averaged across all channels from the given lobe. The mean of these scores across participants is displayed with standard error across participants for each state. The labels on the X-axis give the first three letters of the lobe name (frontal, parietal, temporal, occipital), and for each lobe, there were data from 7, 9, 9 and 7 participants respectively. LZs scores higher for the frontal lobe than for the other lobes. ANOVA pooled across states for difference between regions gives F = 4.5 and P = 0.005; t-tests for differences between pairs of lobes came out non-significant (P  >  0.05) except for between Fro and Tem for states WR and REM (uncorrected for multiple comparison).

In summary, we found evidence for greater dynamical complexity in the frontal lobe compared to the other cortical lobes, and confirmed that complexity changes across states are more pronounced than across regions.

Complexity in different frequency bands

To test whether the relationship between complexity and conscious level is restricted to specific frequency bands, we re-analysed the LZc, SCE and ACE complexity measures on the data after frequency filtering (using Butterworth filters). The results for 18 channels per participant (chosen as above) are shown as average scores across participants in Fig. 10 for 6 different frequency bands and for 6 different high-pass cut-offs. Figure 10a shows that the elevated complexity of ACE, SCE and LZc in the WR state is most pronounced in the delta (1–4 Hz) and alpha (8–13 Hz) ranges, and the general trend is at least weakly present in all bands. Also in Fig. 10b, higher signal complexity in WR than NREMe is visible for all high-pass cut-offs, yet this difference becomes smaller with increasing high-pass cut-off.

Figure 10

ACE, SCE and LZc for frequency filtered data. a) Bandpass filtering (frequency ranges indicated in each column). b) High-pass filtering (frequency cut-offs indicated in each column). Scores were obtained from applying the measures to 10 s segments of frequency filtered data from 18 channels and then averaging across all segments from all participants. WR is in white and NREMe in red. Significant differences between state pairs are shown by a solid line if P < 0.01 and a dotted line if P < 0.05 (uncorrected t-test across participants). Error bars show standard error computed from the mean scores for each of the 10 participants.

In summary, higher signal complexity in WR than NREMe is present for almost all tested frequency bands and high-pass cut-offs, and it is strongest when frequencies smaller than 4 Hz are left in the signal. This shows that the decrease of signal complexity with NREMe is amplified by an increase of low-frequency waves, yet still exists after most of the applied frequency filters (see also Supplementary Fig. S8).

Avalanche statistics

Recent studies of avalanche events have shown that avalanche size distributions follow a power law to good approximation over a wide range of scales, during both wake and sleep states (Ribeiro ; Priesemann ). This maintenance of a diversity of scales in the neural dynamics during NREMe contrasts with the robust decrease in the three flavours of complexity that is the main observation of this article. Since the present data set is similar in structure to the one analysed in Priesemann ), we replicated from that study some analyses of the distribution of avalanche events during wake and sleep states.

We defined events and avalanches of events in the same way as in Priesemann ), in brief as follows. For each positive deflection between two zero crossings of a channel time series, the area under the deflection was calculated. An event was said to have occurred whenever this area exceeded a threshold. The threshold was set for each channel such that there was the same event rate of 13 events per second for all channels. The time series is then binned [we took the bin size to be half the mean inter-event interval as in (Priesemann )], assigning a 1 to a bin if an event occurs in it, and a 0 otherwise. An avalanche is defined as a cluster of events: in each time bin during an avalanche, there is an event occurring in at least one channel. Avalanches are preceded and followed by time bins in which there are no events. Three quantities that characterize avalanches are analysed: (i) avalanche size s, which is the total number of binary events that occur during the avalanche; (ii) avalanche duration d, which is the number of time bins spanned by the avalanche; (iii) the branching parameter σ, which is the average number of events in one time bin divided by the number of events in the preceding time bin, given that there was at least one event in the preceding time bin.

For the avalanche analyses, we used 18–31 channels per participant, selected and pre-processed as described in the ‘Methods’ section. We successfully replicated the findings of Priesemann ) study, that mean size and duration of avalanches, and the mean branching parameter for avalanche events were all greater during NREM sleep (we used NREMe only) than during WR or REM, see Fig. 11a. Further, for all states, the avalanche distribution exhibited no upper size bound, indicative of a maintenance of neural dyamical scale diversity across states (Fig. 11b). [Note however, that this result is possibly biased due to the small number of channels compared to other more detailed studies (Priesemann ; Touboul and Destexhe, 2010; Priesemann )]. We discuss these findings in conjunction with our findings for the complexity measures in the section ‘Discussion’.

Figure 11

Avalanche analysis as introduced by Priesemann , Figs 7 and 4, for the present data set. a) States shown are NREMe, REM and WR. Scores of each measure (mean avalanche size , mean avalanche duration , mean branching parameter across 10 participants) are shown for each state, normalized by mean score across all three states. Error bars indicate standard error across 10 participants. All three indices have substantially higher scores for NREMe than WR or REM (P < 0.01, FDR corrected t-test across participants). b) The avalanche size distribution for each state (not normalized) approximates a power law. Large avalanches are more frequent in NREMe than in WR or REM.

Controls

Given that our analyses depended on a random selection of channels, we checked that our results could be replicated for other selections of channels. We repeated the complexity measure analysis for different sets of 10 channels, each set chosen randomly across all available channels per participant, and found that for each of 10 different channel sets, the scores of all three measures are higher for WR than NREMe for all participants, substantially so (Cohen’s d > 0.8) for all but at most two participants (Supplementary Table S4). This channel-independent decrease in signal complexity for NREMe is in line with the single-channel Lempel–Ziv complexity (LZs) analysis, for which 48 out of a total of 50 channels across participants gave higher LZs scores for WR than NREMe.

We carried out a control test for dependence of the main result on sampling rate, computing ACE, SCE and LZc for data sampled at 10, 50, 150, 350, 500, 750 and 1000 Hz. The segment length and channel number was fixed at 2500 observations and 18 channels, respectively. All three measures scored for all participants higher for WR than NREMe, for a sampling rate of at least 150 Hz. Consistency of the measures was considerably weaker for the sampling rate of 50 Hz, and the measures became fully inconsistent across participants for 10 Hz, see Supplementary Fig. S11.

We also tested for effect of segment length (fixing the number of channels at 18 and the sampling rate at 250 Hz). We found that, for a length of 500 observations (2 s), all measures still scored for all participants higher in WR than NREMe, as was found above in the main analyses for 2500 observations (10 s). For 100, 50 and 30 observations per segment (0.4 s, 0.2 s and 0.12 s, respectively) there were at most two participants showing higher scores for NREMe than WR for each measure. At segments that were 15 observations long (0.06 s), inter-participant consistency was lost for all three measures, see Supplementary Fig. S12. Note that for fewer observations the influence of spatial complexity—as opposed to temporal complexity—on the score increases.

Finally, when omitting bipolar referencing as a pre-processing step, results are not substantially different (Supplementary Fig. S10). This shows that bipolar referencing as a pre-processing step of the data is not crucial for the detection of complexity changes in WR/NREMe.

In summary, our main results are robust to modification of sampling rate and segment length, provided one stays within reasonable limits, respectively to ensure the full physiological range of frequencies is sampled, and to maintain statistical power. This is in keeping with similar controls carried out in our previous study on scalp EEG data (Schartner ).

Discussion

We have analysed the behaviour of three measures of dynamical complexity on spontaneous depth electrode data recorded during WR and diverse sleep states: LZc, ACE and SCE. All three of these measures scored substantially lower (high effect size, d > 0.8) during NREMe sleep than during WR in all 10 participants. This is in spite of the recordings being taken from different sets of regions in each participant (as prescribed by clinical requirements for epileptic focus detection). For the majority of participants, LZc, ACE and SCE scores during NREMl sleep were in between those for NREMe and WR. By contrast, there was no overall significant difference in the scores of any of the measures between WR and REM sleep.

Our results complement the recent finding by Andrillon ) that LZc of scalp EEG decreases during NREM sleep, and obtains its lowest scores during periods of deep sleep when slow-waves are present. However, in contrast with our study, Andrillon et al. did report a small but significant decrease in LZc also during REM sleep compared to the waking state. This could have been due to the fact that in the Andrillon et al. study participants were engaged in a task during the waking state, whereas the participants in this study were sitting at rest with eyes closed.

Scores for the three complexity measures tended to be imperfectly positively correlated across participants, indicating that they are capturing similar yet not entirely equivalent signal changes. This is in line with model simulations in Schartner ), which show SCE behaving in some cases differently from ACE and LZc. Further, both ACE and SCE strongly correlated inversely with delta power. The complexity measures do however capture more than just spectral changes between states. This was explored in detail in Schartner ), and we confirmed this again on the present data set, via re-computation of the measures on surrogate phase-randomized data, see Fig. 5c and Supplementary Fig. S8. When analysing the measures on frequency restricted data, we found that the changes in complexity were most pronounced in the delta and alpha frequency bands.

The pattern of results obtained was in almost all cases preserved when the measures were computed across groups of channels restricted to a single cortical lobe (or specified sub-cortical region, see Fig. 8). Scores for NREMe were almost always lower than for WR, and exceptions to this followed no discernible anatomical pattern.

We also tested whether there were observable differences in local dynamical complexity between different cortical lobes. We found that on average the measures score higher for channels located in frontal cortex compared to parietal, temporal or occipital cortex, irrespective of the state. Given that there are substantial differences in anatomical structure between brain regions—e.g. notably structural connectivity tends to be much more dense within more anterior cortical regions (Modha and Singh, 2010; Van Den Heuvel and Sporns, 2011)—our results are suggestive of the denser connectivity of frontal cortex supporting increased dynamical complexity.

Significant differences in NREM sleep electrophysiology between cortical lobes have been observed by Andrillon ) and Nir ), namely inhomogeneities in the direction of propagation and frequencies of synchronous sleep spindles and slow-waves. However, when comparing different regions in terms of the magnitude of the difference in complexity between WR and NREMe, we found no evidence for differences between regions. Thus, any inhomogeneities in sleep spindle and slow-wave events are unlikely to have significantly affected complexity scores in our data. This is likely due to the complexity measures’ scores being averages across channels and segments, and thus being sensitive to steady state properties and not to transient events like spindles. Other studies on slow-waves did also find regional differences in slow-wave propagation, e.g. waves originating more frequently in frontal regions, (Sheroziya and Timofeev, 2014), but brain areas consistently recruited by the slow oscillation have not been identified (Dang-Vu ). This is in line with our observations on the homogeneity of delta power for single channels in the Supplementary Material (Supplementary Fig. S4) and a study by Cavelli ), showing that sleep-stage-dependent spectral coherence varied equally across the cortex.

Recent studies of avalanche events have shown that avalanche size distributions follow a power law to good approximation over a wide range of scales, during both wake and sleep states (Ribeiro ; Priesemann ), while large avalanches are more frequent in NREM than WR. We replicated the summary statistics from the study in Priesemann ), which utilized a data set similar to the present one. Specifically, we found mean avalanche size, duration and branching parameter all to be greater during NREMe sleep than WR or REM, and for all states no upper size bound for avalanches. The maintenance of a diversity of scales in the neural dynamics during NREMe sleep contrasts with the robust decrease in the three flavours of complexity that is the main observation of this article. The presence of larger avalanches during NREMe sleep contrasts with the finding that the EEG response to TMS is less widespread during NREM sleep than WR (Casali ). However, this holds for low-intensity TMS perturbations only; for higher intensity TMS the response signal spreads as far as during WR , but in a less complex, more stereotypical manner. A possible explanation is that widespread sub-cortical drive is responsible for the increase in large avalanches, while a simultaneous decrease in cortico-cortical effective connectivity is responsible for the drop in responsiveness to TMS stimulation. This would be in keeping with relay models of thalamo-cortical interaction (Coulon ; McCormick ). Further, if the sub-cortical drive results in more stereotypical cortical activity, then that could explain the decrease in dynamical complexity as measured by ACE, SCE and LZc. We are presently exploring such potential mechanisms with computer simulations.

In summary, we have found that three measures of dynamical complexity, capturing distinct aspects of signal diversity in space and time, all robustly decrease during NREM sleep, across local and global brain networks. Our application of these complexity measures represent new contributions to the statistical characterization of cortical signals for different sleep stages, providing well-defined signatures. Importantly, complementing decreased complexity of cortical response activity to TMS stimulation (Casali ), our measures provide direct evidence for a breakdown of differentiation between regions and diversity of brain states explored, during states of unconsciousness, as predicted by integrated information and complexity theories of consciousness (Tononi and Edelman, 1998; Seth ).

Data Availability

Data not publicly available for legal reasons.

Click here for additional data file.