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.
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 et al., 2001; Burioka et al., 2005; Ferenets et al., 2006, 2007; Liu and Sun, 2015) or the response signals to perturbation via transcranial magnetic stimulation (TMS) (Sarà and Pistoia, 2010; Casali et al., 2013; Sarasso et al., 2015). The relevant studies consider levels of consciousness that differ either due to pharmacological action (Zhang et al., 2001; Ferenets et al., 2006, 2007; Casali et al., 2013; Sarasso et al., 2015), chronic disorders of consciousness (Sarà and Pistoia, 2010; Casali et al., 2013) or natural sleep states (Burioka et al., 2005; Casali et al., 2013; Abásolo et al., 2015; Liu and Sun, 2015; Andrillon et al., 2016). 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 et al., 1994; Tononi and Edelman, 1998; Seth et al., 2006; 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 et al., 2015). 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 et al., 2011; Nir et al., 2011)], while other signatures are exhibited more evenly across the cortex [e.g. average power spectra (Cavelli et al., 2015)].
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 et al., 2013). Given the similarity of our data set with the one studied in (Priesemann et al., 2013), 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.
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 et al., 2005). 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 et al., 2015). 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.
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 et al., 2002) 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 et al. (2007), 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 et al. (2007). NREMe corresponds to the first stable NREM (stage N3) episode and NREMl to the last stable NREM (stage N3) episode of the night (Pigorini et al., 2015). 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 et al., 2016) 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 et al., 2009; Gaillard et al., 2009). 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.
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.
We compute LZc following our previous study (Schartner et al., 2015). 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 Ti for the ith channel was chosen as the mean of the absolute value of the analytic signal of the ith 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.
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 is defined as in Schartner et al. (2015) 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 et al. (2015) for further discussion.
Synchrony coalition entropy (SCE) is also defined as in Schartner et al., 2015 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 et al. (2015), 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 N.
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.
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.
|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|
|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 theMassimini et al., 2004; Nir et al., 2011; Tom, 2013; Vyazovskiy and Harris, 2013; Buzsáki et al., 2013).
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 et al., 2015). 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
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
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 et al., 2016), 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.
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
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
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
In summary, all locally applied complexity measures behave predominantly according to the main trend, with just a few anomalies (
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
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 .
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.
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
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 et al., 2010; Priesemann et al., 2013). 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 et al. (2013), 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 et al. (2013), 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 et al., 2013)], 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 et al’s (2013) 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 et al., 2009; Touboul and Destexhe, 2010; Priesemann et al., 2013)]. We discuss these findings in conjunction with our findings for the complexity measures in the section ‘Discussion’.
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 (
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
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
Finally, when omitting bipolar referencing as a pre-processing step, results are not substantially different (
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 et al., 2015).
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 et al. (2016) 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 et al. (2015), 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 et al. (2015), 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
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 et al. (2011) and Nir et al. (2011), 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 et al., 2008). This is in line with our observations on the homogeneity of delta power for single channels in theCavelli et al. (2015), 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 et al., 2010; Priesemann et al., 2013), while large avalanches are more frequent in NREM than WR. We replicated the summary statistics from the study in Priesemann et al. (2013), 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 et al., 2013). 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 et al., 2012; McCormick et al., 2015). 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 et al., 2013), 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 et al., 2006).
Data not publicly available for legal reasons.
We are grateful to the Dr. Mortimer and Theresa Sackler Foundation, which supports the work of the Sackler Centre for Consciousness Science. M.M.S. is funded by the Informatics Department of the University of Sussex. A.B.B. is funded by EPSRC [grant EP/L005131/1]. A.P., S.S. and M.M. are funded by EU [grant HBP-SP3-T32-2] ‘Wavescale’, Swiss National Science Foundation grant ‘Sinergia’ CRSII3_160803/1 and the James S. McDonnell Foundation Scholar Award 2013. Special thanks to Charlotte Rae, who helped with the anatomical classification and thanks to all the neurosurgeons of the ‘C. Munari’ Epilepsy Surgery Centre of the Niguarda Hospital. L.N. is supported by the Italian Ministry of Health [Targeted Research Grant RF-2010-2319316; L.N.].
Conflict of interest statement. None declared.