Real-time Detection of Precursors to Epileptic Seizures: Non-Linear Analysis of System Dynamics.

We propose a novel approach for detecting precursors to epileptic seizures in intracranial electroencephalograms (iEEG), which is based on the analysis of system dynamics. In the proposed scheme, the largest Lyapunov exponent of the discrete wavelet packet transform (DWPT) of the segmented EEG signals is considered as the discriminating features. Such features are processed by a support vector machine (SVM) classifier to identify whether the corresponding segment of the EEG signal contains a precursor to an epileptic seizure. When consecutive EEG segments contain such precursors, a decision is made that a precursor is in fact detected. The proposed scheme is applied to the Freiburg dataset, and the results show that seizure precursors are detected in a time frame that unlike other existing schemes is very much convenient to patients, with sensitivity of 100% and negligible false positive detection rates.

INTRODUCTION

Epilepsy is characterized by occasional, excessive, and irregular discharging of neurons in the brain, which can be detected by its manifestations, i.e., the seizure.[1] A seizure (ictal) is induced when a critical mass of nerve cells is recruited. The process can be modeled by a system that evolves over time, from a non-epileptic state to a pre-seizure state, to a seizure state, to a post-seizure state, and finally back to a non-epileptic state. The stay-time in each state is unknown, can fluctuate, and depends on various factors, including mental and physical conditions of the patient. Moreover, there is enough evidence that transitions from one state to the next is not abrupt,[1] meaning that during each state (e.g., the pre-seizure state), different (unknown) successive events would cause the system to evolve into the next state (e.g., the seizure state). As such, seizure prediction means detection of such events (called seizure precursors) during the pre-seizure state prior to the actual transition to the seizure state. Seizure treatment and control depends on timely prediction of seizures prior to their occurrences.

The basic concept in seizure prediction is that time series analysis on EEG recordings may be used for detecting seizure precursors, which in turn can be used to trigger appropriate actions for safety and well-being of the patient and/or administering proper treatment for prevention and/or control of the forthcoming seizure.[2] Depending on the type and modality of control measures or treatments, desirable time-frames for prediction of a forthcoming seizure may be different. If treatment is in the form of electrical stimulation, only a short interval (less than one minute) is useful, as compared to other forms of treatment or control that require relatively long advance warning (more than five minutes). Seizure prediction in a timely manner is particularly important for drug-resistant epileptic patients, either to avoid harmful situations by the patient, or to initiate appropriate treatments that would prevent the occurrence of the predicted seizure.[2]

Various approaches have been employed to detect seizure precursors in EEG recordings that include detecting changes in Lyapunov exponents (LEs) of the EEG signals,[34] studying spatio-temporal EEG dynamics,[5] using non-linear methods for intracranial EEG analysis,[6] utilizing wavelet-based methods for similarity analysis of consecutive segments of EEG signal with a reference segment of the same EEG signal,[78] combining energy and wavelet transform of the EEG signal,[9] using support vector machines (SVMs) and convolution networks on intracranial EEGs,[10] using an adaptive Weiner algorithm on local field potentials (EEG)[11] and using SVM classifier on the auto-regressive (AR) model of EEG signals.[2] The difficulty with the approaches in[3711] is that some seizures are predicted at such intervals from their occurrences that are either too short for taking precautions to avoid harmful situations, or are not useful for non-electrical treatments. In addition, the approaches in[68] fail to predict many seizures at the expense of fixing the maximum acceptable false predictions; hence eroding the patient's confidence on seizure predications. Moreover, the methods in[5911] suffer from too many false positive alarms. The methods for extracting features in[10] require excessive calculations, meaning that the computational cost would be prohibitive or the processing time would be too long, which is not suited for real-time seizure prediction. Although the proposed scheme in,[2] which is a major step forward, provides a high rate of true predictions and a low rate of false predictions, but as we will show in this paper, it is possible to predict seizures in a time frame that is more convenient for patients, with fewer false positive alarms and failures and with reasonable (not excessive) calculations.

Considering the fact that we are dealing with a system that evolves over time with transitions from one state to the next, whose particulars are patient dependent, our approach is based on studying system dynamics. In doing so, we consider the largest Lyapunov exponent (LLE) of the discrete wavelet packet transform (DWPT) of the segmented EEG signals for each individual as features for detecting seizure precursors. Such features are processed by a SVM classifier whose output is further processed to identify whether the corresponding segment of the EEG signal contains a precursor to an epileptic seizure. When consecutive EEG segments contain such precursors, a decision is made that a precursor is in fact detected. We will show that our chosen features can significantly improve seizure prediction as compared to utilizing general and universal features, without causing inconvenience to patients. We will also show that our scheme predicts epilepsy seizures in a timely and efficient manner from the standard intracranial electroencephalogram (iEEG) recorded invasively from the cortex of patients.[12]

This paper is organized as follows. In Section 2, we state the problem of seizure prediction and describe performance measures. In Section 3, our scheme is explained in detail, followed by the results and conclusions in Sections 4 and 5, respectively.

PROBLEM STATEMENT AND PERFORMANCE MEASURES

Figure 1 shows the time-line of various events that include a time-interval called inter-ictal, during which, no seizure is expected; a time-interval called pre-ictal, during which, an alarm should be set to indicate a seizure is forthcoming after an elapse of time called the seizure prediction horizon (SPH). SPH is followed by a time-interval called the seizure occurrence period (SOP), during which, a seizure is expected to occur. SOP begins when SPH ends and ends when the actual seizure occurs. The exact time of the seizure onset may vary within SOP. The ictal period begins when SOP ends and ends when the actual seizure ends, followed by a post-ictal period, during which, no seizure is expected.

Figure 1

Different intervals in an electroencephalogram signal pertaining to an epilepsy seizure

We wish to predict seizures in a timely manner that would be useful for triggering appropriate actions for safety and well-being of the patient and/or administrating proper treatment for prevention and/or control of seizure. Specifically, based on the survey that was presented in,[13] we wish to detect a precursor to a seizure at least 10 minutes prior to the onset of that seizure, i.e., the prediction latency, defined as tA≜ SPH + SOP, should be at least 5 minutes. However, a relatively longer prediction latency is more desirable, and in this paper, we set to achieve a minimum of 10 minutes (instead of 5 minutes) advance warning.

If a seizure precursor is detected during the pre-ictal period, but there is no seizure within 90 minutes (i.e., when tA is longer than 90 minutes), it is called a false positive alarm.[8] During inter-ictal periods, a predicted seizure is also a false positive alarm. If a seizure precursor is not detected within 10 minutes prior to the seizure's onset, i.e., if a seizure is missed, or if tA is less than 10 minutes, it is classified as a false negative alarm.[13] The number of false predictions (positive and negative) is an indication of the performance of the system for a given patient.

Predication sensitivity for a given patient, denoted by S, is the ratio of the total number of correctly predicated epileptic seizures to the total number of seizures that occurred for that patient. False positive rate, denoted by RFP, is defined as the number of false positive alarms per hour and is an indication of prediction specificity.[26811]

PROCESSING

Figure 2 shows the block diagram of our proposed system. EEG signals recorded over extensive periods that contain precursors to epileptic attacks are inherently non-stationary.[14] In order to avoid complications that are associated with non-stationary signals, we use a sliding window of 4 seconds (each containing 1024 data points) for each epoch, during which, the EEG is assumed to be stationary. Each epoch is processed individually. In what follows, we describe the functionality of each block in Figure 2.

Figure 2

Block diagram of the proposed seizure prediction system

Data Set and Preprocessing

We applied our scheme to invasive EEG recordings from 21 epilepsy patients suffering from medically intractable focal epilepsy. EEG recordings were made during pre-surgical epilepsy monitoring at the Epilepsy Center of the Freiburg University Hospital in Germany,[12] containing a total of 86 labeled seizures for all 21 patients and less than 2 hours of pre-ictal data for each seizure. Besides, for each patient, we used a 24 hours non-seizure recorded data from the database in.[12] EEG recordings for Nc = 6 channels, denoted by s(t) = [s1(t), s2(t), …, s, where N = 6 is number of channels, were sampled at 256 samples per second. To eliminate possible line noise and to preserve the available information in the recorded EEG to its maximum extent, a 50 Hz notch filter is used.

Feature Extraction: DWPT, SWE, and LEs

Each epoch of EEG recordings denoted by x(n) = [x1(n), x2(n), …, x(n)] carries information related to Nc = 6 channels, containing 6 × 1024 samples (corresponding to 4 seconds of data) and is passed through a chain of two blocks for extracting its features that can be used for detecting seizure precursors. The first block involves DWPT and Shannon wavelet entropy (SWE) and the second one involves estimating the LLE. The extracted features that delineate dynamical dissimilarities are then used to discriminate between inter-ictal and pre-ictal patterns. In the sequel, we describe each block.

DWPT and SWE

To measure variations in filtered EEG signal x(n) that may emanate from seizure precursors, we use DWPT. This is a powerful tool for characterizing signals, which decomposes x(n) into its spectral components with a suitable time-frequency resolution.

Wavelet packet decomposition at level i produces 2 wavelet packets, each corresponding to a node. Each wavelet packet at each node in level i is represented by a vector of its coefficients whose elements are C, where m and j are time-localization and scale-parameters, respectively. The energy of each vector is

Where N is the number of coefficients in the vector that corresponds to the wavelet, and ℤ+ is the set of positive real integers. The total energy for each time window m is

In the time interval m and scale level j, the normalized energy is

Obviously, for each m, we have

To obtain a measure of order/disorder in each time interval m, we calculate the SWE[15] defined by

DWPT provides an optimal time-resolution by subband analysis of the original signal, as it uses the lowest cost function in building the library functions.[16]

In,[17] it is shown that an epileptic EEG signals do not have significant frequency components above 30 Hz, and 5-level wavelet decomposition is sufficient for its analysis. Accordingly, we use the Daubechies mother wavelet of order 4 (db4) with 5-level decomposition for each EEG channel, which decomposes each epoch for each EEG channel into 5 unequal frequency bands between 0 and 128 Hz. For the resulting 62 vectors of wavelet coefficients for each epoch in each EEG channel and for the respective epoch itself, we obtain Shannon wavelet entropies that constitute a feature vector a of 63 values. In this manner, each epoch of 1024 samples is reduced to 63 values, which is a significant reduction in the number of points that represent each epoch.

Extracting Features: LLE

As stated earlier, since we are dealing with a system that evolves over time with transitions from one state to the next, our approach is based on studying system dynamics using the concept of state movements[51819] and LEs. In doing so, for each vector a, we obtain 3 LEs and take the largest one (LLE) that pertains to that vector. In what follows, we explain the rationale for our approach.

To investigate the dynamics of transitions between normal and epileptic states, brain is assumed to be a bi-stable system, with two stable states called ictal and inter-ictal that simultaneously co-exist for the same set of system's parameters.[20] These states correspond to the attractors of the underlying dynamical system.[18] After each transition, one or more parameter values may gradually change, which may facilitate or impede an upcoming transition. It is also assumed that transitions between states are relatively fast as compared to the times spent by the system in these states.[20]

The values of LEs for system parameters that describe the dynamics of a system reflect the degree of chaos in that system.[14] In our case, a reduction in the LLE pertaining to the SWE of each epoch as compared to that of the previous epoch is an indication that the system is becoming less chaotic, i.e., a seizure may be forthcoming. To estimate the LEs for a chaotic time series, we use an efficient method in[21] as described below and shown in Figure 3.

Figure 3

Flowchart for extracting largest Lyapunov exponent

Phase-space Reconstruction (Embedding Process)

For the kth epoch with N samples, the embedding process constructs a P × (N − p) matrix A whose ith column is a = [a, …, a] where the values of a, j = 0, …, N − 1 are the Shannon wavelet entropies for each epoch, as shown in Figure 4a and p is the dimension of the embedded space (the number of estimated LEs).

Figure 4

Largest Lyapunov exponents calculation for a time-limited spectrum-constrained chaotic time series (a) Phase space reconstruction and (b) Neighborhood matrix of displacement vectors and estimating the tangent maps

Tangent Map Estimation

To estimate the dynamics of a series of points in phase space around a specific attractor, a small sphere with radius centered at a is considered, and the nearest K vectors to a that are within this sphere are considered as neighbors to a. The neighborhood matrix denoted by B, as shown in Figure 4b, is

In each step τ, the vector a and its neighborhood a proceed to a and a, respectively. The neighborhood matrix B is estimated as a tangent map in a as M = B, where B is the transpose of B. The value of i is upper bounded to N – K.

QR decomposition

In each step τ, QR decomposition uses an upper triangular matrix R(τ) for a, containing information to describe how an orthogonal basis basis Q for a is modified by M(τ) to produce a new orthogonal basis Q for a so that

M Q= Q      (7)

where Q = I is an identity orthogonal matrix. The matrix M is also called the trajectory matrix. At each step, an upper triangular matrix R for a is formed, whose diagonal elements are estimates of p LEs.

Estimating LEs

Having estimated successive upper triangular matrices R at each step τ, we take the average of the estimated LEs by λ = limτ→ , where (R) is the kth diagonal element of R(τ). The number of steps for each epoch is upper bounded to N – p – K.

We use the Cao method described in[22] to determine the embedding dimension p, resulting in P = 3, and set the number of neighboring points K to 20. In this paper, for each EEG channel, 3 LEs (features) for each epoch is estimated. Noticeable and to some extent abrupt and prolonged changes in the steady state (not short-term transitions) in the feature space for successive epochs emanate from loss of coherence due to transitions from non-seizure or inter-ictal to ictal states. Note that in this way, a SWE feature vector of 63 components is mapped into a much smaller feature vector of 3 components for each epoch of sampled EEG data (1024 points). In order to further reduce the number of features, for each EEG channel, we only use the largest of the 3 LEs, denoted by y, for the kth epoch.

Figure 5a–d show typical variations in the largest estimated LE versus time for a randomly selected patient and EEG channel for more than 20 hours of inter-ictal data and for 3 pre-ictal periods, each leading to a seizure. Note that when a seizure is not forthcoming, e.g., Figure 5a, changes in the steady state value of the largest estimated LE are not prolonged. Note also that while we can observe prolonged changes that might be attributable to a forthcoming seizure, e.g., Figure 5b–d, such changes are not similar for all cases, and in some instances, e.g., Figure 5b, there is a need for further processing to denoise the trend. Specifically, Figure 5b–d show variations for the same patient, but for different seizures. Vertical dashed lines to the right of the last 3 latter figures show the onsets of the seizures as indicated in the database.[12] Note that the time from the onset of visible and prolonged changes in the steady state value of the largest estimated LE to the onset of the corresponding seizure is not the same for all seizures. Furthermore, the pattern of changes may be different even for the same patient for different seizures. It is evident that we need to have a pattern discriminating framework for further processing.

Figure 5

Calculated yk(t) for an epileptic electroencephalograms of a patient during (a) More than 20 hours of inter-ictal data and (b-d) The available pre-ictal data for 3 seizures

SVM

A support vector machine (SVM) is a machine-learning method that is based on the statistical learning theory, in which all training points in each class are equally treated. To distinguish between two disjoint feature vectors, the training feature vector u for each y is mapped into a target value l = {−1, +1} via a kernel function and the chosen soft margin between two disjoint feature vectors.[223] The performance of the SVM depends on the choice of its kernel function and its respective parameter value, as well as on the size of the above mentioned soft margin.

In this paper, as in,[2] our choice for the kernel is the Gaussian radial basis function K(u = e-γ||u, with γ > 0 and a soft margin C.[2] The values of γ and C are so chosen to yield optimal performance for the SVM. In doing so, the k-fold cross validation is performed to determine the values of γ and C. In k-fold cross validation, the training set is divided into k subsets of equal size, from which, one subset is sequentially used to validate the classifier's performance, and k − 1 remaining subsets are used for training.[2] As such, each subset is classified once. Next, to obtain optimal values of C and γ for each fold, a diadic coarse grid search in the parameter's space for C = [20, 21, …, 210] and γ = [20, 21, …, 24] (resulting in 55 trials) is performed and for each combination of C and γ, two performance measures of S and R are obtained. The values of C and γ that yield the highest value of S - R are considered as optimal for that fold. For a rigorous evaluation, the optimal values of C and γ for all folds are obtained, and then the corresponding SVM classifier is used to process the incoming test data in each run.

In a data set, some samples may be outliers, and some may be corrupted by noise. Such samples may negatively affect the SVM's performance. In order to alleviate this, in addition to the SVM's label, we also utilize the probability P (class|input) of the corresponding label. In what follows, we explain how to combine the SVM's outputs (labels and probability estimates) to obtain a measure of rapid changes in successive epochs, which can be used to detect a seizure precursor.

The dimension of feature vector fed to the SVM is the number of EEG channels Nc. This is because for each epoch in each EEG channel, we only consider one feature, namely the LLE, representing a negligible computational load compared to that of[2] for training/testing the classifier, making the approach suitable for real time use.

Post-processing

For each epoch of 4 seconds duration (indexed by k), a feature vector consisting of 6 components (the LLE for each channel) is fed to the SVM, for which the SVM's output denoted by b and its probability denoted by p are obtained. The value of b = −1 indicates that the respective epoch belongs to pre-ictal/ictal phase with p > 0.5, and b = +1 corresponds to an epoch not in the pre-ictal/ictal phase with p > 0.5.

For a window containing five consecutive 4-second epochs (a total of 20 seconds of recorded EEG signals for 6 channels), we obtain , where is the epoch index for each epoch in each window, and is the window index. The value of B shows the degree by which the population of epochs in each window belongs to a pre-ictal/ictal mode. In calculating the value of P, weighting the SVM output by its respective probability reduces the impact of outliers or noisy values due to the fact that their corresponding probabilities are smaller. Note that when the probability is either 0 or 1, the values of B and P are the same.

The above is repeated by indexing through k and κ, providing the values of B and P versus time, as shown in Figures 6 and 7 for the same patient as in Figure 5. As can be seen, the patterns for inter-ictal and pre-ictal periods are visibly different, and seem suitable for further processing as compared to the direct binary output of the SVM. Note that −5 < B < +5 and −5 < P < +0.5.

Figure 6

The calculated Bk(t) over successive windows of 5 samples for (a) Inter-ictal and (b) Pre-ictal states. Horizontal lines are in minutes

Figure 7

The calculated Pk(t) over successive windows of 5 samples for (a) Inter-ictal and (b) Pre-ictal states. Horizontal lines are in minutes

At the end of each window indexed by κ, two binary values of χ and χ that correspond to B and P, respectively, are set as follows. The value of χ is set to +1 when the majority of epochs in that window belong to the pre-ictal/ictal phase (i.e., when B < –2). Otherwise, χ is set to 0. Similarly, the value of χ is set to +1 when the majority of epochs in that window belong to the pre-ictal/ictal phase (i.e., when P < 0). Otherwise, χ is set to 0.

Next, the stay time in each state (+1 or 0) for χ and χ, denoted by ξ and ξ, respectively, are obtained for consecutive windows, as shown in Figures 8 and 9, respectively. Note the visible difference between the values of ξ and ξ for inter-ictal, and pre-ictal periods. Note also that a chattering behavior in ξ and ξ indicates that the nervous system is rapidly changing, which can be regarded as an indication that a seizure is forthcoming. However, care must be exercised in setting the alarm, as the nervous system may experience rapid changes as well when no seizure is forthcoming. We observed that variations in stay times in each state for ξ and ξ can be used in developing a simple algorithm for detecting seizure pre-cursers. In doing so, we experimented with different options and parameter values, and chose the one that can identify seizure pre-cursers reliably and in a time frame that is convenient to patients. Our proposed algorithm is described below.

Figure 8

The output ξ(t) of support vector machine for (a) Inter-ictal and (b) Pre-ictal states. Both horizontal and vertical lines are in minutes

Figure 9

The output ξ(t) of support vector machine for (a) Inter-ictal and (b) Pre-ictal states. Both horizontal and vertical lines are in minutes

If the stay time in state 1 for χ (i.e., the value ξ) exceeds 10 minutes, it is an indication that the nervous system is undergoing a noticeable synchronization period, i.e., a seizure is forthcoming, and an alarm is set to indicate that a seizure precursor is detected. If, however, in the meantime, 4 consecutive values of χ are not the same, we focus on stay times in state 1 for χ (i.e., the value of ξ) for the next 20 minutes. During this period of 20 minutes, we begin by looking at a sequence of four non-overlapping windows each containing 3 consecutive non-zero values of ξ. Only when the values of ξ are not the same in each of the four sequential windows, i.e., when the 3 consecutive non-zero values of ξ are not the same 4 times in a row, we take it that the nervous system is reliably evolving into a seizure state, and an alarm is set to indicate that a seizure precursor is detected. Otherwise, we look at the next sequence of four non-overlapping windows each containing 3 consecutive non-zero values of ξ, and repeat the above. This cycle continues until all non-zero values of ξ in each period of 20 minutes are taken into account.

EXPERIMENTAL RESULTS

The Freiburg EEG Database

We apply our scheme to the data set in[12] that includes a total of 86 seizures with 2-5 seizures for each patient and at least 50 minutes of pre-ictal data for each seizures, and 21 × 24 = 504 hours of inter-ictal data for 21 epileptic patients.

Using a Neurofile NT digital video system with 128 channels, EEG recordings were acquired at the sampling rate of 256 samples/s via a 16-bit analog-to-digital converter. Implanted grid, strip, and depth electrodes with six contacts were used, resulting in 6 recording channels.

For the same reason mentioned in,[2] we selected a subgroup of 9 patients from the 21 available ones, each having electrocorticography recordings and more than two seizures.

Numerical Results and Discussion

Numerical results are shown in Table 1. Note that by applying our proposed scheme, all seizure precursors were successfully detected (i.e., similar to,[2] sensitivity of 100% was achieved for all patients) but within the time frame that was very much convenient to patients (from 90 to 10 minutes prior to the onset of a seizure) as compared to the prediction time frame of 115 to 2 minutes prior to the onset of a seizure in.[2]

Table 1

Summary of results

Moreover, on the average, a false alarm rate of 0.08/h for the inter-ictal period (as compared to the false alarm rate of 0.16/h for the inter-ictal period in[2]) and a false alarm rate of 0.10/h for the whole (ictal and inter-ictal) period as compared to the false alarm rate of 0.36/h for the whole period in[2]) were obtained, meaning that in very few cases, other events were falsely detected as seizure precursors. Note that the results obtained in[2] are superior to other existing works (e.g.,[10]). Hence, the fact that our proposed scheme produces better results as compared to[2] indicates that our proposed scheme significantly outperforms other existing schemes. Processing time for implementing our proposed scheme to generate feature vectors for a 24 hours data segment using Matlab 7.0 on a 2.00 GHz Intel dual core CPU with 2 GB of RAM is below 20 minutes. Note that when feature vectors are obtained, incoming EEG signals are processed on-line, and seizure precursors are detected in real time.

It is evident that the proposed scheme yields improved performances for all patients, in terms of significantly lower false alarm rates and more convenient detection times.

CONCLUSIONS

We presented a novel approach for detecting epileptic seizure precursors that is based on the analysis of system dynamics obtained from intracranial electroencephalogram (iEEG). In doing so, we used the discrete wavelet packet transform (DWPT) and the Shannon wavelet entropy (SWE), together with the largest Lyapunov exponent (LLE), and non-linear dynamics to obtain a set of patient-specific discriminating features. We showed that such features can be utilized to detect epileptic seizure precursors in a timely and efficient manner. The sensitivity of 100% and a negligible false positive rate were achieved. The results are significantly better than those of existing methods.

BIOGRAPHIES

Sahar Nesaei received her B.Sc. and M.Sc. degrees from Shahaid Beheshti University and Tarbiat Modares University, Tehran, Iran in 2002 and 2005, respectively both in Biomedical Engineering. She is now a PhD candidate in the Department of Electrical and Computer Engineering, Tarbiat Modares University, Tehran, Iran. Her current research interests are applications of nonlinear dynamics and time series statistical analysis for understanding the brain's state-transitions in human epileptic patients for predicting epileptic seizures.

E-mail: s_nesaei@modares.ac.ir

Ahmad R. Sharafat is a professor of Electrical and Computer Engineering at Tarbiat Modares University, Tehran, Iran. He received his B.Sc. degree from Sharif University of Technology, Tehran, Iran, and his M.Sc. and his Ph.D. degrees both from Stanford University, Stanford, California, all in Electrical Engineering in 1975, 1976, and 1981, respectively. His research interests are advanced signal processing techniques, and communications systems and networks. He is a Senior Member of the IEEE and Sigma Xi.

E-mail: sharafat@modares.ac.ir