Abstract

We investigated the electroencephalogram (EEG) activity in normal and epileptic subjects using two quantifiers adapted from nonlinear dynamics and deterministic chaos theory: (1) central tendency measure (CTM) and (2) Higuchi fractal dimension (HFD). CTM quantifies degree of variability while HFD quantifies complexity. The application of such techniques is justified by ascertaining the presence of nonlinearity in the EEG time series through the use of surrogate test. Both, CTM and HFD are found to be very effective in discriminating seizure state from seizure-free state of epileptic patients and in separating normal subjects (while performing motor/imagery tasks) from epileptic patients, with CTM outperforming HFD. HFD showed 92.5% (90.0%) sensitivity and 95.0% (100%) specificity in discriminating seizure-free EEG from seizures EEG (normal EEG while performing motor tasks from seizures EEG) while CTM showed 99.2% (98.0%) sensitivity and 95.9% (99.0%) specificity in discriminating seizure-free EEG from seizures EEG (normal EEG while performing motor tasks from seizures EEG). These results suggest that there exist differences in the ability to generate random time series between normal and epileptic subjects, and between seizure-free and seizure states, as estimated by CTM and HFD. We conclude that these quantifiers are most promising in providing new insight into the evolution of complexity of underlying brain electrical activity in different states.

Keywordscentral tendency measure, deterministic chaos, electroencephalogram, epilepsy, Higuchi fractal dimension, nonlinear dynamics

Author and Article Information

RecievedSep 30 2014 Accepted Dec 19 2014 PublishedJan 21 2015

Citation Kamath C (2015) Analysis of EEG signals in epileptic patients and control subjects using nonlinear deterministic chaotic and fractal quantifiers. Science Postprint 1(2): e00042. doi:10.14340/spp.2015.01A0003.

Copyright ©2014 The Authors. Science Postprint published by General Healthcare Inc. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs 2.1 Japan (CC BY-NC-ND 2.1 JP) License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

Funding No funders.

Competing interest: The author has no competing interests.

Ethics statement All the data is taken from publicly freely available for download and accessible EEG databases from (1) Department of Epileptology, Bonn University and (2) Physionet.org. The providers of these databases only demand a proper citing of the paper that describes these databases, which we have complied with. We used these databases so that it becomes easy to compare the results of our method with those of other researchers who have used the same databases.
Corresponding authorChandrakar Kamath

Address “ShanthaNilaya”, 107, 1st stage, Ananthnagar, Manipal-576104, India.

E-mail chandrakar.kamath@gmail.com

Peer reviewersD. Easwaramoorthy¹, Revewere A, B.

1. Affiliation: School of Advanced Sciences (Mathematics), VIT University, India

Introduction

Epilepsy is the most common brain disorder that is character-ized by intermittent abnormal firing of neurons in the brain which manifests as recurrent unprovoked seizures in the electroencephalogram (EEG) signal 1. The symptoms depend upon the location of the seizure onset in the brain and the spread of the seizure. A person is diagnosed as epileptic, only if there is a recurrent tendency to have seizures. Epileptic brain, like other nonlinear chaotic system repeatedly, but intermittently, makes abrupt transitions into and out of the seizure state 2-7. The behavior has been attributed to the epileptogenic focus driving the brain into self-organizing phase transitions from chaos to order. The EEG correlate of such a transition is characterized by the abrupt appearance of well organized oscillations out of the ongoing background activity, a characteristic of a seizure. In the waking state, the ongoing background activity is described by low to medium amplitude electrical activity, dominated by irregular waveforms. During seizure, this background activity is replaced by rhythmic, higher amplitude, organized, and self-sustained signal. In fact, the reverse transition (resetting mechanism) from order to chaos is initiated by the occurrence of a seizure 8. This dualism of chaos and order is the key feature of nonlinear dynamics. Since EEG signal may be considered chaotic, deterministic chaos theory seems to be promising to study the EEG dynamics and the underlying chaos in the brain 9. It has been shown in the literature that nonlinear analysis permits an improvised characterization and relative comparison of different physiological/pathological brain states 10, 11. More importantly, nonlinear methods may be applied to linear signals, but the vice-versa is not true. Various entropy and complexity algorithms have been developed and tried in the past two decades to characterize chaotic behavior in the EEG time series 12-17. Ghosh-Dastidar et al. used largest Lyapunov exponent, correlation dimension and standard deviation of the signal as features, after the application of wavelet analysis, to discriminate EEG signals of healthy and epileptic subjects 18. They achieved an accuracy of about 96.7% with their approach. Alam et al. employed largest Lyapunov exponent and correlation dimension in the empirical mode decomposition domain to discriminate EEG signals of healthy and epileptic subjects 19. Lempel-Ziv Complexity (LZC) has been used as a measure to quantify the regularity in the various epochs of epileptic seizure time series data 20. Hu et al. compared LZC with correlation entropy in the context of epileptic seizures detection from EEG data and discussed the advantages of the normalized LZC over the correlation entropy 21. Gao et al. employed several complexity measures, like LZC, permutation entropy, Lyapunov exponent, Kolmogorov entropy, and correlation dimension for characterizing EEG and to indicate epileptic seizures 22. Some of these, like, Lyapunov exponent and correlation dimension require large amount of data to arrive at meaningful results. Also, some require the data to be stationary over the duration of computation of the parameter, which cannot be realized easily with the practical biological data.

Further, the recent advances in the area of chaos and nonlinear dynamics analysis have revealed chaotic properties and fractal patterns of natural phenomena leading to a shift in the paradigm. A strong reason for this paradigm shift is the acknowledgement of variability in natural systems to be healthy, and in fact, is desirable. Such variability allows accommodating external perturbations 23. In this paper, we investigate the efficacy of two complexity measures from nonlinear dynamics and deterministic chaos theory, namely, central tendency measure (CTM), and Higuchi fractal dimension (HFD). The motivation to choose these measures is that they can be used to characterize a signal irrespective of the nature of the underlying dynamics, i.e. whether the signal is chaotic, deterministic, or stochastic 24, 25. This property is important to track phase transitions from chaos to order and from order to chaos. Further, these measures from the chaos theory can reveal subtle changes in brain waves and be related to higher cognitive processes. Moreover, they provide higher sensitivity and specificity for classification at low computational cost. We use these two, above mentioned, biomarkers and evaluate their potential in distinguishing among EEGs of different physio-pathological states, namely, (i) awake and relaxed with eyes open (AR-EO) EEG signals from normal subjects, (ii) normal EEG signals while performing motor/imagery tasks, and (iii) seizure-free (interictal) and seizure (ictal) signals from epileptic patients.

The prime advantages of CTM are the following: (1) It can be used to characterize a signal irrespective of the nature of the underlying dynamics, i.e. whether the signal is chaotic, deterministic, or stochastic; (2) With this approach, rather than defining a time series as chaotic or not chaotic, the degree of variability or chaos is evaluated; (3) We hypothesize that variability measures of EEG allow distinguishing between EEGs of different physio-pathological states and CTM is an appropriate measure for this; (4) CTM is a nonlinear approach to continuous chaotic modeling and hence the results can be used directly as a marker for the disease 26; (5) CTM is a simple parameter to estimate the signal variability with a low computation cost; (6) CTM measures the degree of variability accounting for events with significant slopes while neglecting all jitter (or noise) under a given level; (7) The CTM related second-order difference plot also provides a visual display of how closely the points are clustered around the center. The order of dispersion could help physicians to improve preliminary diagnosis. Scatter plots of first differences and CTM have been used in modeling biological systems, such as, hemodynamics and heart rate variability 27. CTM has also been tried for the analyses of EO and EC EEG signals using second-order difference plots 28. Their study used EEG signals from 33 able-bodied and 17 spinal cord injured participants. They demonstrated a statistical difference between CTMs of EO and EC states, but no quality parameters were evaluated. CTM has also been applied on intrinsic mode functions derived from seizure (ictal) and non-seizure (interictal) EEG signals to quantify their variability 29. They showed a statistical difference between CTMs of rhythms of seizure and seizure-free EEG signals (p <0.01). Their study did not compare normal state EEG with either interictal or ictal state EEG. The first shortfall of all the above CTM methods is that they used a long procedure to arrive at the optimum radius, which is critical in determining the outcome of CTM, since no guidelines are available to find optimum radius. We circumvented this problem employing a simple index, radial distance index (RDI), outlined in Section “Curvilinearity and Radial Distance Index (RDI)” in Materials and Methods. The second drawback is that most of them did not justify the application of their method through surrogate data analysis. The third shortfall is that none of the above CTM methods evaluated the diagnostic performance parameters (sensitivity and specificity). Such parameters are important not only for diagnosis, but also in the application of the proposed method in automated classification systems to separate different EEG states or implantable devices that serve as pacemakers or defibrillators for the brain. The present study is intended to go more deeply into the usefulness of CTM/RDI as a biomarker in distinguishing among EEGs of different physio-pathological states. The application of such a method is also justified by establishing the presence of nonlinearity in the electroencephalogram time series through the use of surrogate test.

When chaos theory has been applied to EEG signals, the dimensional complexity has mainly been used to investigate the changes in the dynamical behavior of the brain 30. Two approaches are available to compute the dynamical complexity. The first is by computing the correlation dimension by reconstructing the attractor in a multi-dimensional phase-space and the second is by computing fractal dimension of the signal in time domain. While the former requires a large quantity of stationary data and reconstruction of attractor in phase space, the latter can be achieved using HFD, a simple, fast, and reliable way to compute fractal dimension the of EEG signal directly in the time domain 31, 32. HFD quantifies complexity and self-similarity of the EEG signal. The main advantages of HFD are the following: (1) It can be used to characterize a signal irrespective of the nature of the underlying dynamics, i.e. whether the signal is chaotic, deterministic, or stochastic 31, 32; (2) Computation of appropriate FD requires, in general, stationary signals. Since computation of correlation dimension demands large number of samples it is difficult to achieve stationarity condition with long intervals of time series. Nevertheless, HFD does not require long duration epochs. Segment lengths of the order of 100 samples provide good time resolution; (3) It is less sensitive to noise 30; (4) As the reconstruction of the attractor phase space is not necessary to implement HFD, the algorithm is simple and fast so that it can be readily applied to real time signals in time domain; (5) HFD effectively compresses long signal epochs into short length epochs without losing diagnostically significant information so that instead of observing the original long record page by page, the physician can get a quick overview at once 32; (6) The running fractal dimension can be used to characterize short-lasting phenomena, while the mean fractal dimension can be used as meaningful feature 25. Klonowski has shown the efficacy of HFD in biomedicine applications, in particular, analysis of posturographic signals, evoked EEG, monitoring depth of sedation 25. Sourina O et al. used HFD to recognize emotions from EEG 32. Accardo et al. 33 and Klonowski 31, 34 found the suitability of HFD in identifying EEGs in different physio-pathological states such as open and closed eyes in wakeful conditions or during epileptic seizures. However, they had not evaluated the performance parameters (sensitivity and specificity). Such parameters, as mentioned above, are important not only for diagnosis, but also in the application of the proposed method in automatic classification systems or implanted devices that serve as pacemakers or defibrillators for the brain. Also, most of them did not justify the application of their method through surrogate data analysis. The aim of this study is to evaluate the performance parameters in distinguishing different physio-pathological states, hypothesizing that CTM and HFD change depending on the EEG signal. The application of such methods is also justified by ascertaining the presence of nonlinearity in the electroencephalogram time series through the use of surrogate test.

Materials and Methods

EEG data

The EEG data used is from three different groups: normal, seizure-free and seizure, provided by University of Bonn EEG database which is available in public domain 35. The database consists of five sets (designated A–E) each containing 100 single channel EEG segments of 23.6 s duration. These segments have been picked from continuous multi-channel EEG recordings after removal of any artifacts, like, muscle activity, or eye movements, making sure that they fulfilled stationarity requirements. Each segment is treated as a separate EEG signal so that in all there are 500 EEGs. This large data set improves statistical significance during comparison of results. Sets A and B contain segments taken from surface EEG recordings in awaken state acquired from five healthy volunteers using a standard 10–20 electrode placement scheme. The subjects were awake and relaxed with their eyes open (AR-EO) for set A and eyes closed (AR-EC) for set B, respectively. The segments for sets C, D, and E were acquired from five epileptic patients undergoing presurgical diagnosis. The diagnosis was temporal lobe epilepsy (epileptogenic focus: hippocampal formation). Sets C and D contained only activity measured during seizure free intervals (interictal epileptiform activity), with segments in set C recorded from hippocampal formation of the opposite hemisphere of the brain and those in set D recorded within epileptogenic zone. On the other hand, set E contained only seizure activity (ictal intervals), with all segments recorded from sites exhibiting ictal activity. All the EEG signals were recorded using the same 128-channel amplifier system using an average common reference. The data were digitized at 173.6 samples per sec. with 12 bit resolution. The bandpass filter setting was at 0.53–40 Hz (12 dB/octave). Each dataset has 4096 samples. In this work, four sets A, B (healthy subjects) and C, E (interictal epileptiform and ictal activity of epileptic subjects), each with 100 single channel EEG segments of 23.6 sec. duration, have been selected for CTM and HFD analysis. The dataset A comprising AR-EO EEG signals from healthy subjects, we designate as Normal Group-1.

To show the efficacy of the proposed nonlinear measures, EEG motor movement/imagery dataset from normal subjects performing a series of motor/imagery tasks, is also included for comparison 36, 37. Each record contains 64 channels recorded using the BCI2000 system and a set of task annotations. Only a few randomly chosen records from this database have been used for comparison. This database we designate as Normal Group-2.

Second-order difference plot and central tendency measure (CTM)

Quantifying signal variability employing CTM begins with second-order difference (SOD) plot. If s(n) represents a time series, a graph of [s(n + 2) - s(n + 1)] plotted against [s(n + 1) - s(n)] produces a scatter plot of first differences of the data, which many times, is called a second-order difference plot of s(n). Such a second-order difference plot, which is centered about the origin, can be useful tool to physicians, who can make a preliminary diagnosis by the visual inspection of these scatter diagrams. For example, consider Figure 1(a) for Normal Group-1 EEG and Figure 1(b) for Normal Group-2 EEG. Though some differences in the variability can be seen not much can be said about diagnosis by visual inspection. Now consider Figure 2, which shows the corresponding SOD plots for signals of Figures 1(a) and 1(b). In Figure 2(a) the dispersion of points is restricted more to quadrants one and three, while in Figure 2(b) the dispersion is spread to all the quadrants. This implies that though the SODs belong to EEGs of normal subjects, based on physiological states the SODs and the corresponding dispersion (or variability) is completely different suggesting a different diagnosis. Thus, SODs can help physicians to improve preliminary diagnoses. With this approach, rather than defining a time series as chaotic or not chaotic, the degree of variability or chaos is evaluated.

The central tendency measure (CTM) is used to quantify the degree of variability/chaos from the SOD scatter plot 26, 27. Given a circular region of radius r, about the origin, the CTM is defined as the ratio of the number of points that fall within this radius r, to the total number of points in the entire plot. This count involves the number of successive rates that include all sign combinations or quadrant regions within radius r. If N is the total number of points in the time series the scatter plot will have N - 2 points. CTM is computed, by selecting a circular region of radius r about the origin in the scatter plot and counting the number of points that fall within this region and normalizing this count by dividing it by the total number of points, as below. The radius r is selected based on the characteristic of the data.

CTM(r) represents a fraction of total count (N - 2), that lies within the circle of radius r. It has been shown that changes in frequency content and degrees of freedom are the prime reasons for any differences in CTM values between the two brain states 38.

Curvilinearity and Radial Distance Index (RDI)

The performance of this method highly depends upon the chosen value of r, which in turn depends on the character of the data. However, there are no specific rules to arrive at an optimum value of r. Instead of arriving at an optimum value of r, we adopt a different approach as below: We plot CTM for radius values from 0.1, 0.2,… 10 times standard deviation of the time series 39. It can be observed that the plot for CTM exhibits a kind of exponential rise from 0 to 1.0 (on a normalized scale) with increasing radial distance, r. The curvilinearity of this plot in the knee region can serve as a visual indicator of variability. The larger the curvilinearity, the smaller is the variability 40. To relate with the curvilinearity in the knee region and hence, variability of the series, we use a simple measure, called radial distance index, RDI, as the value of that radial distance when the CTM(r) just reaches 90% of the corresponding final value (on a normalized scale = 0.9) 40. The value of 90% (of the final value) is arrived at empirically so as to include the knee or curvature portion of the plot and to bring out a proper difference between subjects being compared; AR-EO EEG signals of healthy subjects or seizure-free and seizure EEG signals of epileptic subjects. It is found that the curvilinearity decreases with increasing value of RDI. In other words, the variability increases with increasing RDI. The degree of variability in EEGs of different physio-pathological states was found by comparing the mean values of RDI in respective states.

Higuchi’s fractal dimension (HFD)

The concept of FD was first introduced by Mandelbrot to study either temporal or spatial phenomenon that exhibit correlation over a range of scales. FD primarily quantifies scale independent complexity of a signal over time/space. It is known that the Euclidean dimension for a curve is 1 while that for a plane is 2. Practically, for example, for a sine wave this dimension is barely above 1 (a signal generated by a fully deterministic process) and for random noise this dimension is almost equal to 2 (a signal generated by fully a stochastic process). Any bio-signal, in general, has a dimension in between these two extremes. Treating EEG signal as a fractal curve, the FD can be used to characterize discrete dynamics of the signal. In fact, there is no need to ascertain if a given signal is really chaotic 25. Higuchi proposed a simple, accurate, and fast method to compute the fractal dimension of a signal 41. It represents the FD of a curve representing the amplitude of the signal under investigation on a plane as a function of time, which at the same time serves as a measure of complexity of the signal. It is to be noted that Higuchi’s fractal dimension (HFD) is completely different from the correlation dimension. If HFD is computed using a moving time window, then the resulting HFD characterizes the changes in signal complexity. Lower values of HFD imply dominance of lower frequencies in the spectrum, while higher values of HFD indicate the presence of higher frequencies in the spectrum of the signal. Abrupt changes in HFD represent nonstationarities in the signal 42.

Let x(1), x(2), …, x(N) represent a one dimensional discrete time sequence. The Higuchi fractal dimension is based on a measure of length L(k), of the curve that represents the original time series while using a segment of k samples as a unit of measure 41. If L(k) scales like

where, D is called fractal dimension and is a measure of complexity and irregular characteristics of the time series. For a simple curve D = 1 and for a curve which fills out the entire plane D = 2.

From the original time series: x(1), x(2), …, x(N), k number of new time series are constructed according to Higuchi’s algorithm. The new time series is given by

where, N is the total number of samples, m is the initial time, k is the interval time, and int(r) is the integer part of a real number r. For each time series x_m^k, the absolute differences between each two successive data points are summed to calculate the vertical length L_m(k) of the signal measured with the scale size k, starting from the m^th data point as given by

It is to be noted that L_m(k) is not a distance in the Euclidean sense. Instead it represents the normalized sum of absolute values of difference in ordinates pair of points spaced distant k (with the initial point m). Then the length of the EEG segment L(k) for the time delay k, is computed as the average of the k number of L_m(k) values, as defined by

Finally L(k) is plotted against 1/k on log-log scale for k = 1, 2, …, k_max. The angular coefficient of the linear regression of the plot {ln(L(k)), ln(1/k)} constitutes an estimation of the fractal dimension D of the given time series.

The best value for the maximum number of time series constructed, k_max, is to be found by the following procedure. As mentioned above, for each constructed time series the curve length is computed and plotted against the corresponding 1/k value on a log-log scale. The slope of the resulting curve represents the fractal dimension of the given time series. Next the fractal dimensions are plotted for different values of k_max. The point at which the fractal dimension plateaus corresponds to the best value of k_max beyond which the k_max value does not yield better results.

Surrogate data analysis

While the assumption of nonlinearity for many bio-signals may be acceptable in principle, it is necessary to show explicitly that application of such nonlinear tools is justified. Given a time series with limited number of samples and finite precision, it may not be possible to discriminate between nonlinear dynamics and linear dynamics involving stochastic components. Hence, there is a need to test for nonlinearity of the time series and we employ method of surrogate data for the statistical examination of the nonlinear properties underlying the original data. In this approach, the original data is compared with its surrogates which are derived from the original preserving most of its properties (like mean, variance, and power spectrum) 43. If the data under test is generated by a nonlinear process, the value for the computed statistic would be significantly different from that of the surrogate data, and the null hypothesis that a linear method characterizes the data can be rejected. On the other hand, the surrogate data sets and original data yield almost similar values for the statistics (within the standard deviation of the surrogate data sets) then the null hypothesis that the original data is random noise cannot be ruled out. The method of computing surrogate data sets from the original time series is as follows: First the Fourier transform of the original time series is computed, then the phases are randomized, and the inverse Fourier transform is evaluated. The resulting time series is that of the surrogate data. More details can be found in 43.

Statistical tests and box plots

To evaluate the statistical significance of the two nonlinear analysis methods, CTM/RDI and HFD, we use Kruskal-Wallis tests on the results of following comparisons: (i) awake and relaxed (AR-EO) EEG signals from normal subjects compared with seizure-free (interictal) and seizure (ictal) signals from epileptic patients; (ii) normal EEG signals while performing motor/imagery tasks compared with seizure-free (interictal) and seizure (ictal) signals from epileptic patients. When Kruskal-Wallis test rejects null hypothesis at least one sample is statistically significant from at least one other sample. However, the test does not identify in which sample this statistical significance occurs or for how many pairs of groups statistical significance is obtained. The Kruskal-Wallis test is a nonparametric version of one-way analysis of variance. The prime advantage is that the results are often not severely affected by changes in a small portion of the data. It evaluates whether the population medians on a dependent variable are the same across all levels of a factor. The assumption behind this test is that the measurements come from a continuous distribution, but not necessarily a normal distribution. The test is based on an analysis of variance using the ranks of the data values, not the data values. Parameters are regarded as statistically significant if p <0.01. The Kruskal-Wallis test is usually accompanied by a box and whisker plot which provides a visual summary of the data being analyzed.

Receiver operating characteristic (ROC) Analysis

Firstly, Kruskal-Wallis tests are used to evaluate the statistical differences between the CTM and HFD values of EEG series in healthy subjects (AR-EO) and in epileptic patients (seizure-free and seizure periods). If significant differences between physio-pathological states are found, then the potential of the nonlinear analysis method to discriminate between EEG series in healthy subjects or/and epileptic patients is evaluated using receiver operating characteristic (ROC) plots in terms of area under ROC (AROC) and the following performance parameters: sensitivity, specificity, precision, and accuracy. ROC plots are used to gauge the predictive ability of a classifier over a wide range of values 44. A threshold value is applied such that a feature value below this threshold will be assigned one category while a feature value above the threshold will be assigned other category. ROC curves are obtained by plotting sensitivity values (that represent the proportion of the features of category-1 and test positive) along the y axis against the corresponding (1-specificity) values (which represent the proportion of the correctly identified features of the category-2) for all the available cutoff points along the x axis. Accuracy is a related parameter that quantifies the total number of features (both categories 1 and 2) precisely classified. The AROC measures this discrimination, that is, the ability of the test to correctly classify into categories 1 and 2, and is regarded as an index of diagnostic accuracy. The optimum threshold is the cut-off point in which the highest accuracy (minimal false negative and false positive results) is obtained. This can be determined from the ROC curve as the closet value to the left top point (corresponding to 100% sensitivity and 100% specificity). AROC is a single number summary of performance, the larger the value, the better is the diagnostic test. An AROC value of 0.5 indicates that the test results are better than those obtained by chance, where as a value of 1.0 indicates a perfectly sensitive and specific test. A rough guide to classify the precision of a diagnostic test based on AROC is as follows: If the AROC is between 0.9 and 1.0, then the results are treated to be excellent; If the AROC is between 0.8 and 0.89, then the results are treated to be good; the results are fair for values between 0.7 and 0.79; the results are poor for values between 0.6 and 0.69; If the AROC is between 0.5 and 0.59, then the outcome is treated to be bad.

Results and Discussion

We evaluated the ability of the CTM and HFD to distinguish EEGs of different physio-pathological states, namely, (i) discriminate seizure-free (interictal) and seizure (ictal) states of epileptic group, (ii) discriminate Normal Group-1 in AR-EO state from seizure-free (interictal) and seizure (ictal) states of epileptic group, and (iii) discriminate Normal Group-2 while performing a motor/imagery tasks from seizure-free (interictal) and seizure (ictal) states of epileptic group. We averaged the RDI and HFD values for each state in the different groups. The results of statistical analysis of Kruskal-Wallis tests for both healthy and epileptic groups are shown in the respective box and whisker plots. The ROC analysis was then performed to assess the diagnostic ability of each nonlinear method.

Figures 1(a), 1(b), 1(c), and 1(d) show respectively, represent-ative EEG signals from Normal Group-1 in AR-EO state, Normal Group-2 while performing a motor task, epileptic group in seizure-free (interictal) state, and epileptic group in seizure (ictal) state. Figures 2(a), 2(b), 2(c), and 2(d) show the respective representative SOD plots for the same EEG signals depicted in Figure 1, in the same order. The plots clearly show differences between the different states. The points corresponding to AR-EO state tend to be located close around the origin, while those corresponding to motor task state are widely spread in the diagram. This implies that the EEG variability in brain activity is increased when performing a motor/imagery task compared to that when eyes are open in the normal subjects. This can be attributed to changes in frequency content and degrees of freedom in the corresponding EEG signals as the mental state is switched from one to another. The plots in Figures 2(c) and 2(d) clearly show differences between the two states of brain activity in an epileptic patient. It is found that the EEG during seizure intervals is characterized by greatly increased variability while during interictal intervals it is characterized by considerably decreased variability, even compared to those of normal subjects in Figures 2(a) and 2(b). This again can be attributed to changes in frequency content and/or degrees of freedom in the corresponding EEG signals as the brain state is switched from one to another. As mentioned in the methodology section, CTM are computed for several values of radii and the RDI is evaluated. The averaged CTM plots for different comparisons are shown in Figures 3(a) and 3(b). Figure 3(a) shows the averaged CTM plots for performing (i) discriminate seizure-free (interictal) and seizure (ictal) states of epileptic group, and (ii) discriminate Normal Group-1 in AR-EO state from seizure-free (interictal) and seizure (ictal) states of epileptic group. Figure 3(b) shows the averaged CTM plots for performing (iii) discriminate Normal Group-2 while performing a motor/imagery tasks from seizure-free (interictal) and seizure (ictal) states of epileptic group. The corresponding distributions of RDI values using box-whisker plots (without outliers) for the EEG signals from Normal Group-1 in AR-EO state and epileptic group in seizure-free (interictal) and seizure (ictal) states (comparison-1) is depicted in Figure 4(a), and for the EEG signals from Normal Group-2 while performing motor/imagery tasks and epileptic group in seizure-free (interictal) and seizure (ictal) states (comparison-2) is depicted in Figure 5(a), together with the p-value from Kruskal-Wallis tests. As mentioned above, the box plots provide a graphical visual summary of the data analyzed. Higher RDI values (141.50 ± 19.42) for seizure state compared to lower RDI (21.80 ± 4.57) for seizure-free state implies that variability in EEG is significantly enhanced when seizure occurs. The Normal Group-1 shows an intermediate RDI range (34.80 ± 4.21) and so also Normal Group-2 (44.61 ± 1.27). The values are expressed as mean ± SD. From the plots in Figure 4(a) or Figure 5(a) it is clear that the boxes for the different states do not overlap. However, the box of Normal Group-1/Normal Group-2 overlaps with whisker of the seizure-free state and vice versa. Nevertheless, the smallest p-value (p = 0) of Kruskal-Wallis test, in either case, shows that the RDI values for the physio-pathological states are significantly different. This implies that CTM method is able to separate well the EEGs of different physio-pathological states.

Next, we evaluated the diagnostic performance parameters of CTM method using ROC plots. The ROC curves of RDI for seizure-free and seizure states of epileptic subjects in Figure 4(b). The parameters at the cutoff points in the respective cases are shown in Table 1. In Figure 4(b) the AROC = 0.9984. This indicates that a randomly selected subject in seizure-free (interictal) state from the epileptic group has an RDI value smaller than that of a randomly chosen subject with seizure (ictal) state from the epileptic group in 98.40% of the time. In Figure 4(c) the AROC are 0.7898 and 0.9994, respectively, for separating Normal Group-1 in AR-EO state from seizure-free (interictal) and seizure (ictal) states of epileptic group. In Figure 5(b) the AROC are 0.9050 and 0.9996, respectively, for separating Normal Group-2 while performing motor/imagery tasks from seizure-free (interictal) and seizure (ictal) states of epileptic group. From the higher values obtained for the diagnostic parameters, as seen from Table 1, it is evident that in all the five cases the CTM has a potential to separate one physio-pathological state from the other.

Before proceeding further, we assess the appropriateness of applying the above nonlinear techniques through surrogate data analysis. We generated an ensemble of 20 surrogates for each EEG signal, as mentioned in section “Surrogate data analysis” in Materials and Methods. RDI and HFD were computed for each surrogate and averaged to obtain the corresponding mean RDI and HFD. We compared the RDI and HFD results from the original data with the mean RDI and HFD results from the ensemble surrogate data using paired t-tests. The statistical results are shown in Table 2, which reveal significant differences between original data and their surrogate in each dataset. This implies that variation inherent in the EEG data was not due to random fluctuations, but instead is more likely due to the consequences of some deterministic process. These findings confirm the appropriateness of applying the above CTM and HFD nonlinear analysis to the EEG data.

Now we investigate the efficacy of HFD in distinguishing different physio-pathological states. Figure 6(a) shows representative plots of log(L(k)) versus log(1/k) for Normal Group-1 in AR-EO state, and epilepsy group in seizure-free and seizure states (comparison-1), from which HFD values are computed as mentioned in section “Higuchi’s fractal dimension (HFD)” in Materials and Methods. Figure 6(b) shows representative plots of log(L(k)) versus log(1/k) for Normal Group-2 while performing a motor/imagery task and epilepsy group in seizure-free and seizure states (comparison-2). The corresponding box plots for HFD values in different physio-pathological states pertaining to the two discriminations are depicted in Figures 7(a) and 8(a), respectively. As mentioned above, the box plots provide a graphical visual summary of the data analyzed. Higher HFD values (1.936 ± 0.009) for seizure (ictal) state compared to HFD (1.839 ± 0.015) for seizure-free (interictal) state implies that fractal dimensionality in EEG is enhanced when seizures occur. Intermediate values of HFD are obtained for Normal Group-1 (1.901 ± 0.012) and Normal Group-2 (1.898 ± 0.013). All the values are expressed as mean ± SD. From the plots in Figures 7(a) and 8(a) it can be observed that the boxes for different physio-pathological states do not overlap. However, the box of one state overlaps with whisker of the other. Nevertheless, the smallest p-value (p = 0) of Kruskal-Wallis test indicates that the HFD values for different physio-pathological states are significantly different. Further, higher values of HFD during ictal intervals correspond to dominance in fast waves (alpha and beta) and reduced slow waves (delta and theta) of EEG. Similarly, lower values of HFD during interictal intervals indicate an increase in slow frequencies (delta and theta waves) and a decrease in fast frequencies (alpha and beta waves) of EEG. This implies that HFD method is also able to separate the EEGs of different physio-pathological states.

Finally, we evaluate the different performance parameters of HFD method using ROC analysis. The ROC curves of HFD for seizure-free and seizure states of epileptic subjects are depicted in Figure 7(b). The parameters at the cutoff points in the respective cases are shown in Table 3. In Figure 7(b) the AROC = 0.9668 which means that a randomly selected epileptic subject in ictal state has an HFD value larger than that of a randomly chosen epileptic subject in interictal state in 90.80% of the time. In Figure 7(c), the AROC are 0.9037 and 0.7986, respectively, for separating Normal Group-1 in AR-EO state from seizure-free (interictal) and seizure (ictal) states of epileptic group, while in Figure 8(b) the AROC are 0.8872 and 0.8093, respectively, for separating Normal Group-2 while performing a motor/imagery task from seizure-free (interictal) and seizure (ictal) states of epileptic group. From the discussion it is very evident that in all the cases the HFD has a potential to separate one physio-pathological state from the other. Further, comparison of Tables 1 and 3 reveals that both the chaotic modeling quantifiers perform extremely well, with CTM outdoing the HFD.

Now we carry out a comparison of the results of our method and few other methods that have used the same database. Ghosh-Dastidar et al. used largest Lyapunov exponent, correlation dimension and standard deviation of the signal as features, after the application of wavelet analysis, to discriminate EEG signals of healthy and epileptic subjects and achieved an accuracy of about 96.7% 16. Entropies have been applied to discriminate normal and epileptic subjects. The ANFIS classifier used showed an accuracy of 90% in the classification 45. Najumnissa et al. employed wavelet coefficients and ANFIS to distinguish normal and epileptic subjects with an average sensitivity of 97%, specificity of 99% and an accuracy of 98% 46. For detecting normal and epileptic seizure signals, S.H. Lee used wavelet coefficients in two-dimensional phase space and extracted 24 features, which were classified by a neural network with weighted membership functions 47. The classifier performance showed a sensitivity of 98.28%, a specificity of 95.16%, and an accuracy of 96.66%. Recently, S.H. Lee et al. proposed a method based on minimum number of features extracted from wavelet transform of normal and epileptic EEG signals using fuzzy neural network with weighted membership functions 48. Their new classifier showed an improved sensitivity, specificity, and accuracy of 99.60%, 100.0%, and 99.83%, respectively. F.S. Bao et al. used 38 features with a probabilistic neural network classifier 49 and reached an accuracy of 99.5% in separating normal and interictal groups, an accuracy of 98.3% in distinguishing normal and ictal groups, and an accuracy of 96.7% in distinguishing interictal and ictal groups with their approach. Our method, not using any neural network classifier, using CTM/RDI (considering average of Normal Group-1 and Normal Group-2) achieved an average accuracy of 80.65% in separating normal and interictal groups, an average accuracy of 99.40% in separating normal and ictal groups and an accuracy of 98.4% in distinguishing interictal and ictal groups.

To summarize, the chief findings of this study are the following: (1) CTM, in epileptic patients, during seizure (ictal) state is considerably increased while that during seizure-free (interictal) state is significantly decreased relative to each other. This implies that variability in EEG significantly increases during ictal state and significantly decreases during interictal state. This can be attributed to changes in frequency content or degrees of freedom of the EEG signals with a change in brain activity; (2) Normal Group-1 and Normal Group-2 show an RDI range in between those of seizure-free and seizure. This implies that variability in normal EEG is in between those of seizure-free and seizure; (3) During ictal state HFD is increased compared to that during interictal state. This implies that fast waves (alpha and beta) are enhanced and slow waves (delta and theta) are reduced in EEG during ictal state; (4) Normal Group-1 and Normal Group-2 show an HFD range in between those of seizure-free and seizure; (5) Both the chaotic modeling quantifiers performed extremely well. Between CTM and HFD, CTM showed much better performance; (6) CTM and HFD are valuable not only for detection of seizure and epilepsy, but also in the application of the proposed method in automated seizure detection systems and implanted devices that serve as pacemakers or defibrillators for the brain.

Future efforts will be focused to extend the results to other physiological and pathological states.

Conclusions

There is strong evidence that the mechanisms generating EEG obey nonlinear deterministic laws and that these processes are chaotic. This study on chaotic analysis of EEG time series from healthy and epileptic subjects using CTM and/or HFD shows promise not only in the discrimination of different physiological and pathological brain states, but also in providing new insight into the evolution of complexity of underlying brain activity. Both the methods are characterized by simplicity and adaptive capability. These results obtained using well known publicly available EEG databases and without any advanced classifiers, are better than or comparable to other approaches employing more signal analysis, more number of features, and complicated classifiers. The modification of the proposed method in order to automatically detect highly suspicious segments (regardless of their length) into long time EEG records and classify them into different physio-pathological states is an aspect that will be addressed in a future communication.

Author Contributions

The author, who is also the corresponding author, is the sole contributor to this work.

References

1. Bethesda MD (2004) Seizures and epilepsy: Hope through research. National Institute of Neurological Disorders and Stroke (NINDS). Available from: http://www.ninds.nih. Gov/health and medical/pubs/seizures and epilepsy htr.htm. (cited Mar 2004).
2. Iasemidis LD, Sackellares JC, Zaveri HP, Williams WJ (1990) Phase space topography and the Lyapunov exponent of the electrocorticogram in partial seizures. Brain Topogr. 2(3): pp.187−201. doi: 10.1007/BF01140588.
3. Iasemidis LD, Shiau DS, Chaovalitwongse W, Sackellares JC, Pardalos PM, Principe JC et al. (2003) Adaptive epileptic seizure prediction system. IEEE Trans. Biomed. Eng. 50 (5): pp. 616−627. doi: 10.1109/TBME.2003.810689.
4. Iasemidis LD (1991) On the dynamics of the human brain in temporal lobe epilepsy [Doctoral Dissertation]. Ann Arbor: University of Michigan.
5. Iasemidis LD, Sackellares JC (1991) The evolution with time of the spatial distribution of the largest Lyapunov exponent on the human epileptic cortex. In: Duke D and Pritchard W (eds.) Measuring Chaos in the Human Brain. Singapore: World Scientific. pp. 49−82.
6. Iasemidis LD, Sackellares JC (1996) Review: Chaos theory and epilepsy. The Neuroscientist 2(2): pp. 118−126. doi: 10.1177/107385849600200213.
7. Iasemidis LD, Principe JC, Sackellares JC (1999) Measurement and quantification of spatiotemporal dynamics of human epileptogenic seizures. In: Akay M (ed.) Nonlinear biomedical signal processing in Medicine. IEEE Press. ISBN: 978-0780360129.
8. Sackellares JC, Iasemidis LD, Shiau DS, Gilmore RL, Roper SN (2000) Epilepsy–When chaos fails. In: Lehnertz K, Arnhold J, Grassberger P, Elger CE (eds.) Chaos in the brain? Singapore: World Scientific. pp. 112−133.
9. Fiasché M, Schliebs S, Nobili L (2012) Integrating neural networks and chaotic measurements for modeling epileptic brain. In: Villa AEP, Duch W, Érdi P, Masulli F, Palm G (des.) Artificial Neural Networks and Machine Learning – ICANN 2012. Proceedings of the 22nd International Conference on Artificial Neural Networks Part I, Lausanne, Switzerland, Sep 11–14 2012. Springer Berlin Heidelberg. pp. 653−660. doi: 10.1007/978-3-642-33269-2_82.
10. Schliebs S, Fiasché M, Kasabov N (2012) Constructing robust liquid state machines to process highly variable data streams. In: Villa AEP, Duch W, Érdi P, Masulli F, Palm G (eds.) Artificial Neural Networks and Machine Learning. Proceedings of the 22nd International Conference on Artificial Neural Networks (CANN 2012) Part I, Sep 11–14 2012. Lausanne, Switzerland: Springer Berlin Heidelberg. pp. 604–611. doi: 10.1007/978-3-642-33269-2_76.
11. Klonowski W, Jernajczyk W, Niedzielska K, Rydz A, Stepien R (1999) Quantitative measure of complexity of EEG signal dynamics. Acta Neurobiol. Exp. 59(4): pp. 315−322.
12. Theiler J, Rapp PE (1996) Re-examination of the evidence for low-dimensional, nonlinear structure in the human electroencephalogram. Electroenceph Clin. Neurophysiolol. 98(3): pp. 213–222. doi:10.1016/0013-4694(95)00240-5.
13. Jeong J (2004) EEG dynamics in patients with Alzheimer’s disease. Clin. Neurophysiol. 115(7): pp. 1490–1505. doi:10.1016/j.clinph.2004.01.001.
14. Shaw R (1981) Strange attractors, chaotic behavior, and information flow. Z Naturforsch. A 36a(1): pp. 80−112.
15. Li XL, Ouyang GX, Richards DA (2007) Predictability analysis of absence seizures with permutation entropy. Epilepsy Res. 77(1): pp. 70–74. doi:10.1016/j.eplepsyres.2007.08.002.
16. Kumar SP, Sriraam N, Prabhu BG (2008) Automated detection of epileptic seizures using wavelet entropy feature with recurrent neural network classifier. Proceedings of the 2008 IEEE Region 10 Conference of TeNCON, Nov 19–21 2008, Hyderabad, India. pp. 1−5. doi:10.1109/TENCON.2008.4766836.
17. Ocak H (2009) Automatic detection of epileptic seizures in EEG using discrete wavelet transform and approximate entropy. Expert Syst. Appl. 36(2): pp. 2027−2036. doi:10.1016/j.eswa.2007.12.065.
18. Ghosh-Dastidar S, Adeli H, Dadmehr N (2007) Mixed-band wavelet-chaos-neural network methodology for epilepsy and epileptic seizure detection. IEEE Trans. Biomed. Eng. 54(9): pp.1545–1551. doi: 10.1109/TBME.2007.891945.
19. Alam SMS, Bhuiyan MIH, Aurangozeb, Shahriar ST (2012) EEG signal discrimination using non-linear dynamics in the EMD domain. IJCEE 4(3): pp. 326−330. doi: 10.7763/IJCEE.2012.V4.505.
20. Radhakrishnan N, Gangadhar BN (1998) Estimating regularity in epileptic seizure time series data - A complexity measure approach. IEEE Eng. Med. Bio. Mag. 17(3): pp. 89−94. doi: 10.1109/51.677174.
21. Hu J, Gao J, Principe J (2006) Analysis of biomedical signals by the Lempel-Ziv complexity: the effect on finite data size. IEEE Trans. Biomed. Eng. 53(12): pp. 2606−2609. doi: 10.1109/TBME.2006.883825.
22. Gao J, Hu J, Tung W (2011) Complexity measures of brain wave dynamics. Cogn. Neurodyn. 5: pp. 171−182. doi: 10.1007/s11571-011-9151-3.
23. Doyle TLA, Dugan EL, Humphries B, Newton RU (2004) Discriminating between elderly and young using a fractal dimension analysis of centre of pressure. Int. J. Med. Sci. 1(1): pp. 11−20. doi:10.7150/ijms.1.11.
24. Hornero R, Abásolo D, Jimeno N, Sánchez CI, Poza J, Aboy M (2006) Variability, regularity, and complexity of time series generated by schizophrenic patients and control subjects. IEEE IEEE Trans. Biomed. Eng. 53(2): pp. 210−218. doi: 10.1109/TBME.2005.862547.
25. Klonowski W (2009) Everything you wanted to ask about EEG but were afraid to get the right answer. Nonlinear Biomed. Phys. 3(2): p. 2. doi: 10.1186/1753-4631-3-2.
26. Cohen ME, Hudson DL (2001) Inclusion of ECG and EEG analysis in neural network models. In: Engineering in Medicine and Biology Society 2001, Proceedings of the 23rd International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), Oct 25–28 2001, Istanbul, Turkey. vol. 2: pp. 1621−1624. doi: 10.1109/IEMBS.2001.1020524.
27. Cohen ME, Hudson DL, Deedwania PC (1996) Applying continuous chaotic modeling to cardiac signals. IEEE Eng. Med. Biol. 15(5): pp. 97−102. doi: 10.1109/51.537065.
28. Thuraisingham RA, Tran Y, Boord P, Craig A (2007) Analysis of eyes open, eye closed EEG signals using second-order difference plot. Med. Bio. Eng. Comp. 45(12): pp. 1243−1249. doi: 10.1007/s11517-007-0268-9.
29. Bajaj V, Pachori RB (2012) Separation of rhythms of EEG signals based on Hilbert-Huang transformation with application to seizure detection. In: Lee G, Howard D, Kang JJ, Ślęzak D (eds.) Convergence and hybrid information technology, Lecture notes in computer science. Proceedings of the 6th International Conference, ICHIT 2012, Aug 23–25 2012. Daejeon, Korea. vol. 7425: pp 493−500. doi: 10.1007/978-3-642-32645-5_62.
30. Klonowski W (2011) Application of new non-linear dynamics methods in biosignal analysis. In: Rivas-Echeverria, Allegaert K, Wainstein DE (eds.) Proceedings of the World Medical Conference, Sep 26–28 2011, Prague, CzechRepublic. ISBN: 978-1618040367.
31. Klonowski W (2007) From conformons to human brains: an informal overview of nonlinear dynamics and its applications in biomedicine. Nonlinear Biomed. Phys. 1(1): pp. 1−19. doi:10.1186/1753-4631-1-5.
32. Sourina O, Liu Y (2011) A fractal-based algorithm of emotion recognition from EEG using arousal-valence model. In: Proceedings of the International Conference on Bio-inspired Systems and Signal Processing (BIOSIGNALS 2011), Jan 26–29 2011, Rome, Italy. Rome: Springer. pp. 209−214.
33. Accardo A, Anito M, Carrozzi M, Bouquet F (1997) Use of the fractal dimension for the analysis of electroencephalographic time series. Biol. Cybern. 77(5): pp. 339−350. doi: 10.1007/s004220050394.
34. Klonowski W (2011) Personalized neurological diagnostics from biomedical physicist’s point of view. In: Rivas-Echeverria, Allegaert K, Wainstein DE (eds.) Proceedings of the World Medical Conference, Sep 26–28 2011, Prague, CzechRepublic. ISBN: 978-1618040367.
35. Andrzejak RG, Lehnertz K, Rieke C, Mormann F, David P, Elger CE (2001) 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 64(6):061907. doi:10.1103/PhysRevE.64.061907.
36. Schalk G, McFarland DJ, Hinterberger T, Birbaumer N, Wolpaw JR (2004) BCI2000: A general-purpose brain-computer interface (BCI) system. IEEE Trans. Biomed. Eng. 51(6): pp. 1034−1043. doi: 10.1109/TBME.2004.827072.
37. Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG et al. (2000) PhysioBank, PhysioToolkit, and PhysioNet: Components of a new research resource for complex physiologic signals. Circulation 101(23): pp. e215−e220. doi: 10.1161/01.CIR.101.23.e215.
38. Tong F (2003) Cognitive neuroscience: Primary visual cortex and visual awareness. Nat. Rev. Neurosci. 4(3): pp. 219–229. doi:10.1038/nrn1055.
39. Thuraisingham RA (2009) A classification system to detect congestive heart failure using second-order difference plot of R–R intervals. Cardiol. Res. Pract. 2009: pp. 1−7. doi: 10.4061/2009/807379.
40. Kamath C (2012) A new approach to detect congestive heart failure using Teager energy nonlinear scatter plot of R-R interval series. Med. Eng. Phys. 34(7): pp. 841−848. doi: 10.1016/j.medengphy.2011.09.026.
41. Higuchi T (1998) Approach to an irregular time series on the basis of the fractal theory. Physica D- Nonlinear phenomena 31(2): pp. 277−283. doi: 10.1016/0167-2789(88)90081-4.
42. Hinrikus H, Bachmann M, Karai D, Klonowski W, Lass J, Stepien P et al. (2011) Higuchi’s fractal dimension for analysis of the effect of external periodic stressor on electrical oscillations in the brain. Med. Biol. Eng. Comput. 49(5): pp. 585–591. doi: 10.1007/s11517-011-0768-5.
43. Theiler J, Eubank S, Longtin A, Galdrikian B, Farmer JD (1992) Testing for nonlinearity in time series: the method of surrogate data. Physica D-Nonlinear phenomena 58(1–4), pp. 77–94. doi: 10.1016/0167-2789(92)90102-S.
44. Zweig MH, Campbell G (1993) Receiver-operating characteristic (ROC) plots: a fundamental evaluation tool in clinical medicine. Clin. Chem. 39(4): pp. 561–577.
45. Kannathal N, Choo ML, Acharya UR, Sadasivan PK (2005) Entropies for detection of epilepsy in EEG. Comput. Meth. Prog. Bio. 80(3): pp. 187–194. doi:10.1016/j.cmpb.2005.06.012.
46. Najumnissa D, Rangaswamy TR (2012) Detection and classification of epileptic seizures using wavelet feature extraction and adaptive neuro-fuzzy inference system. Int. J. Comp. Eng. Res. 2(3): pp.756−761.
47. Lee S-H (2013) Detection of normal and epileptic seizure EEG signals using a fuzzy neural network. In: 3rd International Conference on IT and Intelligent Systems (ICITIS'2013), Dec 25–26 2013, Bangkok, Thailand. pp. 84−86.
48. Lee S-H, Lim JS (2014) Minimum feature selection for epileptic seizure classification using wavelet-based feature extraction and a fuzzy neural network. Appl. Math. Inf. Sci. 8(3): pp. 1295−1300. doi:10.12785/amis/080344.
49. Bao, FS, Lie DYC, Zhang Y (2008) A new approach to automated epileptic diagnosis using EEG and probabilistic neural network. In: Tools with Artificial Intelligence 2008, 20th IEEE International Conference, Nov 3−5 2008, Dayton, Ohio. vol. 2: pp. 482−486. doi: 10.1109/ICTAI.2008.99.