ARTICLES Chinese Science Bulletin 2006 Vol. 51 No. 24 3059—3064

DOI: 10.1007/s11434-006-2213-y

A new measure to characterize multifractality of sleep electroencephalogram MA Qianli1,2, NING Xinbao1, WANG Jun2 & BIAN Chunhua1 1. State Key Laboratory of Modern Acoustics, Department of Electronics Science and Engineering, Institute for Biomedical Electronic Engineering, Nanjing University, Nanjing 210093, China; 2. Province Key Laboratory of Image Processing & Image Communication, College of Telecommunications and Information Engineering, Nanjing University of Posts & Telecommunications, Nanjing 210003, China Correspondence should be addressed to Ning Xinbao (email: xbning@ nju.edu.cn) Received January 23, 2006; accepted April 7, 2006

Abstract Traditional methods for nonlinear dynamic analysis, such as correlation dimension, Lyapunov exponent, approximate entropy, detrended fluctuation analysis, using a single parameter, cannot fully describe the extremely sophisticated behavior of electroencephalogram (EEG). The multifractal formalism reveals more “hidden” information of EEG by using singularity spectrum to characterize its nonlinear dynamics. In this paper, the zero-crossing time intervals of sleep EEG were studied using multifractal analysis. A new multifractal measure Δasα was proposed to describe the asymmetry of singularity spectrum, and compared with the singularity strength range Δα that was normally used as a degree indicator of multifractality. One-way analysis of variance and multiple comparison tests showed that the new measure we proposed gave better discrimination of sleep stages, especially in the discrimination between sleep and awake, and between sleep stages 3 and 4. Keywords: sleep, electroencephalogram, multifractality.

The research of nonlinear dynamics in electroencephalogram (EEG) has made much headway in recent years. Nonlinear analysis methods have been successfully applied to the studies of brain functions and ― pathological changes in EEG[1 3]. These studies also proved that EEG exhibited at least partly chaotic characteristics. The detrended fluctuation analysis (DFA)[4], which has been widely used recently, revealed the longwww.scichina.com

www.springerlink.com

range power-law correlation in EEG, indicating time scale invariant and fractal structure[5,6]. However, EEG is also rather noisy, displaying short-term decorrelation like white noise, and consequently, the EEG has been traditionally considered as a linear stationary random process. The paradoxical combination of short-term decorrelation and long-range correlation, stochastic and deterministic suggests that a single nonlinear parameter, such as largest Lyapunov exponent, correlation dimension, fractal dimension, scaling exponent, etc., may not be able to fully characterize the “stochastic chaos” (as named by Freeman[7]) nature of EEG. The long time behavior of chaotic, nonlinear dynamic systems can often be characterized by (mono) fractal or multifractal measures. Monofractals are homogeneous in the sense that they have the same scaling properties, characterized by a single singularity exponent throughout the entire signal. In contrast, multifractals can be decomposed into many (possibly infinite) sub-sets characterized by different exponents. Multifractal signals are intrinsically more complex and inhomogeneous than monofractals. Multifractal models have been used to account for scale invariance properties of various objects in very different domains ranging from the energy dissipation or the velocity field in turbulent flows[8] to underlying hierarchical structure in proteins[9]. Physiologic signals generated by complex self-regulating systems, such as heartbeat interval, electrocardiogram, gait etc. have been proven to be multifractal, and the degree of multifractality often relates to pathological state or natural aging process[10–12]. The multifractal formalism has been applied to EEG analysis recently[13,14]. EEGs taken from healthy subjects and epilepsy patients, during imaginary and real visual-motor tracking tasks, and from different recording sites have been studied, and the results suggest that the multifractal formalism might be a good method for characterizing EEG dynamics. In most of the present studies, nonlinear analysis was applied directly to the EEG amplitude fluctuation time series acquired by discrete-time sampling. However, EEG is a typical signal which is particularly susceptible to biological and experimental artifacts that manifest themselves as large amplitude deviations, and this leads to a questioning of the validity of techniques based on generalizations of geometric box counting, such as various approaches of the multifractal analysis. In this paper, we used the zero-crossing time intervals (CTI), 3059

ARTICLES which were derived from zero-crossings of EEG time series, to study the multifractality of EEG. 1

Zero-crossing time interval series of EEG

Zero-crossing series is a frequency-based series, defined as the set of zero-crossing events. Zero-crossings are robust to amplitude artifacts. And in addition, because it contains only frequency information, it is useful to distinguish the contributions from frequency and/ or amplitude separately and to understand whether the frequency profile alone has fractal characteristics, independent of amplitude changes. Zero-crossings have been used historically to automate EEG analysis over long time-periods, because they can clearly distinguish low frequency and high-frequency activity by measuring the time between each crossing with respect to the frequency band of interest[15]. However, despite the long research history, using zero-crossings to study nonlinearity in EEG dynamics appears to be scarce. Watters and Martins[6] used DFA to analyze the “EEG walk” constructed from the zero-crossings of EEG. They observed a mean scaling exponent across all subjects and sites of α = 0.67, which indicates time scale invariance and fractal structure of EEG zero-crossings. In the present study, to adapt to multifractal analysis, we used another series which was also constructed from zero-crossings. Let the EEG be x(t). The zero-crossing time is the level set {ti, x(ti) = 0} where the index i registers the order of the zero crossing event. In practice, {ti} is first determined by linear interpolation and then used to define the set of crossing time intervals (CTI) C = {τi = ti+1 − ti}. Fig. 1 shows the amplitude fluctuations, zero-crossings and CTI series of a segment of EEG signal. Frequency activities are evidently exhibited in

the graph of zero-crossings (Fig. 1(b)) and CTI series (Fig. 1(c)). The CTI of the fractional Brownian motion has been proven to follow a power law distribution[16]: p(τ) ~ τ-v, where p(τ) is the probability density function (PDF). Fig. 2 shows the PDF estimate of EEG CTI. It can be seen from the figure that the PDF also follows a power law distribution, suggesting fractal structure in EEG CTI. 2

Multifractal singularity spectrum

The multifractal singularity spectrum was first proposed by Halsey et al. [17] to describe normalized distributions (measures) lying upon possibly fractal sets, for example those arising from dynamical systems theory. They focus upon the scaling properties of such measures, by considering their singularities, which are characterized by two indices: α (Lipschitz-Holder exponent), which determines the strength of their singularities; and f, which describes how densely they are distributed. The spectrum of singularities is described by giving the possible range of α values and the function f(α). The measures mentioned above correspond to varied physical quantities according to the specific problem under study, such as CTI (time intervals) and amplitude fluctuations (voltage) of EEG in this paper. Suppose that we cover the support of the measure with boxes of size L, and the number of the boxes is N. Then N~L−1. Define Pi (L) = Mi(L)/MT as the probability in the ith box, where Mi(L) is the integrated measure in the ith box and MT is the total amount of the measures. The singularity strength α is defined as Pi ( L) ~ Lαi . (1) The smaller the exponent α, the more singular the

Fig. 1. Amplitude fluctuations (a), corresponding zero-crossings (b) and CTI series (c) of a segment of EEG signal.

3060

Chinese Science Bulletin

Vol. 51 No. 24 December 2006

ARTICLES amplifies the more singular region of P, while for q<1 it accentuates the less singular regions, and for q=1 the measure μ(1) replicates the original measure. Then the Hausdorff dimension of the measure theoretical support of μ(q), i.e. f , is given by 1 f (q) = − lim ∑ i μi (q, L) log[μi (q, L)] N →∞ log N = lim

(6) log L In addition, we can compute the average value of the singularity strength αi=logPi/logL with respect to μ(q) by evaluating 1 α (q) = − lim ∑ i μi (q, L) log Pi ( L) N →∞ log N L →0

Fig. 2. Log-log plot of estimated (○) p(τ) of EEG CTI series. The solid line shows least square fit to the PDF.

measure, and the “stronger” the singularity. It is a characterization of the scaling in the ith region or spatial location. If we count the number of the boxes N(α) where the probability Pi has singularity strength between α and α+dα, then f (α) can be loosely defined as the fractal dimension (box dimension or Hausdorff dimension, to describe the density of the distribution) of the set of boxes with singularity strength α by

N (α ) ~ L− f (α ) .

(2) Large errors will be induced if we estimate α and f directly from eqs. (1) and (2) or by Legendre transform from the generalized dimensions Dq. Chhabra and Jensen[8] proposed an algorithm for determining the f(α) directly from experimental data. We briefly introduce their algorithm in the following. Please refer to ref. [8] for more details. First, the entropy S of the process from which the measures arise is given by

S = −∑ Pi log Pi ,

(3)

i

and the Hausdorff dimension dh of the measure theoretical support of the measure can be related to the entropy by a theorem by Billingsley which gives 1 N (4) d h = − lim ∑ Pi log Pi . N →∞ log N i =1 We can use these results to evaluate f (α) for a multifractal measure P(x). This is done by first constructing a one-parameter family of normalized measures μ(q) where the probabilities in the boxes of size L are

μ (q, L) = [ Pi ( L)]q

∑ [ Pj ( L)]q .

(5)

j

The parameter q provides a microscope for exploring different regions of the singular measure. For q>1, μ(q) www.scichina.com

www.springerlink.com

∑ i μi (q, L) log[μi (q, L)] .

= lim

∑ i μi (q, L) log Pi ( L) .

(7) log L Eqs. (6) and (7) provide a relationship between a Hausdorff dimension f and an average singularity strength α as implicit functions of the parameter q. The range of the singularity strength is defined as Δα = α max − α min , (8) L →0

where α max = α (q → −∞) and α min = α (q → +∞). Δα is a measure of the deviation from monofractal, indicating the degree of multifractality. Smaller Δα indicates that the measure tends to be monofractal, and contrarily, larger Δα indicates multifractality[12]. 3

Application to sleep EEG analysis

The experimental data used in this study come from PhysioBank “The Sleep-EDF Database”[18,19], which is a collection of sleep recordings from 8 healthy subjects. The recordings were obtained from Caucasian males and females (21―35 years old) without any medication. Each record contains horizontal EOG, FpzCz and PzOz EEG, sampled at 100 Hz, and is accompanied by a hypnogram. Hypnograms are sleep staging results manually scored for every 30 s epoch according to Rechtschaffen & Kales. We applied multifractal analysis to CTI series of sleep EEG. Fig. 3 shows an example of the linear fit to Σμilogμi vs. log L with varied q from the bottom in the order of q = 1, 3 and 5. It can be seen from the figure that there is no ambiguity in determining the slopes, and different q correspond to different slopes. This confirms multifractality in EEG CTI series. Then the multifractal singularity spectra were estimated for every sleep stage. For convenience, we will 3061

ARTICLES defined as

Δ asα ≡ Δα − − Δα + ,

Fig. 3. Σμilogμi vs. logL with q = 1 (◊), 3 (□) and 5 (○). The solid lines show linear square fits correspondingly.

use “Wake” for denoting the state of wakefulness, “S1”―“S4” for sleep stages 1―4 and REM for rapid eye movement sleep in the following description and figures. Fig. 4 shows the singularity spectra of varied sleep stages. It can be observed from the figure that the spectra are asymmetric, and there are differences in the spectra’s asymmetry among sleep stages. In the state of wakefulness, the region of q>0 (corresponding to smaller α) is larger than the region of q<0 (corresponding to larger α), i.e. the singularity spectrum is more concentrated on the “stronger” region of singularity. And with the transition into deeper sleep, the region of q>0 decreases, and contrarily the region of q<0 increases, i.e. the singularity spectrum moves towards the “weaker” region of singularity. The asymmetry of REM is between S1 and S2.

where Δα− = αmax –α0 is the range of the singularity strength with q<0, Δα+ = α0 –αmin is the range of the singularity strength with q>0, and α0 = α (q = 0) is the singularity strength where the singularity spectrum obtains its maximum. We studied 24-hour EEG records of four subjects (sc4002e0, sc4012e0, sc4102e0 and sc4112e0). The signals were divided into segments of 30 s for being consistent with the hypnograms. The total number of validated segments is 11314 (states distribute as, Wake: 7722; S1: 286; S2: 2036; S3: 289; S4: 240; REM: 741). Then CTI series was constructed for each segment, and was analyzed with multifractal formalism. Δasα and Δα were estimated for quantifying changes of multifractal singularity spectra. The hypnogram and the corresponding Δasα and Δα of subject sc4012e0 are plotted in Fig. 5 as an illustration of their variation with sleep state transition. We grouped analysis results of all the four subjects according to sleep stages given in the database. Fig. 6 shows the means and standard deviations of the two measures at various sleep stages. We did one-way analysis of variance (ANOVA) to test the significance of the difference among the stages. To balance the group sizes, for each group of Wake, S2 and REM, 300 segments were randomly selected, and combined with all the available segments of S1, S3 and S4, and then tested with ANOVA. The test of Δasα came out with F(5, 1709) = 394.41, and Δα with F(5, 1709) = 287.96. To further study the discrimination of the states given by each measure, we also did multiple comparison test. The results show that, at a significance level of P<0.05, both Δasα and Δα cannot distinguish between the states of S1 and REM, and in addition, Δα also cannot distinguish between the states of S3 and S4. 4

Fig. 4. stages.

Multifractal singularity spectra (f(α) vs. α) of varied sleep

We introduce a measure Δasα to describe the asymmetry of EEG multifractal singularity spectrum and its variation with sleep state transition which is 3062

(9)

Discussion and conclusion

In this paper, we studied the CTI series of EEG using multifractal formalism. We found evidence of multifractal structures in EEG CTI series, suggesting multifractality of EEG. Multifractal analysis was applied to EEG CTI series of varied sleep stages. We observed that the multifractal singularity spectrum of EEG CTI series was asymmetric, and there are differences in the spectra’s asymmetry among sleep stages. In the state of wakefulness, the singularity spectrum is more concenChinese Science Bulletin

Vol. 51 No. 24 December 2006

ARTICLES

Fig. 5. The hypnogram (a) and the corresponding Δasα (b) and Δα (c) of subject sc4012e0. The arrows denote short-term wakefulness in sleep.

Fig. 6. Means and standard deviations of Δasα (a) and Δα (b) of all the four subjects at various sleep stages.

trated on the “stronger” region of singularity. And with the transition into deeper sleep, the singularity spectrum moves towards the “weaker” region of singularity. So we proposed a new measure Δasα to describe the asymmetry of EEG multifractal singularity spectrum. We compared the new measure with the singularity strength range Δα which was normally used as a degree indicator of multifractality. This was done by applying them to the discrimination of sleep stages. The results of ANOVA and multiple comparison test indicates that, Δasα gives better discrimination of sleep stages overall, and also it gives good discrimination of the states S3 and S4 that Δα fails to do. However both of the two measures cannot distinguish between the states of S1 and REM, and it is difficult to automatically discriminate these two states at all times[20]. www.scichina.com

www.springerlink.com

Carefully examining the curves of two measures in Fig. 5 and the statistic results in Fig. 6, we can find that the value’s fluctuation of Δasα in the state of wakefulness is smaller than that of Δα. And Δasα is obviously more sensitive to short-term wakefulness in sleep (denoted as the arrows in Fig. 5). In fact, both of the two measures can characterize the variation of EEG’s nonlinear properties with sleep state transition, and both have some limitation on that. Multifractal singularity spectrum has much information on nonlinear properties of the system. It is worth to be further studied that how to find effective parameter from the spectrum to characterize brain functions and pathological changes in EEG. Biological processes of human are extremely sophisticated. Sleep is not simply a state in which the brain is 3063

ARTICLES resting, but a dynamic, complicated condition during which the brain is quite active. The sleep stages can hardly be discriminated accurately by any single analysis method or any single measure at the present time. The new method and measure we proposed in this paper is not an exception that can distinguish sleep states accurately by itself. The R&K criteria which have been considered as a “golden standard” for sleep stage discrimination, also needs multiple bio-signals, and multiparameters for manual scoring. Nevertheless, results of our method show that there are significant differences among sleep stages at the significant level of P<0.05. The purpose of our study is to find new measures by utilizing latest advances in nonlinear dynamics. We hope that by combining our new measure with linear/ nonlinear parameters that have been proven to be effective, the automatic discrimination accuracy of sleep stages will be improved and meet the requirements of clinical practices finally. Acknowledgements This work was supported by the National Natural Science Foundation of China (Grant No. 60501003).

References 1

Thakor N V, Tong S. Advances in quantitative electroencephalo-

2

Acharya U R, Faust O, Kannathal N, et al. Non-linear analysis of

6

7 8

9 10 11 12 13

14

15 16

17

gram analysis methods. Annu Rev Biomed Eng, 2004, 6: 453―495 18

EEG signals at various sleep stages. Comput Methods Programs Biomed, 2005, 80(1): 37―45 3

Jiang Z H, Feng H Q, Liu D L, et al. Analyzing sleep EEG using correlation dimension and approximate entropy. J Biomed Eng (in

19

Chinese), 2005, 22(4): 649―653 4

Peng C K, Buldyrev S, Goldberger A L, et al. Long-range correlations in nucleotide sequences. Nature, 1992, 356: 168―170

5

Watters P A. Time-invariant EEG power laws. Intern J Syst Sci, 2000, 31: 819―826

3064

20

Watters P A, Martin F. A method for estimating long-range power law correlations from the electroencephalogram. Biol Psychol, 2004, 66(1): 79―89 Freeman W J. A proposed name for aperiodic brain activity: Stochastic chaos. Neural Netw, 2000, 13(1): 11―13 Chhabra A B, Meneveau C, Jensen R V, et al. Direct determination of the f(alpha) singularity spectrum and its application to fully developed turbulence. Phys Rev A, 1989, 40(9): 5284―5294 Balafas J S, Dewey T G. Multifractal analysis of solvent accessibilities in proteins. Phys Rev E, 1995, 52(1): 880―887 Ivanov P C, Amaral L A, Goldberger A L, et al. Multifractality in human heartbeat dynamics. Nature, 1999, 399(6735): 461―465 West B J, Scafetta N. Nonlinear dynamical model of human gait. Phys Rev E, 2003, 67(5): 051917 Wang J, Ning X, Ma Q. Multiscale multifractality analysis of a 12-lead electrocardiogram. Phys Rev E, 2005, 71(6): 062902 Popivanov D, Jivkova S, Stomonyakov V, et al. Effect of independent component analysis on multifractality of EEG during visual-motor task. Signal Processing, 2005, 85(11): 2112―2123 Wang W, Ning X B, Wang J, et al. Interleaving distribution of multifractal strength of 16-channel EEG signals. Chin Sci Bull, 2003, 48(16): 1700―1703 Smith J. Automated EEG analysis with microcomputers. Med Instrum, 1980, 14: 319―321 Ding M, Yang W. Distribution of the first return time in fractional Brownian motion and its application to the study of on-off intermittency. Phys Rev E, 1995, 52(1): 207―213 Halsey T C, Jensen M H, Kadanoff L P, et al. Fractal measures and their singularities: The characterization of strange sets. Phys Rev A, 1986, 33(2): 1141―1151 Goldberger A L, Amaral L A N, Glass L, et al. PhysioBank, physioToolkit, and physioNet: Components of a new research resource for complex physiologic signals. Circulation, 2000, 101(23): e215―e220 Mourtazaev M S, Kemp B, Zwinderman A H, et al. Age and gender affect different characteristics of slow waves in the sleep EEG. Sleep, 1995, 18(7): 557―564 Grozinger M, Roschke J. The automatic recognition of REM sleep: A challenge and some answers. Methods Find Exp Clin Pharmacol, 2002, 24(Suppl D): 33―35

Chinese Science Bulletin

Vol. 51 No. 24 December 2006