Characterizing Sleep Stages by the Fractal Dimensions of Electroencephalograms
Sy-Sang Liaw1 and Jiann-Yeu Chen2*
1Department of Physics, National Chung Hsing University, Taiwan
2Research Center for Sustainable Energy and Nanotechnology, National Chung Hsing University, Taiwan
Submission: April 03, 2017; Published: July 11, 2017
*Corresponding author: : Jiann-Yeu Chen, Research Center for Sustainable Energy and Nanotechnology, National Chung Hsing University, Taiwan, Email: firstname.lastname@example.org
How to cite this article: Sy-Sang L, Jiann-Yeu C. Characterizing Sleep Stages by the Fractal Dimensions of Electroencephalograms. Biostat Biometrics Open Acc J. 2017;2(2): 555584. DOI: 10.19080/BBOAJ.2017.02.555584
Fourier analysis has shown that human electroencephalography (EEG) signals mainly consist of frequencies below 80 Hz. Waves with a frequency below 13Hz are dominant during sleep. Here we find that the fractal dimensions of sleep EEG results in a frequency range of 2Hz to 13Hz have potential to be used as a continuous parameter to characterize sleep status. Our results show that during sleep, the fractal dimension is high during REM stage and decreases with subsequent sleep stages according to the standard scheme . The lowest value of the fractal dimension for any subject's sleep EEG occurs at N3 (deep sleep).
Electrical potential measured from the scalp of a sleeping subject varies from several to a hundred micro-volts in a non- stationary way. It is known that the time sequence of these signals, called electroencephalograms (EEG), exhibits distinct features during different sleep stages. For example, the spindle and K-complex waves occur during light-sleep stages, while low-frequency (<4Hz) waves account for more than 20% of the frequency distribution during deep-sleep stages. The R & K standard  has long been used to help scorers plot subject hypograms. The hypogram of a typical healthy subject shows that a full night’s sleep comprises quasi-periodical stages including wakefulness, light sleep (N1 and N2), deep sleep (N3), and rapid eye movement (REM, stage 5). A trained scorer typically takes more than an hour to plot the hypogram of an 8-hour sleep session. With hospitals, health centers and research labs producing increasing amounts of EEG data, a fast and automated method for characterizing sleep stages is needed to increase efficiency and reliability.
Several methods have been developed to automatically characterize sleep stages. One method  uses pattern recognition techniques based on the distinct EEG patterns in different stages described in the R&K standards, while another  uses average EEG signal frequencies. Here we find that fractal dimensions might be useful for characterizing sleep EEG signals. Fractal dimensions have been studied in the past to analyze EEG including sleep patterns [4-6]. However, as we will see in our calculations below, no part of the sleep EEG time sequence is a simple monofractal, and their fractal properties vary at different scales. Thus the use of a single fractal dimension as the only parameter for characterizing sleep EEG requires careful interpretation. In particular, the scale-dependent nature of the sleep EEG should be taken into consideration. Previous fractal analyses on the sleep EEG did not discuss the scale-dependent nature of the signals, and should be re-examined to validate their applicability.
We use the modified Inverse Random Midpoint Displacement (mIRMD) method to calculate the fractal dimensions of the sleep EEG. The mIRMD considers the fractal dimension as a quantity to characterize the degree to which a curve turns [7,8]. We have shown  that the fractal dimension of a time sequence calculated by the mIRMD is equivalent to that obtained by the Higuchi method (which calculates the Hurst exponent H1) or DFA (which calculates the Hurst exponent H2) for many monofractals including white noise, Weierstrass functions, and random walk trajectories. A particular merit of the mIRMD is that a time sequence needs not to be detrended for any time interval before actually being used to obtain the Hurst exponent. The trend effect of the first order is automatically eliminated in mIRMD . The algebraic formulation of mIRMD (Eq. (1)) is simple and easily coded. Let us denote a time sequence by f (i), i =0,1,2,3N. f (i) is the value of the signal (in units of H V for sleep EEG) at time i△t, where △t is the constant time interval given by the inverse of the sampling frequency (which is typically 100~250Hz for sleep EEG). The mIRMD calculates the average midpoint displacement of any two values f(i - s) and f(i + s) separated by the time scale 2s△t:
It has been shown  that <G(s)> varies like S2-D if f(i) is a monofractal of dimension D. Thus we can easily find fractal dimension D from the slope of a log-log plot of versus the scale s by the relation D = 2-slope. For example, the mIRMD results for the time sequence of positions of a one-dimensional randomwalker has slope 0.5 (Figure 1), showing that D is 1.5 for a random walk.
average midpoint displacement <G(s)> varies as the scale (s) to the power 2-D. Thus the fractal dimension (D) of a time sequence can be obtained from the slope of the plot: D = 2- slope.
Real time sequences are not monofractals in general. That is, the data points [log(s), log(< G(s) >)] for a real sequence are generally not well fitted by a line. For many cases, they are apparently better fitted by piecewise lines with a crossover (Figure 2), such in heart rate signals [9,10], in fluctuations of fatigue crack growth , in wind speed data , in precipitation and river runoff records , and in stock indexes at one minute intervals [8,13]. We call these sequences crossoverfractals  which have one fractal dimension for small scales below some crossover point and the other fractal dimension for large scales. For other time series, such as those found in DNA sequences  and ECG data , the slope of the log-log plot varies with scale. In general, the concept of fractal dimension is not appropriate for these cases.
On the other hand, previous studies indicate that a real series can still be self-affine within a certain range of scales. Quantitatively, in terms of the mIRMD, this means that the average midpoint displacement of a real sequence <G(s)> varies as S2-D approximately with a constant D for s within a range. We will simply call this D the fractal dimension of the sequence within this range, and regard the series as a representation at the specific range of a fractal of dimension D. Here we use mIRMD to calculate of sleep EEG data and find that the EEG at different sleep stages all behave approximately as a fractal at the same range of time scales. Thus, for practical purposes, we are able to characterize the sleep stages by the value of fractal dimensions within that scale range.
Our sleep EEG data were downloaded from Physionet. com . There are three groups of data. The sc-group is for 4 healthy subjects: sc4002e0, sc4012e0, sc4102e0, and sc4112e0. The sampling frequency is 100 Hz. The second group (st-group) consists of 4 subjects who have slight difficulty in falling sleep: st7022j0, st7052j0, st7121j0, and st7132j0. The sampling frequency is also 100Hz. The third group (slp-group) contains 12 subjects, slp03 and 11 others listed in Table 1, who experience apnea during sleep. The sampling frequency for the third set is 250Hz. Since the signals for apnea occupy only a small portion of whole series (there are few occurrences of apnea during one night's sleep and each occurrence lasts for only 10-15 seconds), they do not significantly affect the average value of G(s). That is, the fractal dimensions obtained by the mIRMD method are still useful in differentiating sleep stages of the subjects regardless of apnea symptoms. We have also started to cooperate with the Psychology Lab of CCU  to collect more sleep EEG for future work. The EEG data from CCU are sampled at 200 Hz. In this study, we have shown one result from CCU data to check the consistence with what we have found from Physionet data.
The sleep EEG data are normally displayed in 30 second increments per page called an epoch. We use mIRMD to calculate and plot its log value versus log(s) to determine the fractal property of the EEG data. Figure 3 shows a typical plot for an epoch of data for the healthy subject sc4112e0 in deep sleep. The general feature of the plot is that it has slopes with values between 0 and 1 for small scales and turns into a horizontal line with small oscillations at large scales. In terms of the concept of the scale-dependent fractals , this plot shows that the sleep EEG behaves like fractals with dimensions between 1 and 2 at small time scales but like white noise at large scales.
We next make three arrangements before carrying out the mIRMD calculations for all EEG data we had. (1) Frist, the fractal property of sleep EEG is scale-dependent, as we have seen in Figure 3. The sleep EEG apparently has good self-affine property in the range 0.08sec to 0.5sec in small time scales, or equivalently in terms of frequency, 13Hz to 2Hz. To support this observation, we use the least-square-method to linearly fit the mIRMD results of all epochs for the subject sc4112e0 EEG data. Good quality of the fitting is supported by the results shown in Figure 4a that the distribution of the coefficient of determination R2 is concentrated near 1. As a comparison, a linear fitting for a larger range of scales, from 0.5Hz to 50Hz, does not produce as good result its distribution of the coefficient of determination R2, shown in Figure 4b, is diverse with an average value about 0.6. Therefore, in what follows we will use the slope of the fitted line in time scales between 0.08 and 0.5 sec to determine the fractal dimension D of the sleep data. Note the frequency components corresponding to this time scale range are between 2 and 13Hz which match the key frequency spectra used by scorers to determine a subject's sleep stages based on the R&K standard. (2) Second, the EEG data downloaded from Physionet actually include values from two different channels, Fpz-Cz and Pz-Oz channels. We find that the fractal properties of the time series of both channels are similar. In Figure 5 we plot the fractal dimension D of each epoch sequentially for both channels against the hypnogram determined by scorers. We see though the values of D for these two channels do not agree with each other, their trends are similar. In particular, many of their local minima occur at the same epochs. In below we find one might characterize sleep stages from their relative rather than absolute D values, thus the EEG data from either one channel is equally good for this purpose. (3) Third, since the mIRMD calculates the average value of the midpoint displacements of a series, the longer a series is used for calculation, the better the results. Thus, to find possibly distinct fractal properties of the sleep EEG at each stage, we connect all epochs signals of the same stage (manually scored according to the R&K standard) into one sequence and calculate its Figure 6. For all sleep stages from 0 to 5, the values of log(Figure 6 we also show the mIRMD results for a subject sleep EEG measured by the Psychology Lab of CCU in Taiwan. Qualitatively, their fractal properties are in consistence with that of the data collected from Physionet as we will describe below.
Let D be the fractal dimension the sleep EEG in time scales between 0.08 and 0.5sec. We calculate the values of D for all epochs and find the average and standard deviation for each stage of data. Results for all subjects are given in Table 1. The subject age and the duration of each sleep stage (in percentage of the total sleep time) are also indicated in the table. From Table 1, we see that, for each subject, we have
D(Wake) ≈ D(REM) ≈ D(N1) >≈ D(N2) > D(N3). (2)
We calculate the average D value of all subjects for every stage. The results are shown in Figure 7 which demonstrates clearly the decreasing order shown in Eq.(2). Although it is still difficult to differentiate the REM stage, light sleep stage 1, and wakeful stage based on the fractal dimension of the EEG only, we see that, for each person, the fractal dimension D diminishes as the subject enters deeper sleep. Of course, the fractal dimensions of the sleep EEG differ from subject to subject, depending on various factors including age, sex, health, etc. A value of D = 1.6, may indicate one subject is in light sleep, but another in deep sleep. But once the fractal profile for a subject is established by a monitoring few full-night's sleep, the value of D can easily indicate whether the subject is in deep sleep.
Our analysis also shows that the subjects in the slp-group generally have larger fractal dimensions than those in the sc- group at any sleep stage (Table 1). It remains unclear whether this is due to the fact that the subjects in the slp-group (age 35- 44) are older than those in the sc-group (age 26-33) or due to their apnea symptoms. Data from more healthy subjects of all ages are needed before one can draw conclusions about the relationship between subject age and their fractal dimensions.
A monofractal time sequence can be constructed with afrequency distribution that satisfies a power law [16,17], or a crossover-fractal with two power exponents .The relation between the fractal dimension D and frequency distribution power β is commonly held to be
β = —(2H 1) = 2D — 5, 1 < D < 2, 0 < H < 1 (3)
where H is the Hurst exponent. Numerical testing of the relation (3) finds that the equality fails significantly for H < 0.1. Nevertheless, the monotonic relation between p and D is good for the whole range 0 < H < 1. Due to the monotonic relation between p and D, one would think the information contained in the fractal dimension can also be revealed in the power spectrum . However, this is not true in practice. For example, a fractal analysis on the S&P500 stock index [8,13] shows that the index time sequence is a crossover fractal: at a small time scale of less than 20 minutes, the stock index has a much smaller fractal dimension than that of its long time variation (Figure 2a). This crossover behavior at the small time scale can hardly be discovered using the power spectrum method by looking at the high frequency components (Figure 2b).
Figure 8 shows the power spectrum of the sleep EEG at different stages for subject sc4012e0 and compares the results with our fractal-dimension analysis. We see that the power spectra exhibit large fluctuations at all frequencies. However, the average power strength decreases roughly linearly between the frequencies 2 to 13 Hz. When a linear-square fit is made for the power spectra within the range 2-13Hz, we obtain the values of the slope p equal to -0.41, -1.24, -1.93, -2.83, -1.88 respectively for sleep stages N0 (awake), N1, N2, N3, and REM. The results are roughly consistent with the values of D (or H) we obtained according to relation (3) except for the sleep stage 0 where the Hurst exponent is close to 0: H = 0.06 < 0.1. The comparisons between the power spectra and the fractal-dimension analysis for other subjects have similar results, showing that the general rule for fractal dimensions given by Eq. (2) is supported by the power spectrum analysis. Namely, the power p between the frequencies 2 to 13Hz in the power spectrum has a relation similar to Eq.(2):
(Wake)>≈ (REM)>≈ (light sleep) > (deep sleep) (4)
However, due to large fluctuations that occur in the power spectra, β is a good parameter for characterizing sleep EEG, but not as good as the fractal dimension D.
In conclusion, we have used the mIRMD method to calculate the fractal dimensions of sleep EEG data sourced from physionet.com. We found that at any sleep stage the EEG is not a monofractal, which was not surprising because real time sequences are generally multifractal and are not expected to satisfy the strict requirements of a monofractal. Nevertheless, we found an interesting property of the sleep EEG, namely that it behaves approximately like a monofractal within the time scale range corresponding to frequency interval 2-13Hz. A value D similar to fractal dimension can be defined reasonably well for sleep EEG patterns within this scale range. We then calculated these D values for the sleep EEG of all 20 subjects and found that the value of D is always the lowest at deep sleep for all subjects. One can thus easily tell whether or not a subject is in deep sleep by monitoring his or her sleep EEG and interactively calculating its fractal dimension D in certain appropriate time scale range. Our calculations also indicate that the fractal dimension D alone cannot distinguish light sleep from REM. To automatically generate a full hypnogram, we need to find other characteristics hidden in the sleep EEG that can be easily detected interactively, and this is the focus of future work. In addition, our calculations show that the fractal dimension D differs between subject groups, most likely due to variation in attributes such as age, sex, or sleep disorders. Future work will collect additional sleep EEG data from selected subjects to study how the fractal dimension D is affected by these attributes.
This work was supported by the National Science Council under the grant number NSC101-2112-M-005-001-MY3, and by the National Center for Theoretical Sciences of Taiwan. We are indebted to the Psychology Lab of National Chung-Cheng University, ChiaYi, Taiwan, which offers us some of its sleep EEG data.