The calculation of the first positive Lyapunov exponent in sleep EEG data

Electroencephalography and clinical Neurophysiology, 86 (1993) 348- 352 © 1993 Elsevier Scientific Publishers Ireland, Ltd. 0013-4649/93/$06100

348

E E G 92647

Short communication

The calculation of the first positive Lyapunov exponent in sleep EEG data J. R6schke, J. Fell and P. B e c k m a n n a Department of Psychiatry and a Department of Physics, University of Mainz, Mainz (Germany) ( A c c e p t e d for p u b l i c a t i o n : 2 F e b r u a r y 19931

Summary To help determine if the E E G is quasiperiodic or chaotic we performed a new analysis by calculating the first positive Lyapunov exponent L l from sleep E E G data. Lyapunov exponents measure the m e a n exponential expansion or contraction of a flow in phase space. L~ is zero for periodic as well as quasiperiodic processes, but positive in case of chaotic processes expressing t h e sensitive dependence on initial conditions. We calculated L t for sleep E E G segments of 15 healthy male subjects corresponding to sleep stages !, 1I, III, IV and R E M (according to Rechtschaffen and Kales). Our investigations support the assumption that E E G signals are neither quasiperiodic waves nor simple noise. Moreover, we found statistically significant differences between the values of L i for different sleep stages.

Key words: Sleep; EEG; Lyapunov exponent

The theory of nonlinear dynamic systems provides .some new methods to handle complex systems. In recent years the application of concepts from this theory to sleep E E G data has produced encouraging results. U n d e r selected conditions nonlinear dynamic systems depending on at least 3 state variables are able to generate so-called deterministic chaos. In this case the dynamics show a sensitive dependence on initial conditions, which means that different states of a system being arbitrarily close initially will become macroscopically separated after sufficiently long times. This property can be quantified by calculating the system's Lyapunov exponents, which estimate the mean exponential divergence or convergence of nearby trajectories in phase space. Regardless of the description of a system's dynamics in terms of differential equations, the behavior of such chaotic systems is not predictable over long time periods. The calculation of the correlation dimension D 2, which is a measure of the complexity and which also estimates the degrees of freedom of a signal under study, has shown that in normal healthy subjects the deeper the sleep the lower the E E G ' s dimensionality (Ehlers et al. 1991; R6schke and Aldenhoff 1991, 1992a,b). Moreover, in psychosis and under the influence of drugs alterations of the E E G ' s dimensionality have also been reported (Koukkou et al. 1992a,b; R6scbke 1992a,b). Nevertheless, the interpretation of these results is not obvious. There remain many open questions in the evolving field of chaotic brain dynamics. O u r investigations of the first positive Lyapunov exponents for sleep E E G provide further evidence that the h u m a n electroencephalogram is a chaotic signal.

Correspondence to." Joachim R6schke, M.D., Ph.D., Department of Psychiatry, University of Mainz, Untere Zahlbacher Str. 8, 6500 Mainz (Germany).

Methods To each E E G segment we applied the reconstruction procedure by embedding the signal into a 10-dimensional phase space following the proposal of Takens f1981): y(t) = ( x j ( t ) , x : ( t + r ) . . . . . x i ( t + 2 n r ) ) For the time mcremenl "r we used the first zero crossing of the autocorrelation function. This choice of r insures that two successive delay coordinates are as independent as possible. Excessively short delays would compress the attractor to the diagonal of the reconstruction space, whereas for excessively large values of ~- the structure of the attractor would disappear. We calculated the first Lyapunov exponent L~ applying a modified version of the Wolf algorithm (Wolf et al. 19851 following a proposal of Frank et al. 11990). Essentially, the Wolf algorithm iteratively computes the vector distance L of two nearby points and evolves its length for a certain propagation time. After m propagation steps the first positive Lyapunov exponent results from the sum over the logarithmical ratios of the vector distances. kt=(1/(tm-t,~))

~

Iog2(L'(tk)/L(tk

111

k-I

The search for a replacement point extends on a cone centered about the previously evolved vector. Points lying outside the region over which dynamics are assumed to be linear and points lying closer than the average noise level are discarded. If no point has been found, the angle of the cone Is expanded and the search goes on. Frank et al. (1990) propose another displacement technique introducing a priority function, which depends on the replacement length r and the orientation change 0: p(r,0) = ( ~x +/3 × ( ( b -

r)/(b-a))'

) xcos 0

with ct = 0.1. /3 = 0.9 and y = 3.0 (b: scalmax, a: noise level).

C A L C U L A T I O N O F L Y A P U N O V E X P O N E N T IN SLEEP E E G In our experience (with sleep E E G data sets of 16,384 points) the modified algorithm is less time consuming, converges faster to the final value and yields systematically higher results than the original Wolf algorithm. A crucial problem concerning algorithms similar to the one proposed by Wolf and coworkers is to find reasonable input parameters. Fig. 1 shows how the calculation of L 1 for a single E E G segment (sleep stage IV) depends on the embedding dimension, the noise scale, the m a x i m u m scale and the evolving time. In practice, the minimal embedding dimension needed for the estimation of L 1 is unknown (Grassberger et al. 1991). According to Takens it should be in the range of 2 d + 1 , where d is the dimension of the dynamic system. Realistic values for the average noise levels were extracted from graphs of the correlation function, which were used to calculate the correlation dimension. Intermediate knees for large embedding dimensions are related to noise contamination (Principe and Lo 1991). T h e m a x i m u m scale is the upper bound for the distance of points considered for replacement. It depends on the m a x i m u m amplitude of the signal. We computed L 1 for several values of scaimax (for a single E E G segment) to estimate a region where the calculations are approximately stationary. W e expressed scalmax in terms of the m a x i m u m possible distance (maxdist) in n-dimensional phase space corresponding to the maximal amplitude of the E E G epoch under study. In order to estimate the evolving time the low frequency part of the power spectra can serve for an approximation• We chose an evolving step of 40 data points, corresponding to a frequency of 2.5 Hz. Altogether, we have selected the following input

349 parameters for our calculations: scalmin (noise level)= 10/~V, scalmax = 15% maxdist, evolv = 40 samples, embedding dimension = 10. For each calculation we performed 200 replacement steps.

Materials

W e investigated 15 healthy male subjects aged between 23 and 63 years. Subjects were volunteer recruits from the university student population and the general public. All were in self-reported good health with regular sleep-wake patterns. There was no evidence of hypnotic drug use or above average alcohol or caffeine c o n s u m p t i o n All were free of a past history or current symptoms of psychopathology as well as of any medical condition known to influence sleep. T h e registration of the sleep E E G was started at 11:00 p.m. and was finished at 7:00 a.m. next day. Surface electrodes were placed on the scalp at Pz, Cz, C3 and C4 (10-20 system) and mastoid to record E E G activity, at the outer canthi of the left and right eyes to record eye movements and on the chin to record submental electromyographic activity. Interelectrode impedances were all below 5 kO. For visual analysis according to Rechtschaffen and Kales (1968) the sleep E E G was recorded using a Schwarzer E E G machine (ES 12000). Additionally, the E E G data were digitized by a 12-bit analog-digital converter, sampled with a frequency of f = 100 Hz (50 Hz low pass filter, 48 d B / o c t a v e ) and stored on the disk of a Hewlett Packard computer (A 900). According to Rechtschaffen and Kales

dependence on embedding dim. 6

dependence on scalmin

PrkxdlNdLyupunov-mtponontL1

Principal Lyaptmmt..expcment

$

LI

2.6

2 ~~ 32-'

o 4

"i" I

1,5

lilii'r ,if f a

lo 12 embedding

14 is dimension

is

""

. . . . . . . . .

!

.... " ....

I.

1 0,6

•-I-.

.............

i. . . . . . . . . . . .

0 2o

0

2

4

dependence on scalmax

6

8

10

12 14 16 18 s c a l m l n (uV)

20

22

24

26

28

30

dependence on evolving time

I~rklcipaJ Lympunov-exponen! L1

Princlaal L y s p u f l o v - e x p o n e n t L1

6

lo

I! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

t 4 ...... l ............................................................

/.

o

o 0

10

15 ecalmax

20

26

(% maxdtet)

30

35

40

o

0.2

0.4

o,s

o.s

1

1.2

t.4

~e

evolving time (e)

Fig. 1. D e p e n d e n c e to L 1 on the embedding dimension, the noise scale (scalmin), the m a x i m u m scale (scalmax) and the evolving time exemplarily for sleep stage IV.

350

J. ROSCHKE ET AL.

the sleep E E G was scored by two independent judges. They determined 5 artifact-free time periods each of n = 16,384 data points unambiguously corresponding to one of the sleep stages I, II, III, IV and REM. The representative collections of the sleep stages were from the first half of the night, except the R E M periods.

Results

Before systematically studying sleep E E G segments we calculated L t for synthetic signals composed of 16,384 data points (same data length as the EEG segments and same input parameters). In Fig. 2

periodic

signal

L1 2,5 2

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1,5

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

095

-0,5

~

.

.

.

.

.

' 0

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

' 10

.

.

20

.

.

.

.

.

.

.

.

.

30

.

.

.

.

.

.

.

.

.

40

.

.

.

.

.

.

.

.

.

.

50

60

70

80

t (S)

running

Lyap.-exp.

quasiperiodic

signal

L1 0,5 0

-

-0,5

-1 -1,5 -2 -2,5 -3 -3,5

I

i

10

0

i

20

I

30

i

40

,I

50

,,

I

70

60

80

t (s) running

Lyap.-exp.

white noise L1 6 --

5

-

4

.

3

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

1

.

.

.

.

.

.

.

0

.

.

2

0

f

:----- "~-~ . . . . . . . . . . . . .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

" .... .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

40

.

.

.

30

.

.

.

20

.

.

.

10

.

.

.

I

.

.

.

I

.

.

.

I

I

: ....

.

.

I

ii

: .............

I

so

I

60

i

70

so

t (s)

running

Lyap.-exp.

Fig. 2. Running exponents for a periodic signal (sum of two sine waves with identical amplitudes, fi = 4 Hz and f2 = 8 Hz), a quasiperiodic signal (fl = 4 Hz and f2 = (30)1/2 H z ) and white noise.

CALCULATION OF LYAPUNOV EXPONENT IN SLEEP EEG First positive Lyapunov-exponent L1 4

3,5

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

I

I

II

II

III

IV

III

IV

RI=M Rt=M

sleep s t a g e m

lead position Cz

~

lead position Pz

Fig. 3. Mean values and standard deviations of L t during sleep stages I, II, III, IV and REM for 15 healthy male subjects.

the running exponents for a periodic signal, a quasiperiodic signal (frequencies in the range of dominant EEG frequencies) and white noise are demonstrated. The periodic as well as the quasiperiodic signal are composed of two sine waves with identical amplitudes. For the periodic and the quasiperiodic signal, the values shrink rapidly to L l < 0.01 (theoretically L 1 = 0). For white noise the algorithm yields L 1 = 5.26. For an infinite sequence of data points the theoretical value for L 1 is infinity. The results for all subjects and for lead positions C z and Pz a r e shown in Fig. 3. There are no statistically significant differences between the average L l for the lead positions C Z and Pz (Student's t test). The individual values of L 1 vary between L 1 = 1.23 and L 1 = 3.61. The average L 1 are highest during sleep stage I and diminish during stages II and III. For stage IV the average L 1 lie close to the value of stage III. The average value for paradoxical sleep (REM) is slightly lower than that for stage I. The relative standard deviations (interindividual variations) vary between 8% and 19% of the average values. They are lowest for sleep stages I and REM. For lead position C z the interindividual variations are lower than for Pz, except for sleep stage REM. The average L 1 during sleep stage I and REM for both lead positions are statistically significantly higher (t test, P < 0.05) than during stages II, III and IV. Additionally, for lead position C z the average L 1 during sleep stage II is statistically significantly higher (P < 0.05) than during stages III and IV. We may say that the results for lead position C z exhibit the tendency that the deeper the sleep the lower L 1 with the restriction that the values during sleep stage III are not statistically significantly different from those during stage IV. At position Pz we cannot find any statistically significant differences between the stages II, III and IV.

Discussion

Our results support the hypothesis that EEG signals are neither pure quasiperiodic waves nor simple noise. The values of L 1 lie far away from the results for synthetic quasiperiodic signals and white noise. Nevertheless, quasiperiodicity plus noise and colored noise still remain competing theories. Skarda and Freeman (1987) postulated that behavior can best be modeled as a sequence of ordered, stable states in an evolutionary trajectory. They suppose that chaotic mechanisms enable the neural network to learn new behavior, since without such mechanisms the system could only converge to behavior it has already learned. Our results also demonstrate that there are statistically significant differences between the first positive Lyapunov exponents L 1 of different sleep stages (classified according to Rechtschaffen and

351 Kales). Earlier results concerning the evaluation of the correlation dimension D 2 exhibited similar tendencies. Essentially, the deeper the sleep the lower D 2. That means the degrees of freedom of the s l e e p EEG - in other words the complexity - are reduced continuously the more the sleep moves to slow wave sleep. We would like to mention that the similar behavior of D 2 and L I depending on the deepening of sleep is an absolutely non-trivial outcome. But in contrast to the results for D 2, we did not detect statistically significant differences of L I between all sleep stages. It is an open and challenging question if and how parameters from dynamic system theory like D 2 and L 1 can be interpreted in terms of information processing. Since correlation dimensions mathematically reflect the complexity of signals whereas first positive Lyapunov exponents describe the divergence of trajectories starting at nearby initial states, the question in what way D 2 and L 1 do correspond to complexity and flexibility of information processing is an important one. Compared to conventional spectral analysis, the information gained from a signal is reduced in a drastic way by the application of the analytical methods of nonlinear system theory. This kind of analysis leads to a single value, which reflects a property of the whole dynamic system. We have already pointed out that the numerical results for L~ depend on the choice of the input parameters, the data length and, of course, on the implemented algorithm. Consequently, the absolute values of L 1 should be taken with care. A reasonable interpretation should be based on statistically significant differences of L I. This work was supported by the Deutsche Forschungsgemeinschaft (Ro 809/4-1). We would like to thank an unknown reviewer for his meaningful suggestions.

References

Ehlers, C.L., Havstad, J.W., Garfinkel, A. and Kupfer, D.J. Nonlinear analysis of EEG sleep states. Neuropsychopharmacology, 1991, 5: 167-176. Frank, G.W., Lookman, T., Nerenberg, M.A.H., Essex, C., Lemieux, J. and Blume, W. Chaotic time series analyses of epileptic seizures. Physica D, 1990, 46: 427-438. Grassberger, P., Schreiber, T. and Schaffrath, C. Nonlinear time sequence analysis. Int. J. Bifurc. Chaos, 1991, 1: 521-547. Koukkou, M., Lehmann, D., Wackermann, J., Dvorak, I. and Henggeler, B. The dimensional complexity of the EEG in untreated acute schizophrenics, in persons in remission after a first schizophrenic episode and in controls. Schizophr. Res., 1992a, 6: 129. Koukkou, M., Lehmann, D., Wackermann, J. Dvorak, I. and Henggeler, B. Dimensional complexity of EEG brain mechanisms in schizophrenia. Biol. Psychiat., 1992b, 242: 191-196. Principe, J.C. and Lo, P.C. Towards the determination of the largest Lyapunov exponent of EEG segments. In: W.S. Pritchard and D.W. Duke (Eds.), Measuring Chaos in the Human Brain. World Scientific, Singapore, 1991. Rechtschaffen, A. and Kales, A. A Manual of Standardized Terminology, Technics and Scoring System for Sleep Stages of Human Subjects. NIH Publication No. 204. US Government Printing Office, Washington, DC, 1968. R6schke, J. Aspects of the chaotic structure of sleep EEG data. Neuropsychobiology, 1992a, 25: 61-62. R6schke, J. Strange attractors, chaotic behavior and informational aspects of sleep-EEG data. Neuropsychobiology, 1992b, 25: 172176. R6schke, J. and Aldenhoff, J.B. The dimensionality of human's electroencephalogram during sleep. Biol. Cybern., 1991, 64: 307313.

352 R6schke, J. and Aldenhoff, J.B. A nonlinear approach to brain function: deterministic chaos and sleep EEG. Sleep, t992a, 15: 95-101. R6schke, J. and Aldenhoff, J.B. Estimation of the dimensionality of sleep-EEG data in schizophrenics. Eur. Arch. Psychiat. Neurol. Sci., 1992b, 242: 191-196.

J. ROSCHKE ET AL. Skarda, C. and Freeman, W. How brains make chaos in order to make sense of the world. Behav. Brain Sci., 1987, 10: 161-195. Takens, F. Detecting Strange Attractors in Turbulence. Lecture Notes in Mathematics. Springer, Berlin, 1981:898 pp. Wolf, A., Swift, J.B., Swinney, H.L. and Vastano, J.A. Determining Lyapunov exponents from a time series. Physica D, 1985, 16: 285.