Although treatment for epilepsy is available and effective for nearly 70 percent of patients, many remain in need of new therapeutic approaches. Predicting the impending seizures in these patients could significantly enhance their quality of life if the prediction performance is clinically practical. In this study, we investigate the improvement of the performance of a seizure prediction algorithm in 17 patients with mesial temporal lobe epilepsy by means of a novel measure. Scale-free dynamics of the intracerebral EEG are quantified through robust estimates of the scaling exponents—the first cumulants—derived from a wavelet leader and bootstrap based multifractal analysis. The cumulants are investigated for the discriminability between preictal and interictal epochs. The performance of our recently published patient-specific seizure prediction algorithm is then out-of-sample tested on long-lasting data using combinations of cumulants and state similarity measures previously introduced. By using the first cumulant in combination with state similarity measures, up to 13 of 17 patients had seizures predicted above chance with clinically practical levels of sensitivity (80.5%) and specificity (25.1% of total time under warning) for prediction horizons above 25 min. These results indicate that the scale-free dynamics of the preictal state are different from those of the interictal state. Quantifiers of these dynamics may carry a predictive power that can be used to improve seizure prediction performance.

Funding: This work was supported by the Canadian Institutes of Health Research (CIHR) grants MOP-10189/102710 ( http://www.cihr-irsc.gc.ca ), by the Royal Society of Canada and the Natural Sciences and Engineering Research Council of Canada (NSERC) grant CHRPJ 323490-06 ( http://www.nserc-crsng.gc.ca ), and by the joint NSERC/CIHR grant CHRP-CPG-80098. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability: Data is confidential and may be available from the Montreal Neurological Institute provided approval by the Research Ethics Board of the Montreal Neurological Institute and Hospital. Requests may be submitted to Jean Gotman ( jean.gotman@mcgill.ca ), Montreal Neurological Institute, 3801 University Street, room 767, H3A 2B4 Montreal.

Copyright: © 2015 Gadhoumi et al. This is an open access article distributed under the terms of the Creative Commons Attribution License , which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited

Using scaling analysis based descriptors, we investigate the discriminability of preictal and interictal EEG epochs in a group of patients with intractable mesial temporal lobe epilepsy. We then assess the prediction performance of our previously proposed patient-specific classifier using these descriptors individually and in combination with our previously introduced state similarity measures.

In this study, we examine whether combining measures of two a priori unrelated properties of the EEG, the so-called thermodynamic and scale invariance properties, would result in an improvement of the prediction performance. Descriptors of the thermodynamic property of the intracerebral EEG, namely the wavelet energy and entropy, have been introduced in our previous study [ 1 ] and estimated in the high frequency range (50–450 Hz). We demonstrated that measures derived from these descriptors were useful to achieve significant and practical seizure prediction performance results [ 17 ]. Scale invariance is a property related to temporal irregularity in signals, first introduced in the study of turbulence [ 18 ]. Evidence of scale invariance in the intracranial EEG and its significance to brain electrophysiological function has been demonstrated in [ 19 ]. Here, we analyze the scale invariance property of the intracerebral EEG via descriptors extracted through a novel scaling analysis framework: the bootstrap and wavelet leader based multifractal analysis [ 20 , 21 ].

Many studies in seizure prediction adopted the approach of combining EEG features that capture different properties of the signal [ 11 – 16 ]. Whether reported prediction performances are superior to those achieved with single feature algorithms is yet to be verified. It remains intuitive though that an increase in the prediction performance is expected when combining the predictive power of different methods in a complementary manner.

The existence of a preictal state distinguishable from an interictal state of the epileptic brain is supported by growing evidence of electrophysiological changes preceding the ictal phase [ 1 – 6 ]. Transitions from interictal to ictal state are likely to be governed by various mechanisms [ 7 – 9 ] with different electrophysiological manifestations. Revealing changes in dynamical properties of the electroencephalogram (EEG) related to state transition is probably a multifaceted problem and solving it requires a combination of tools each adapted to reveal changes in a distinct aspect of signal properties. For example, preictal changes are probably better detected using a combination of intra- and inter-channel EEG information rather than using either information alone [ 10 ].

Materials and Methods

Ethics statement This study was approved by the Montreal Neurological Institute and Hospital Research Ethics Board (REB). The REB acts in conformity with standards set forth by the Tri-Council Policy Statement and Ethical Conduct for Research Involving Humans (Canada), and by the Rules and Regulations of the Department of Health and Human Services and the Food and Drug Administration (US) governing human subjects research and functioning in a manner consistent with internationally accepted principles of good clinical practice. All adult subjects signed an REB-approved written informed consent form.

Materials To ensure an unbiased comparison between the performance of scaling analysis based features and originally proposed state similarity measures, we use the same dataset as in our previous work [17]. Long-lasting intracerebral EEG (iEEG) recordings, low-pass filtered at 500 Hz and sampled at 2000 Hz, from 17 consecutive patients who attended the Montreal Neurological Institute, Montreal, Canada, between 2004 and 2011 and who underwent depth electrode investigation of refractory epilepsy, were analyzed. The selected patients had a mesial temporal lobe epilepsy (MTLE) diagnosis. They had at least five seizures recorded at 2000 Hz during the investigation. Depending on the number of implanted electrodes, nine to 18 bipolar iEEG channels from the four deepest contacts of depth electrodes targeting the bilateral mesial temporal structures (amygdala, hippocampus and parahippocampus) were selected. A total of 1565 hours and 175 seizures were analyzed. Electrographic onsets of seizures were determined by an experienced neurologist. Only seizures that are at least two hours apart were considered in the analysis in order to minimize the postictal effect. The iEEG data was split into training and testing subsets for each patient. The training subset has two or three preictal epochs of uninterrupted iEEG recordings lasting between six and 22 min across patients, and five interictal epochs lasting almost one hour each and separated by at least one hour from each other and by four hours from any seizure. The training data were selected from the beginning of the patient’s iEEG recordings. The remaining iEEG data were used for testing. They are made of continuous multiday iEEG recordings. The number of test seizures per patient ranged between three and 24. In total, 119 hours of iEEG data were used in training and 1446 hours in testing (see S1 Table for a list of training and testing data duration and number of test seizures for each patient). Patient clinical characteristics, surgery outcomes and iEEG data are summarized in Table 1. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Table 1. Summary of iEEG dataset and seizure onset. https://doi.org/10.1371/journal.pone.0121182.t001

Theoretical background of scaling analysis Scale invariance or “scaling” is defined as the absence of a particular time scale playing a characteristic role in the process [22]. Such a process is called a “scale free” process. For stochastic processes such as in the case of EEG, scale invariance implies that the statistical properties at different time scales (e.g., hours versus minutes versus seconds) effectively remain the same [23]. Implicitly, scaling analysis identifies and characterizes the rules describing the relation between different scales. These rules are principally governed by power laws. Practical scaling analysis aims mostly at estimating the exponents that characterize these power laws: If X designates an EEG signal from one channel and T X (a, t) a multi-resolution measure of the content of X around a time t and a scale a, then the scale invariance property of X is described by the power law behaviors of the time average of the qth power of T X (a, t) with respect to the analysis scale a for a given (large) range of scales a ∈ [a m , a M ], a M /a m >>1. (1) with n a the number of samples of X at the scale a. For a range of order q values, the set of estimated exponents ζ(q) fully describes the statistical distribution of the signal and provides a powerful means of characterizing the scale invariance of the signal.

Models of scaling analysis Various processes exhibiting scale invariance properties are used as models for scaling analysis. The following summarizes some of the most studied models. - Self-similarity: the framework of “self-similarity” [24] is one of the first models that fulfills the scale invariance property in a mathematically simple and precise way [25]. A random process X(t) is said to be self-similar if: (2) Where designates equality of statistical properties. The parameter H, called Hurst exponent, defines and controls the scale invariance of the process X(t). - 1/f processes: for many real-data processes, scale invariance may exist only for a range of scales and/or the scaling exponents cannot be described with a single parameter. For such processes, self-similarity is not a suitable model. The so-called 1/f processes are more flexible models of scale invariance that address the limitations of the self-similarity model. These are processes for which the spectral density obeys a power law with a sufficiently large range of frequencies: (3) where α is a scalar referred to as the “index of long-range dependence”. It is related to the Hurst exponent by: (4) Fig 1 illustrates a preictal and an interictal iEEG epochs exhibiting a 1/f process. Two particular 1/f processes define two interesting models, namely long-range dependency [26] and monofractality. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 1. Power spectra of 1-min interictal (left) and preictal (right) epochs from 18 iEEG channels (patient 15). Power spectrum in red is the average power spectrum of all channels. Both epochs show a power law behavior (P(f) ∼ 1/fα) illustrated with the fitted dotted red line for the range of frequencies delimited with a dashed blue line (corresponding to the range of scales j between 3 and 7). The exponents α obtained by a least-square fit are different in the preictal and interictal epochs. https://doi.org/10.1371/journal.pone.0121182.g001 - Long-range-dependency (LRD): LRD is a model defined for scaling observed in the limit of small frequencies (equivalently, large scales), f m → 0 in Eq (3). It is usually defined in terms of covariance properties relevant to second-order stationary processes. In long-range dependent processes, the temporal correlation between values decays very slowly, allowing for a dependency between past and future values. - Fractals and multifractals: fractal processes are models of scale invariance described through the local (Hölder) regularity of the process sample path measured by the Hölder exponents h(t): (5) The local regularity is the degree of smoothness of the signal in a small time interval. Strong local regularity indicates a smooth signal and weak local regularity indicates less smoothness. When the Hölder exponents h(t) are a deterministic and stationary function, the process exhibits a slow and smooth regularity fluctuation over time and is called fractal. When the Hölder exponents h(t) are the same for all t, the process exhibits constant regularity along its sample paths; it is referred to as monofractal process. An example of monofractal process in the fractional Brownian motion (Fig 2.A). When the fluctuation of h(t) is random and highly undeterministic, the process is said to be multifractal. Instead of Hölder exponents, the irregularity of such processes is described rather through a spectrum D(h) called the singularity (multifractal) spectrum [27,28]. An example of a multifractal process is the multifractal random walk (Fig 2.B). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 2. Estimation of cumulants c 1 and c 2 using wavelet leader and bootstrap based scaling analysis in signals with different scale invariance properties. Each signal is 4096 samples. 100 bootstrap wavelet leader resamples were used in cumulant estimation. (a) Realization of a self-similar signal (fractional Brownian motion). ζ(q) is linear in q. c 1 ≠ 0 and c 2 ≈ 0. (b) Realization of a multifractal random walk. ζ(q) is nonlinear in q. c 1 , c 2 ≠ 0. (c) EEG channel recording showing nonlinear relation of ζ(q) in q. c 1 , c 2 ≠ 0. (Each column, top to bottom: Signal plot, regression plot of ζ(q) exponent estimates and boxplot of cumulant estimates). https://doi.org/10.1371/journal.pone.0121182.g002

Wavelet leader based scaling analysis Wendt et al. [21] proposed a scaling analysis framework for a robust estimation of the scaling exponents using cumulants as surrogate measures of the exponents. In the proposed framework, wavelet leaders [29]—quantities derived from the discrete wavelet analysis of the signal—are used as a multi-resolution measure in Eq (1). Using a mathematical formalism (see S1 appendix), it was established that the scaling exponents ζ(q) can be estimated through the cumulants c p of the Taylor series polynomial expansion: (6) The advantage of using the cumulants as surrogate measures of the scaling exponents lies in the fact that they emphasize the departure from linear behavior of ζ(q) in q: there exists p ≥ 2: c p ≠ 0 [30]. From a single path (single observation) of the signal, the cumulants c p can be estimated through statistical models and bootstrap approaches. At each scale, R bootstrap resamples are generated from the wavelet leader coefficients. R bootstrap cumulant estimates are calculated and then averaged to obtain the final estimates. In practice, the first few cumulants are sufficient to gather most of the scaling property from the signal [21,31]. For c 1 ≠ 0, knowing whether c 2 ≠ 0 is practically equivalent to choosing between monofractal and multifractal models of scale invariance. In this study, we use the first and second cumulants, c 1 and c 2 . Fig 2 shows examples of cumulant estimates from realizations of simulated signals with different scale invariance properties and from a bipolar recording of an intracranial EEG channel.

Cumulant estimation The first two cumulants c 1 and c 2 are calculated using the Wavelet Leader and Bootstrap based MultiFractal analysis (WLBMF) toolbox (http://www.irit.fr/~Herwig.Wendt/software.html) which implements the empirical multifractal analysis formalism proposed by Wendt et al. [21]. Cumulants were estimated in the training and testing data for scales j between 3 and 7 (corresponding to a frequency range between 11.8 Hz and 187.5 Hz for which the power spectrum exhibits a 1/f-like behavior, see Fig 1) and moments of order q between -5 and 5 in 2s consecutive non overlapping sliding windows using 100 bootstrap resamples of wavelet leaders. A Daubechies wavelet with 3 vanishing moments (Daubechies 3) was used as a mother wavelet. Then the average cumulant was calculated in a 1-min window sliding at a 15s step. The length and overlap of the latter window were selected to match those used in our previous work [17]. This is to warrant an unbiased evaluation of the cumulants when comparing the performance of the seizure prediction method using cumulants to that using state similarity measures.

Statistical comparison of preictal and interictal cumulants In each patient independently, we compared the values of cumulant observations in 5-min preictal and 5-min interictal epochs of the training dataset. We hypothesized that the median of the cumulant observations does not statistically change when sampled from a 5-min preictal or from a 5-min interictal epoch. To verify this hypothesis, we compared average cumulant observations from the 5-min preictal and the 5-min interictal epochs using an unpaired two-tailed Wilcoxon rank-sum statistical test of the equality of medians. Preictal observations from available 5-min preictal training data were averaged. Similarly, interictal observations from consecutive 5-min epochs of training interictal data were averaged. The statistical comparison was carried in each channel independently and a Bonferroni correction was performed to control for multiple (channels) comparisons.

Cumulant versus spectral power A change in the state of vigilance (sleep/wakefulness) is generally accompanied by a significant change in the spectral power. This could affect the values of cumulant estimates as they indirectly measure inter-band frequency relations. To control for the state of vigilance as a possible confounder in the preictal and interictal cumulant comparison, we assessed whether any difference between preictal and interictal cumulant observations could only be due to a difference in the state of vigilance. Using the same sliding window we used to compute the average cumulant estimates (i.e. 1-min length, 15s step), we calculated the mean spectral power in the preictal and interictal epochs of the training dataset in the conventional EEG frequency bands: 0.5–3.9 Hz (delta), 4–7.9 Hz (theta), 8–14 Hz (alpha), 15–30 Hz (beta), 30–58 Hz (gamma), 80–250 Hz (ripples) and 250–450 Hz (fast ripples). The average 5-min preictal and interictal spectral power observations were then calculated the same way the average cumulant observations were obtained. By averaging the spectral power across multiple epochs, the confounding effect of the state of vigilance on the comparison between cumulant observations is minimized. The discriminability between preictal and interictal states using cumulants is tested by evaluating the difference between average preictal and interictal cumulant observations and comparing this difference to that observed between average spectral powers.

Seizure prediction using cumulant features Our original patient-specific seizure prediction method [17] uses a set of three measures, namely the persistence, distance and inclusion, referred to as state similarity measures, which quantify the similarity between the states underlying iEEG epochs and a reference state underlying the immediate preictal (90s) iEEG epoch. The state similarity measures are derived from the two dimensional thermodynamic profiles of the iEEG epochs obtained from time-series of wavelet energy and entropy calculated using the wavelet transform modulus maxima method (WTMM). We adapted our prediction method to use cumulants individually and in combination with state similarity measures as feature sets to predict seizures, with no change to the parameters used in training and testing procedures of the original algorithm. Fig 3 depicts a flow chart of the training and testing steps of the adapted version of the method. The outline of these steps is presented below. A detailed description of the concepts and procedures involved in each step can be found in the aforementioned reference. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 3. Flow chart of the seizure prediction method. Steps 1 to 3 are training procedures and 4 to 6 are testing procedures. Boxes in orange show procedures used instead when combination of features are used. https://doi.org/10.1371/journal.pone.0121182.g003

Training procedures All training procedures for the optimization of the classifier are performed on training data. Step 1- Feature computation: Cumulant estimates c 1 and c 2 (and state similarity measures) are calculated for each iEEG channel of the training dataset as previously described. Step 2- Discriminant analysis: A 10-fold leave-one-out cross validation is performed across preictal and interictal feature observations. A linear discriminant analysis is performed in each iteration. The classifier performance in each channel is iteratively assessed using a score combining the sensitivity and the specificity of the classification. Step 3- Performance assessment: Final channel scores are ranked and the top three channels with the highest scores are selected for testing procedures. When using state similarity measures, four frequency bands are independently assessed in each channel. In addition to channel selection, the frequency band yielding the best performance is also selected for testing procedures.

Testing procedures All testing procedures to assess the classifier are performed on test data in a quasi-prospective manner. Step 4- Feature computation: All features (cumulant estimates c 1 and c 2 and state similarity measures) are calculated for top ranked EEG channels (and selected frequency bands in the case of state similarity measures) in the sliding window. Step 5- Classification: Windows are labeled as ‘preictal’ or ‘interictal’ by the classifier in each channel separately. Channel decisions are then aggregated using majority voting rule for an ultimate classification of the window. Step 6-: Preictal change detection: Preictal changes are signaled through ‘preictal’ classification of five consecutive windows, corresponding to 2 min of iEEG continuously classified as ‘preictal’. Warnings are then raised for a variable period of time (referred to as period of active warning) controlled with the persistence-τ parameter. Raised warnings remain active for as long as a preictal change is detected, according to the persistence of warning lights protocol [32]. The parameter persistence-τ is equivalent to sum of the seizure occurrence period (the period during which the seizure is to be expected) and the seizure prediction horizon (the minimum time between the warning and the beginning of the seizure occurrence period) as defined in [33]. Warnings followed by iEEG interruptions lasting above three minutes are discarded since it is unknown if the patient was in ictal or interictal phase during this time and the classifier outcome could not be assessed.

Statistical validation The performance of the algorithm is compared to that of a chance predictor to assess the statistical significance of its prediction power. The improvement over chance is evaluated through a statistical comparison between the sensitivity of the algorithm and the sensitivity of the chance predictor which is demonstrated to be equal to the proportion of time spent in warning ρ [32]. More explicitly, we calculate the p-value of the significance of improvement over chance, given analytically in Equation (7), to evaluate the superiority of the prediction performance to chance at the 5% significance level. (7) where n is the number of seizures correctly predicted and N the total number of seizures.

Comparison of prediction performances Six combinations of cumulant and state similarity measures (feature sets FS1 to FS6) were used in the seizure prediction method: FS1: c 1 , FS2: c 2 , FS3: c 1 combined with c 2 ; FS4: State similarity measures combined with c 1 ; FS5: State similarity measures combined with c 2 and FS6: State similarity measures combined with c 1 and c 2 . Using the testing dataset, the prediction performance of the algorithm was evaluated for each feature set. The sensitivity, the proportion of time spent under warning, and the warning rate were assessed for ten persistence-τ values between five and 60 min. The false positive rate is also assessed for comparability with existing seizure prediction methods which reported the specificity in terms of false positive rate. The improvement over chance is assessed for each persistence-τ value separately. The sensitivity is defined as the proportion of test seizures that occurred within the period of active warning. The proportion of time spent under warning is the total time of active warning periods divided by the total time of test data. The warning rate is the number of all warnings divided by the total time of test data. The false positive rate is defined according to the recommendation of Mormann et al. [34]. It is the number of warnings not followed by a seizure occurrence within the period of active warning, divided by the total time of test data excluding the time of assumed preictal periods. The latter is defined as the sum of the period of time preceding each test seizure for which a false prediction cannot occur by definition. The time of assumed preictal periods can be obtained by multiplying the number of test seizures by the value of persistence-τ.