Multiscale Lempel–Ziv complexity for EEG measures

Clinical Neurophysiology 126 (2015) 541–548

Contents lists available at ScienceDirect

Clinical Neurophysiology journal homepage: www.elsevier.com/locate/clinph

Multiscale Lempel–Ziv complexity for EEG measures Antonio J. Ibáñez-Molina a,⇑, Sergio Iglesias-Parro a, María F. Soriano b, José I. Aznarte b a b

Department of Psychology, University of Jaén, Spain Hospital San Agustín, Linares, Jaén, Spain

See Editorial, pages 435–436

a r t i c l e

i n f o

Article history: Accepted 8 July 2014 Available online 18 July 2014 Keywords: EEG Lempel–Ziv Multiscale method Complexity

h i g h l i g h t s  We showed that the classical Lempel–Ziv complexity measure neglects rapid rhythms of the EEGs.  We proposed a method in which the median based binarization previous to LZC calculation, was

replaced with multiple binarizations at different scales.  Our results indicated that the new approach better characterized the complexity of simulated signals

and EEGs from eyes-closed and eyes-open conditions.

a b s t r a c t Objective: To demonstrate that the classical calculation of Lempel–Ziv complexity (LZC) has an important limitation when applied to EEGs with rapid rhythms, and to propose a multiscale approach that overcomes this limitation. Methods: We have evaluated, both with simulated and real EEGs, whether LZC calculation neglects functional characteristics of rapid EEG rhythms. In addition, we have proposed a procedure to obtain multiple binarization sequences that yield a spectrum of LZC, and we have explored whether complexity would be better captured using this computation. Results: In our simulated signals, classical LZC did not capture modulations of a rapid component when a slower component of more amplitude was included in the signal. In real EEGs from healthy participants with eyes closed and eyes open, classical LZC calculation failed to show any difference between these two conditions. However, a multiscale LZC showed that complexity was lower for eyes closed than for eyes open conditions. Conclusions: As hypothesized, our new approximation captures the complexity of series with fast components masked by slower rhythms. Significance: The method we introduce significantly improves LZC calculation, and it allows a better characterization of complexity of EEG signals. Ó 2014 International Federation of Clinical Neurophysiology. Published by Elsevier Ireland Ltd. All rights reserved.

1. Introduction Complexity is one of the more explored properties of the EEG, especially when it is considered the result of a nonlinear dynamical ⇑ Corresponding author. Address: Departamento de Psicología, Universidad de Jaén, Paraje las Lagunillas s/n, 23009 Jaén, Spain. Tel.: +34 953211988; fax: +34 953211881. E-mail address: [email protected] (A.J. Ibáñez-Molina).

system. Although there is not agreement about the conceptual meaning of complexity, it has been related with terms as irregularity, desynchrony or randomness (see Costa et al., 2005; Stam, 2005, for a review). In the study of the EEG’s structure, one of the most relevant conceptualizations of complexity was given by Kolmogorov (1965) and Solomonoff (1964). According with these authors, the complexity of a sequence is the number of bits of the shortest computer program which can generate this sequence. Lempel and Ziv (1976) proposed an estimation of this complexity

http://dx.doi.org/10.1016/j.clinph.2014.07.012 1388-2457/Ó 2014 International Federation of Clinical Neurophysiology. Published by Elsevier Ireland Ltd. All rights reserved.

542

A.J. Ibáñez-Molina et al. / Clinical Neurophysiology 126 (2015) 541–548

that has been named Lempel–Ziv complexity (LZC), and Kaspar and Schuster (1987) provided an easy algorithm to calculate it. When this algorithm is used with EEG data sets, the signal is binarized by means of a threshold (Td) which is usually the median or the mean of the voltage values of long EEG segments. The resulting binary string is analyzed with the LZC algorithm. This process scans the binary sequence and parses it into subsequences with different patterns. A complexity counter C(n) increases by 1 unit when a new subsequence appears during the process and the immediately next symbol is regarded as the beginning of the next subsequence. The scan continues until the end of the sequence. For example, C(n) equals 6 for the sequence 11001010001111 because different patterns are 11001010001111 (a detailed explanation of the algorithm can be found in Kaspar and Schuster 1987; Zhang et al., 2001). LZC is important in neuroscience because it has been widely used to characterize the EEG of several mental and neurological disorders (Nagarajan, 2002; Aboy et al., 2006). It has been shown that schizophrenics and patients with depressive disorders exhibit a high LZC (Li et al., 2008). The level of complexity changes with patients’ age (Fernández et al., 2011, 2012). Alzheimer patients show an abnormally low LZC (e.g., Hornero et al., 2009; Dauwels et al., 2011). It has been demonstrated that LZC is a reliable measure for detection of early seizure onsets (Jouny and Bergey, 2012). LZC has also been used to detect different levels of consciousness. For example, Zhang et al. (2001) calculated LZC to detect levels of consciousness during deep anesthesia. The main reason why LZC is frequently used in the EEG analysis is because it has several advantages: (a) LZC can be applied to any type of time series. (b) LZC can be used with very short time series. (c) LZC can be used with non-stationary signals. (d) LZC is a nonlinear measure that is non-parametric and very easy to compute. However, from our point of view, there is one limitation in LZC computation when applied to EEG signals. This limitation will be the main focus of our study. As it was mentioned above, LZC values are indirectly calculated from a binary sequence obtained from the original signal. The binary sequence is constructed using a Td. Values greater than Td are assigned a one and values equal or lower than Td are assigned a zero in the binary sequence. Although the binarization of the signal simplifies LZC calculation, it also produces a loss of information from the original signal, and this information could have a great value to better characterize brain functions. This is precisely the point we would like to emphasize in this study. As it is widely known from the beginning of the study of the EEG, the amplitude of neural rhythms that give rise to EEGs is inversely related with its frequency (Niedermeyer, 1999). This relationship is quasi-linear when we consider the log–log of power and frequency of the signals (e.g., Buzsáki and Draguhn, 2004). In general, slow rhythms have greater amplitude than faster rhythms (see Sanei and Chambers, 2007, for a complete review on EEG signals). This dependence between amplitude and frequency is important in LZC calculation because binarization can be guided by slow rhythms; that is, variations of a slow EEG rhythm (high amplitude) will have more influence in the changes of its median or mean than faster rhythms (low amplitude). While high amplitudes will deviate the signal values above and below the median, low amplitudes will frequently oscillate far from the median. Hence, these far-from-the-median oscillations will not change the binary sequence (see Fig. 1 for a graphical explanation). In this study, we have demonstrated that binary sequences obtained from the EEG neglect functional characteristics of rapid EEG rhythms. In addition, we propose a procedure to obtain multiple binarization sequences that yield a spectrum of LZC. We suggest constructing multiple binarization sequences with different smoothed versions of the original EEGs as Tds. In our proposal,

binary sequences would be obtained by comparing each data point of the EEG with its smoothed versions. If the original series are smoothed using a long window and these values are used as Tds, the binary sequence will characterize oscillations of high amplitude that are present in the original signal (similarly to when the median is the Td). On the contrary, a short windowed smoothing will produce binary sequences characterizing fast rhythms with low amplitude (see Fig. 1 for an example). As each length of window will be sensitive to a range of amplitude-frequency, we believe complexity is better captured using multiple window lengths and then, multiple LZC values. Hence, the method we introduce does not replace but extends the classical LZC computation. Thus, the main goals of our study are (1) to show how the binarization process in classical (median based) LZC calculation neglects fast EEG rhythms and (2) to introduce a multiscale binarization obtained from a smoothing procedure with multiple window lengths. In order to achieve our goals we compared the classical median based LZC with a new method on series of simulated signals and real EEGs. 2. Signal simulations We created signals with fast and slow components. Our hypothesis was that when the amplitude of slow components was much larger than the one of fast components, the median based binarization will be guided by the slow component. This type of signal shares a basic aspect with the EEGs which is that slow rhythms with high amplitude as delta or even alpha would guide binarization more than fast rhythms with lower amplitude as beta or gamma. We hypothesized that the LZC from multiscale binarizations would be much more sensitive to fast components than median based binarizations, and that this would be more evident when the Td comes from short windowed smoothed signals. As we stated in the previous section, it is clear that when a short window is used in a smoothing procedure, the resulting signal is a low pass filtered version of the original. If this signal is used as a Td, the binary sequence would reflect variations in the fast component and then, the LZC would be its characterization. 2.1. Methods 2.1.1. Characteristics of simulated signals A total of 100 signals of length N = 2.2001  104 were created. Each one was composed of 2 components and Gaussian white noise. One component of the signals was a band-pass filtered noise centered in 8 Hz (slow rhythm), and another component was a band-pass filtered noise centered around 30 Hz (rapid rhythm). The third component was created with Gaussian distributed noise. Each time point of the signal represented one millisecond (sampling rate of 1000 Hz), and each unit of amplitude represented one microvolt (lV), the unit in which the EEG signals are quantified. In order to control amplitude differences between components, each one was scaled as explained below. In order to obtain signal variability for further statistics, we made 40 independent realizations of each simulated signal. Formally, each simulated signal {x(n)} was the sum of three components:

xðnÞ ¼ A Noise1 þ B Noise2 þ C Noise3 ;

n ¼ 1; . . . ; N

ð1Þ

where Noise1 is a 5–11 Hz band-pass filtered Gaussian white noise (l = 0; r = 2), Noise2 is a 27–33 Hz band-pass filtered Gaussian white noise, and Noise3 is unfiltered Gaussian white noise. All sources of noise were selected using probability distributions of

A.J. Ibáñez-Molina et al. / Clinical Neurophysiology 126 (2015) 541–548

543

Fig. 1. Simulated signal with two components and noise; one component of 8 Hz with high amplitude and another component of 30 Hz with low amplitude. In order to compare binary sequences obtained from the median and a smoothed signal, Td31 and TdN are plotted with the signal. Both sequences are shown at the top part of the figure. See how the median based binarization (PN(n)) tends to be 1 on the left part of the sequence and 0 on the right part. It responds to the 8 Hz oscillatory component. The binary sequence constructed with Td31 (P31(n)) showed more frequent alternation between series of 1s and 0s which corresponds to the faster rhythm (30 Hz).

l = 0 and r = 2. Filters were applied using a two-way least squares FIR filter as it is implemented in the EEGLAB toolbox (Delorme and Makeig, 2004). A, B, and C quantify amplitude variations of components. A and B were calculated by multiplying a random normally distributed number (rnd) from the interval (0.8, 1) by each fixed amplitude. Values of A were calculated by multiplying rnd by 35, 40, 45 or 50; and B by multiplying rnd by 5, 10, 15, 20 or 25. Values of C were 5, 10, 15, 20 or 25 (see Fig. 2b to see the mean power density spectra of signals). 2.1.2. Multiple binarization with smoothed versions of the signals As mentioned above, a 0–1 sequence {P(n)} = s(1), s(2), . . .s(N), was created by comparison of each data point x(n) in the series with its Td in the following way:

SðnÞ ¼



0

if xðnÞ < T d

1 if xðnÞ P T d

ð2Þ

Twenty-one binary sequences were created for each signal. The first binary sequence was constructed using the median of the entire signal as Td (TdN). The other 20 binarizations were created using smoothed versions of the signal as Tds. Each data point x(n) had a unique Tdw(n) which was calculated by:

     wk  1 wk  1 ; . . . ; xðnÞ; . . . ; x n þ ; T dw ðnÞ ¼ median x n  2 2 wk  1 wk  1 n¼1þ ;...;N  ð3Þ 2 2 where W = [wk,. . ., wm], k = 1,. . ., m is the vector that contains window lengths of the smoothing procedure. Note that according to Eq. (3) a total of wk21 points from the sides of the series are not included in the computation of Tdw(n). For simplicity, all binary sequences were made of the same length (N  wmax  1). In our study, W = [201, 191,. . ., 21, 11]. This wide range is adequate to study a large variety of oscillatory components. In the case of the present simulation, since we know beforehand the oscillatory properties of the components in the signal, we will pay more attention to some window lengths of interest .They can be explored according to the wavelengths of the slow and fast components introduced in the simulations. Specifically, if we divide the sampling rate by a given frequency we will obtain the number of points of one oscillation (wavelength). For example, a window length that corresponds

to 30 Hz would be 1000/30  33.3. Then, in our selected windows, lengths around 31 and 41 will capture changes in the 30 Hz component. In general, a slightly larger window than the fast component and sensibly shorter than the one of the slow component will produce a signal without the fast component (low pass filter). The window has to be a bit larger than the oscillatory component of interest due to small random variations in the centered frequency of the component. Importantly, if the smoothed signal obtained with this window is used as a Td, then, any difference between the original and the new signal will be due to the rapid rhythms of about 30 Hz contained in the original. Hence, if we set the low pass filtered version of the signal as the Td, the resulting binary sequence will reflect rapid variations of the signal. On the contrary, a window length larger than the slow component will produce a nearly flat signal with values much closer to the median. In our simulations, 201 was a window length longer enough to capture 8 Hz oscillations (it will neglect oscillations slower than 5 Hz because 1000/ 5 = 200). Then, this window length was of interest in the study as a reference for comparison with the shorter ones. For each {x(n)}, twenty {Pw(n)} were created, where the sub index w denotes the window length applied in the smoothing. We estimated LZC from P201(n) to P11(n) where wk+1 = wk  10. We also estimated PN(n) which is the binary sequence generated from the median. As window lengths smaller than 101 are not influenced by the 30 Hz rhythm (it will capture frequencies close to 10 Hz), we expected that at least LZC calculated from P101(n) to P11(n) will be more sensitive to fast complexity fluctuations than LZC from PN(n) and larger windows. Another possible outcome in our work would be an effect of the amplitude of the noise added to the signals. However, we did not expect a significant contribution of noise in the effects when its amplitude or strength was similar to the mean amplitude of the fast component. 2.1.3. Estimation of multiscale LZC In order to obtain the LZC spectrum of all {x(n)}, each Pw(n) was explored according to the following steps: (a) We designed a moving window procedure with a 2  103 data points of window length and 90% of overlapping. Even though in this study we do not explore time changes in the signals, a sliding window procedure is adequate as a potential method to obtain a good time resolution and thus

544

A.J. Ibáñez-Molina et al. / Clinical Neurophysiology 126 (2015) 541–548

Fig. 2. Multiscale LZC means (a) with standard error bars, and power plots (b) for simulated signals of different amplitudes of the fast component. The shadowed area in the graph indicates statistically reliable windows (equal and smaller than w141). Differences in the fast component reached a maximum in the range of w51 to w31 (approximately the wavelengths of a 30 Hz rhythm). In order to help comparisons between power plots and LZC plots, Window lengths decreased from left to right in the x axis.

it can be used as a tool to detect dynamical changes. The total amount of scanned points on each Pw(n) was 2.18  104. The total number of LZC computations for each signal was 91. (b) LZC was calculated for each window by means of a complexity counter Cw(n). During a left to right scan of a given binary sequence, Cw(n) increased by one unit every time a new subsequence of consecutive characters was encountered. (c) Each LZCw was obtained when Cw(n) values were normalized with

C w ðnÞ LZCw ¼ n=log2 n

ð4Þ

where the sub index w indicates the window length of the smoothing that produced Pw(n). Note that LZCN will refer to the median based LZC.

(d) The final LZCw value of each signal was calculated by the average of all 91 values obtained with the moving window procedure. (e) In order to study the possible effects of moving window length, we repeated steps (a)–(d) for window lengths of 5  103, 1  104 and finally, we calculated the LZC spectrum for a single window of 2  105 points. 2.2. Results of simulations In order to test the reliability of the differences in complexity across window lengths (see Fig. 2a), we designed a 5  21 ANOVA with factors: amplitude of fast components (a5, a10, a15, a20 or a25) and window length (wN or w201,.., w11). Note that wN will produce LZCN. In order to introduce variability in the LZCw values, we conducted 40 realizations of simulations. Then, complexity was evaluated according to the amplitude of the fast component, and

A.J. Ibáñez-Molina et al. / Clinical Neurophysiology 126 (2015) 541–548 Table 1 ANOVA results depending on window lengths (wN, w251, w241. . . w21, w11), amplitude of fast components (fast components in the Table) (a5, a10, a15, a20, a25) and their interaction (interaction) with level of noise with 4 scaling factors (5–10, 10– 15, 15–20 and 20–25). Overall noise

F

d.f.

p

g2p

Fast components Window lengths Interaction

79.663 28401.60 110.070

4, 156 20, 780 80, 3120

<.0001 <.0001 <.0001

.671 .999 .738

Noise 5–10 Fast components Window lengths Interaction

39.935 16331.19 56.578

4, 156 20, 780 80, 3120

<.0001 <.0001 <.0001

.506 .998 .592

Noise 10–15 Fast components Window lengths Interaction

25.290 6923.06 36.805

4, 156 20, 780 80, 3120

<.0001 <.0001 <.0001

.393 .994 .486

Noise 15–20 Fast components Window lengths Interaction

31.461 8108.70 32.349

4, 156 20, 780 80, 3120

<.0001 <.0001 <.0001

.447 .995 .453

Noise 20–25 Fast components Window lengths Interaction

18.911 8436.93 45.642

4, 156 20, 780 80, 3120

<.0001 <.0001 <.0001

.327 .995 .539

amplitudes of the slow component were collapsed in the analyses. Results of this ANOVA are summarized in the Table 1. A main effect of amplitude of fast components was found showing that LZCw was different depending on the amplitude of the fast component in the signal. A main effect of window length was also found indicating that LZC was higher as the window length shortened (see Fig. 2a). More importantly here, the amplitude of fast components and window length factors interacted. To qualify this interaction, we conducted separated ANOVAs for each window length. p-Values from multiple comparisons were corrected using the Bonferroni–Holm familywise error rate controlling procedure (Holm, 1979). These analyses showed significant effects of amplitude of fast components in windows equal and shorter than 121 (the shadowed area of Fig. 2a indicates significant effects of fast at a given window length). As can be seen in Fig. 2a, binarizations with window lengths of 31 or 41 points easily captured variations of the 30 Hz component. Longer windows (201–141 points of length) and LZCN did not yield to different complexities with rapid rhythm modulations. Results for moving windows of length 5  103, 1  104 and 2  104 data points yielded to equivalent effects, indicating that this factor was not relevant in our simulated time series (only components slower than 0.5 Hz would have made a difference). Moreover, the time of computation was much higher for longer segments. We also wanted to explore the effect of noise in the LZCw calculation. With this aim, we compared the effect of the fast component at each window length for four levels of Gaussian noise. Each condition of noise was the result of two scaling factors (C in Eq. (1)). All signals were divided in four groups with the following scaling factors of noise: 5 and 10; 10 and 15; 15 and 20; 20 and 25. Then, we obtained the mean values of LZCN and LZCw for the fast component in each group. As can be seen in Fig. 3, the general pattern of results was not different for each level of noise (see also Table 1 for results of ANOVAs on each level of noise). Significant differences between the fast components were found at similar window lengths (see the shadowed areas of Fig. 3 for significant effects). It is noteworthy that as C increased, LZCw from long windows and LZCN, increased as well. This means that the background level of noise in the signals modulated the range of LZC values in the study.

545

These analyses demonstrate that LZCN alone may not capture complex properties of a fast rhythm when another slow rhythm is also present and the amplitude of the latter is much higher. In the next section we will exhibit results from multiscale LZCw tested with real EEGs. 3. EEG experiment It is known that there is a general increase in the power of the alpha rhythm when participants are relaxed with EC relative to EO. In addition, the remaining rhythms are also augmented when EC, but not present across the whole scalp (Barry et al., 2007). This is probably due to a disruption of information in the visual input when patients are resting with EC (Barry et al., 2007, 2009). We believe that these EC vs. EO changes in the power density spectra could be a good context to test the proposed multiscale LZC. First, power changes for EC and EO are greater and distributed across the scalp in the alpha band (8–12 Hz). This effect might be much more evident with window lengths sensitive to this wavelength, while complexity values obtained with the classic LZC for EC and EO could be similar if there are minor differences in slower rhythms. Second, if, in general, all frequency bands diminish when EO, it may reflect a general desynchronized oscillatory pattern that could be captured at different scales using our method. Third, it is not clear from previous studies with LZC of EEG data if LZCN will be different for EC and EO (see however Watanabe et al., 2003 for a multivariate approach). 3.1. Method 3.1.1. EEG recording and preprocessing Two sets of 100 EEG segments were selected from a public database provided by Andrzejak et al. (2001). One set was registered when participants rested with eyes closed and the other set when they were resting with eyes open. All EEG segments were artifact free recordings from five healthy participants at electrodes FP1, FP2, F7, F3, Fz, F4, F8, T3, C3, Cz, C4, T4, T5, P3, Pz, P4, T6, O1 and O2 (10–20 system). The length of each segment was 23.6 s digitized at a sampling rate of 173.61 Hz and band-pass filtered at 0.53–40 Hz. 3.1.2. EEG binarizations and LZC calculation Each EEG segment was binarized with thresholds corresponding to the median and 124 smoothing versions with window lengths: W = [251, 249,. . ., 7, 5]. As mentioned in the previous section, these window lengths can be related with different frequency bands by dividing the EEG sampling rate with the approximate band of interest. For example, delta oscillations of 2 Hz (173.61/2  86.8), theta of 5 Hz (173.61/5  34.7), alpha of 10 Hz (173.61/10  17.3), beta-1 of 19 Hz (173.61/19  9.1) and beta-2 of 35 Hz (173.61/35  4.9). As mentioned somewhere above, as each window length will be sensitive to similar and faster rhythms, it will be appropriate to choose a slightly larger window than the wavelength of the frequency of interest. As in the simulations, LZCw measures were obtained for each sequence from the means of LZCw series. These series were calculated with a moving window procedure in which segments of 2 s (348 samples) were analyzed with a 90% overlap. This length was selected taking into account that it is sensitive to rhythms down to 0.5 Hz, which is very close to the high-pass filter at 0.53 Hz applied to the EEGs. Although, larger segments in the LZC computation would not capture further oscillatory patterns, an effect of segment size was still possible. Gómez et al. (2006) analyzed MEG signals from controls and patients with Alzheimer disease

546

A.J. Ibáñez-Molina et al. / Clinical Neurophysiology 126 (2015) 541–548

Fig. 3. The effect of Gaussian white noise in the LZCw means. Four levels of noise were explored. The effects of the fast components were similar for all levels. It can be observed that high levels of noise yielded to higher levels of LZC at longer windows.

and found that classical LZC estimations stabilized at segment sizes of about 2000 sample points. Given that these authors used a similar sampling rate as the one in this study, we repeated the LZCw estimations at 20 segment sizes with samples ranging from 348 to 2001 with increasing steps of 87 points (500 ms). Analyses of variance on this data set demonstrated that LZCw differences for EO and EC conditions were not different across all segment sizes. Hence, for simplicity, in the next section we report analyses only for segments of 2 s. 3.2. Results of the EEG experiment We compared the LZCN and 124 different LZCws for W = [251, 249,. . ., 7, 5] during EC and EO conditions. We conducted a within participants ANOVA design with two factors: state (EC or EO) and window length (wN, w251, w249,. . ., w7, w5). The results of this ANOVA are summarized in Table 2. A main effect of state indicated that LZCw values were lower for EC than EO (see Fig. 4a). The main effect of window length showed that LZC increased as the window length decreased. This inverse trend between window length and LZ values was qualified by a state  window length interaction showing that the inverse trend only holds for the EO condition. To better explore this interaction, we conducted separated planned comparisons of EC vs EO on each window length. As in the simulations, p-values were corrected using the Bonferroni–Holm familywise error rate controlling procedure (Holm, 1979). Windows equal and smaller than 113 clearly exhibited lower LZC for EC when compared with EO (the shadowed area in Fig. 4a indicates significant differences between EC and EO). However, there were Table 2 ANOVA results depending on window lengths (wN, w251, w249. . . w7, w5), experimental condition (state: EC vs. EO) and their interaction (interaction). Overall noise

F

d.f.

p

g2p

Window lengths State Interaction

12.865 29.042 150.840

1, 198 124, 24,552 124, 24,552

<.0001 <.0001 <.0001

.06 .598 .432

no reliable differences between the complexity of EC and EO when the window length was more than 113 points or N. These results indicated that LZCN itself did not capture complexity differences between the EEG of EC and EO conditions in healthy participants. When the binarization was obtained with the median as Td, LZC was equal for EC and EO conditions. Complexity differences appeared only for multiscale LZCw estimations when windows lengths were smaller than 113. These differences strongly suggest that the EC condition produced a less complex signal than the EO condition. In Fig. 4b, the power plot of EC and EO conditions can be compared with LZCw. As it can be seen, the power increase peaking at alpha corresponded to a decrease in complexity that reached a minimum around w = 19 and w = 21. These window lengths are surprisingly close to the window lengths of an 8–10 Hz frequency rhythm. 4. Discussion 4.1. LZCN ignores fast components in the EEG signals In this paper we present a method that improves LZC estimations for EEG series. We demonstrated that LZCN may neglect rapid rhythms of EEGs. As indicated above, it is known that EEG’s frequency bands and power show an inverse linear progression in a natural logarithmic scale. Hence, any EEG binarization using the median as a Td will be led by slow components of the signal (high amplitude). We have demonstrated this by means of a series of signal simulations. LZCN did not capture modulations of a rapid fast component when a slower component of more amplitude was included in the signal. Aboy et al. (2006) made an important contribution in the understanding of LZC calculated to biomedical signals. They claimed that LZCN was not modulated by the number of harmonics of signals, and that amplitude modulations did not result in a different LZCN. However, they did not explore different interactions between frequency and amplitude. This is an important issue in the case of EEGs because, as shown in Fig. 1, amplitudes of components can interact with their frequencies in the

A.J. Ibáñez-Molina et al. / Clinical Neurophysiology 126 (2015) 541–548

547

Fig. 4. Multiscale LZC means (a) with standard error bars, and power plots (b) of the EEGs from healthy participants at EC and EO conditions. The scaling procedure revealed less complex signals in the case of EC vs EO. The shadowed area in the graph indicates statistically significant differences. In order to help comparisons between power plots and LZC plots, Window lengths decreased from left to right in the x axis.

process of binarization. We have shown these dependencies when we found that rapid fluctuations of low amplitude were not accounted by binary sequences constructed with the median. LZCN values did not change when the amplitude of the fast component was modulated.

4.2. Multiscale LZC extends complexity estimations to fast rhythms To improve the accuracy of LZCN for EEG, we proposed a multiscale method in which multiple binary sequences are created from a single signal and then scanned with the LZC algorithm. To test this method, LZCws were calculated for simulated signals and real EEGs. Our results suggested that the use of smoothed versions of the signal as Td for binarization easily captures the complexity of series with fast components masked by slower rhythms. As shown in Fig. 2, short LZCw differed significantly when a frequency component centered in 30 Hz changed its amplitude, while long LZCw and LZCN failed to characterize this modulation. These results demonstrate that scaling in binarization may be a necessary exten-

sion of the LZC computation when we are interested in all oscillatory spectra of a given EEG. Another interesting result obtained in our study is that as the window length of the smoothing procedure decreased, the complexity of resulting binary sequences increased. This finding agrees with the results obtained by Aboy et al. They found that an increase in the frequency rate of signals yielded to more complex signals. Similarly, short windows in our smoothing procedure captured rapid rhythms or frequencies, and then, an increase in complexity of these binary sequences was also expected (see Figs. 2 and 3). This effect was modulated when the amplitude of Gaussian noise was increased, and this is also similar to the results provided by Aboy et al. where an increase in this type of noise produced LZC estimations that approached to 1. One relevant question here might be if a scaling procedure is, in general, an adequate methodological strategy in the case of nonlinear measures of complexity. To answer this, we will mention that a similar scaling procedure has been developed in the case of the sample entropy calculation (see Pincus, 2006 for details of this measure). Costa et al. (2002, 2005) developed a multiscale

548

A.J. Ibáñez-Molina et al. / Clinical Neurophysiology 126 (2015) 541–548

entropy analysis (MSE) to obtain multiple estimations of the sample entropy from time series. From original series, they constructed new series by averaging the data points within windows of increasing length. Then, the sample entropy was calculated for each new signal. This procedure has been applied to EEG signals showing interesting interactions at different scales (remarkable examples were provided by Takahashi et al., 2009, 2010; Catarino et al., 2011; Liang et al., 2014). Despite the similarities between our method and the one proposed by Costa et al., it is noteworthy that in the case of the multiscale LZC, the scaling is aimed to capture different oscillatory properties while in MSE, the scaling is applied to explore temporal dependencies in the signal. In general, MSE studies have shown that scaling is an appropriate procedure to extend the accuracy of non-linear measures. Additionally, as we will show in the next section, the EEG study we present here reinforces the use of scaling in non-linear LZC computations. 4.3. Classical LZCN and multiscale LZCw in the EEGs of EC and EO Complexity analyses on the EEGs from healthy participants with EC and EO showed that LZCN failed to show any difference between EC and EO conditions. However, LZCw estimations from w 6113 were lower for EC than EO conditions. LZCw results agree with linear based analyses of resting EC vs EO EEGs. Barry et al. (2007) found increased power values at all frequency bands when EC. In addition, non-linear estimations using non-linear detrended fluctuation analysis (Gao et al., 2008) suggested complexity changes at multiple scales. Our results with LZCw show the same pattern, indicating that complexity was reduced at many scales (multiple rhythms). It is worthy to remark that LZC from EC showed a rapid decay reaching a minimum about LZC21 and LZC19 which corresponds to alpha-like rhythms (see Fig. 4a). As mentioned above, this significant decrease in complexity is closely related with the classical alpha power increase during EC, and it strongly indicated that LZCw can be accurate at capturing signal properties from different scales. Interestingly, LZCN did not change with EC and EO, probably because slow oscillations did not change in complexity, and therefore, they masked the variations at other scales. The main discussion here, however, it is not why there are not changes in complexity at slow rhythms. We believe that the fact that LZCN was not powerful enough to discriminate between EC and EO in this particular dataset, while LZCw was able to obtain a complexity spectrum that parallels the main characteristics found with classical frequency analyses. This is consistent evidence that support the multiscale LZCw. Acknowledgements The manuscript was improved thanks to the comments by Fer Mesman and two anonymous reviewers. Conflict of interest: None of the authors have potential conflicts of interest to be disclosed. References Aboy M, Hornero R, Abásolo D, Álvarez D. Interpretation of the Lempel–Ziv complexity measure in the context of biomedical signal analysis. IEEE Trans Biomed Eng 2006;53:2282–8. Andrzejak RG, Lehnertz K, Mormann F, Rieke C, David P, Elger CE. Indications of nonlinear deterministic and finite-dimensional structures in time series of brain

electrical activity: dependence on recording region and brain state. Phys Rev E 2001;64:061907. Barry RJ, Clarke AR, Johnstone SJ, Magee CA, Rushby JA. EEG differences between eyes-closed and eyes-open resting conditions. Clin Neurophysiol 2007;118:2765–73. Barry RJ, Clarke AR, Johnstone SJ, Brown CR. EEG differences in children between eyes-closed and eyes-open resting conditions. Clin Neurophysiol 2009;120:1806–11. Buzsáki G, Draguhn A. Neuronal oscillations in cortical networks. Science 2004;304:1926. Catarino A, Churches O, Baron-Cohen S, Andrade A, Ring H. Atypical EEG complexity in autism spectrum conditions: a multiscale entropy analysis. Clin Neurophysiol 2011;122:2375–83. Costa M, Goldberger AL, Peng CK. Multiscale entropy analysis of complex physiologic time series. Phys Rev Lett 2002;89:068102. Costa M, Goldberger AL, Peng CK. Multiscale entropy analysis of biological signals. Phys Rev E: Stat Nonlin Soft Matter Phys 2005;71:021906. Dauwels J, Srinivasan K, Reddy MR, Musha T, Vialatte FB, Latchoumane C, Jeong J, Cichocki A. Slowing and loss of complexity in Alzheimer’s EEG: two sides of the same coin? Int J Alzheimers Dis 2011;2011:539621. Delorme A, Makeig S. EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics. J Neurosci Methods 2004;134:9–21. Fernández A, López-Ibor MI, Turrero A, Santos JM, Morón MD, Hornero R, Gómez C, Méndez MA, Ortiz T, López-Ibor JJ. Lempel–Ziv complexity in schizophrenia: a MEG study. Clin Neurophysiol 2011;122:2227–35. Fernández A, Zuloaga P, Abásolo D, Gómez C, Serra A, Méndez MA, Hornero R. Brain oscillatory complexity across the life span. Clin Neurophysiol 2012;123:2154–62. Gao T, Wu D, Yao D. EEG scaling difference between eyes-closed and eyes-open conditions by detrended fluctuation analysis. In: Wang R, Shen E, Gu F, editors. Advances in cognitive neurodynamics. The Netherlands: Springer; 2008. p. 501–4. Gómez C, Hornero R, Abásolo D, Fernández A, López M. Complexity analysis of the magnetoencephalogram background activity in Alzheimer’s disease patients. Med Eng Phys 2006;9:851–9. Holm S. A simple sequentially rejective multiple test procedure. Scand J Stat 1979;6:65–70. Hornero R, Abásolo D, Escudero J, Gómez C. Non-linear analysis of EEG and MEG in patients with Alzheimer’s disease. Philos Trans R Soc A 2009;367:317–36. Jouny CC, Bergey CK. Characterization of early partial seizure onset: frequency, complexity and entropy. Clin Neurophysiol 2012;123:658–69. Kaspar F, Schuster HG. Easily calculable measure for the complexity of spatiotemporal patterns. Phys Rev A 1987;36:842–8. Kolmogorov AN. Three approaches to the quantitative definition of information. Probl Inf Transm 1965;1:1–7. Lempel A, Ziv J. On the complexity of finite sequences. IEEE Trans Inf Theory 1976;IT-22:75–81. Li Y, Tong S, Liu D, Gai Y, Wang X, Wang J, et al. Abnormal EEC complexity in patients with schizophrenia and depression. Clin Neurophysiol 2008;119:1232–41. Liang WK, Lo MT, Yang AC, Peng CK, Cheng SK, Tseng P, Juan CH. Revealing the brain’s adaptability and the transcranial direct current stimulation facilitating effect in inhibitory control by multiscale entropy. Neuroimage 2014;90:218–34. Nagarajan R. Quantifying physiological data with Lempel–Ziv complexity – certain issues. IEEE Trans Biomed Eng 2002;49:1371–3. Niedermeyer E. The normal EEG in the waking adult. In: Niedermeyer E, Lopes da Silva E, editors. Electroencephalography: basic principles, clinical applications and related fields. Baltimore MD: Lippincott Williams & Wilkins; 1999. p. 149–73. Pincus S. Approximate entropy as a measure of irregularity for psychiatric serial metrics. Bipolar Disord 2006;8:430–40. Sanei S, Chambers J. Introduction to EEG. In: Sanei S, Chambers J, editors. EEG signal processing. New York: Wiley; 2007. p. 1–35. Solomonoff RJ. A formal theory of inductive inference. Inf Control 1964:224–54. Stam CJ. Nonlinear dynamical analysis of EEG and MEG: review of an emerging field. Clin Neurophysiol 2005;116:2266–301. Takahashi T, Cho RY, Murata T, Mizuno T, Kikuchi M, Mizukami K, Kosaka H, Takahashi K, Wada Y. Age-related variation in EEG complexity to photic stimulation: a multiscale entropy analysis. Clin Neurophysiol 2009;120:476–83. Takahashi T, Cho RY, Mizuno T, Kikuchi M, Murata T, Takahashi K, Wada Y. Antipsychotics reverse abnormal EEG complexity in drug-naïve schizophrenia: a multiscale entropy analysis. Neuroimage 2010;51:173–82. Watanabe TAA, Cellucci CJ, Kohegyi E, Bashore TR, Josiasssen RC, Greenbaun NN, Rapp PE. The algorithmic complexity of multichannel EEGs is sensitive to changes in behavior. Psychophysiology 2003;40:77–97. Zhang XS, Roy RJ, Jensen EW. EEG complexity as a measure of depth of anesthesia for patients. IEEE Trans Biomed Eng 2001;48:1424–33.