Real-time Inference and Detection of Disruptive EEG Networks for Epileptic Seizures

Recent studies in brain science and neurological medicine paid a particular attention to develop machine learning-based techniques for the detection and prediction of epileptic seizures with electroencephalogram (EEG). As a noninvasive monitoring method to record brain electrical activities, EEG has been widely used for capturing the underlying dynamics of disruptive neuronal responses across the brain in real-time to provide clinical guidance in support of epileptic seizure treatments in practice. In this study, we introduce a novel dynamic learning method that first infers a time-varying network constituted by multivariate EEG signals, which represents the overall dynamics of the brain network, and subsequently quantifies its topological property using graph theory. We demonstrate the efficacy of our learning method to detect relatively strong synchronization (characterized by the algebraic connectivity metric) caused by abnormal neuronal firing during a seizure onset. The computational results for a realistic scalp EEG database show a detection rate of 93.6% and a false positive rate of 0.16 per hour (FP/h); furthermore, our method observes potential pre-seizure phenomena in some cases.

models were then developed and trained with a large number of extracted features using state-of-the-art machine learning techniques (e.g., support vector machine, logistic regression, neural nets). The results generated by these methods promised a significant improvement of prediction or detection accuracy for various benchmark databases [17][18][19] .
However, many methods in the literature did not account for synchronous neuronal interactions across the brain 20 , and also were unable to characterize the dynamic nature of non-linearity and non-stationarity properly 21 . To this end, nonlinear dynamic learning approaches were proposed to present the epileptic brain as a dynamic network and investigate its spatiotemporal synchronization properties 12,22,23 . The underlying dynamics of global or local epilepsy can be assessed through the change of spatiotemporal synchronization (or dynamic network) patterns 24,25 . Some studies have shown great promises in seizure onset detection and prediction with network-based features, e.g., degree, shortest path, clustering coefficients, and the algebraic connectivity (Fiedler eigenvalue 26 ), characterizing the network property [27][28][29] .
Still, accurately capturing the time-varying dynamics of interconnected neurons (recorded by EEG signals) remains a challenging question due to the difficulty to infer networks from low signal-to-noise ratio data in a more effective way. Pioneering work on automatic recognition of epileptic seizures has demonstrated that it is possible to detect seizures by decomposing the EEG signals into elementary waves and analyzing rhythmic activities in the frequency band 3 to 20Hz 30 . Subsequently, a variety of graph-theoretic based approaches were proposed to construct the epileptic brain network using different techniques. For instance, some researchers computed the correlation 31,32 between EEG signals pair-wise and determined the existence of connectivity between channels by thresholding the coefficient values 32,33 , while others used a phase lag index which was computed using the angles obtained from the Hilbert transform of the EEG signals 33 .
In this study, we propose a new computational method that processes EEG signals in real-time, which in turn can be used for automated seizure control treatments in clinical practice. We utilize a technique for inferring the connectivity of networks (ICON) 34 , which reveals synchronization patterns during seizure onsets by estimating the coupling functions between EEG channels. The seizure onset is detected on a basis of the synchronization property of the recovered time-varying network inferred from multi-channel EEG data. To accelerate the seizure detection algorithm, we introduce a more computationally efficient approach derived from the ICON method and infer the time-varying network by computing the power of the coupling functions using Fourier transform (FT) applied on short epochs (0.25-1 second) of multi-channel EEG data (referred as FT method). Several methods in the literature have used the discrete wavelet transform (DWT) 35,36 in place of the Fourier transform that is considered to suffer from large noise sensitivity, but by combining the power and phase spectrum to infer the epileptic brain network, our approach significantly improves the robustness of the detection algorithm that effectively rejects artifacts. Furthermore, machine learning (ML) techniques are often used for seizure detection, and in the literature, various features are used as input to the ML model, e.g., the energy in EEG signals computed using wavelet transforms 17 . However, the performance of the ML models depends in part on the goodness of the feature(s) used. It is then expected that our graph-theoretic approach, which offers an excellent noise reject ratio, will improve the performance of seizure detectors.

Results
It is believed that epileptic seizures are characterized by widespread synchronous firing of neurons, and this synchronization often appears in a certain frequency band 33,37 , e.g., theta band (4-8 Hz) 38,39 . Hence, monitoring the brain electrical activity using EEG and analyzing its spectral information in the frequency domain can reveal seizure occurrences. However, EEG signals are often corrupted by noise, artifacts and/or brain activities (see Fig. 1a) that can be confused with seizure activities, and this can render accurate detection of seizure onsets very difficult. It is common practice to use either low-pass or band-pass filters to remove unwanted frequency components from EEG signals (see the right panel of Fig. 1a for a ten-second segment of non-seizure followed by seizure activities, respectively) before proceeding to its analysis [40][41][42] . The FT method presented herein, however, eliminates the need to use any such pre-processing of data. ft method. The performance of the proposed FT method was evaluated on the pediatric scalp EEG data from the CHB-MIT database available on the PhysioNet website 43 . The database contains 24 cases collected from 23 subjects. Note that in this paper we refer to each case by a patient number (from #1 to #24). We employed a graph-theoretic approach to fuse the information contained in the 23 EEG channels. This resulted in a better noise and artifacts rejection, even though short-time Fourier transform (STFT) was used (see Fig. 1). Indeed, one can see from Fig. 1b that there are significant peaks in the power spectrum in the time interval of 0 to 700 seconds, however, these power surges do not correspond to any seizure activity. Furthermore, from the second-smallest eigenvalue, λ 2 , of the Laplacian matrix of the inferred network, one can clearly distinguish the noise and artifacts from seizure activities (Fig. 1d). The Fiedler eigenvalue also defines the algebraic connectivity of a graph; in other words, it characterizes the graph synchronizability.
We tested our method on 17 patients, and excluded patients that had seizures characterized by amplitude depression (AD), and the results are summarized in Table 1. Overall, we achieved a sensitivity of 93.6% with a false detection rate of 0.16 FPs/h. The average seizure detection time was 16 seconds from the indicated seizure onset time in the dataset (this was obtained by considering a 6 second decision time). A closer look at the EEG signals, for those patients with high detection delay or latency, revealed that it was due, in most cases, to seizures starting with AD before the occurrence of synchronization that is characterized by coherent high amplitude oscillations. Since this algorithm is not designed to detect amplitude depression (because the detection is based on λ 2 going above a threshold), for these cases, we observed large latencies.
The FT method presented here achieves comparable performance for seizure onset detection with a lower false detection rate than most algorithms proposed in the literature 16,42,44 . However, our method presents some computational advantages. First, by using a short sliding window of size between 0.25 to 1 second, the algorithm requires less memory storage and has a very short runtime, hence making it a good candidate for implementation on portable devices for real-time seizure detection. The second advantage is that we achieved a very good performance with only one extracted feature (λ 2 ); however, extracting more features from the inferred graph, and using more sophisticated machine learning algorithms can significantly improve the sensitivity and the false detection rate. Note that the results presented here were obtained without using any machine learning algorithms, and instead a simple thresholding technique was used as shown in Fig. 1c. icon method. The performance of the ICON method was evaluated using the same dataset. Figure 2 shows two examples, case #16 for patient #1 and case #29 for patient #18, where ICON captured seizure onset through www.nature.com/scientificreports www.nature.com/scientificreports/ the abnormalities in synchronizability of the constructed dynamic brain networks. The synchronizability 45 , given by the exponential moving average (EMA) of λ 2 45,46 , can be derived directly from the dynamic interactions between brain regions inferred by ICON. Figure 2a illustrates a case (for patient #1, case #16) where ICON accurately captured seizure onset and offset times as recorded in the dataset, whereas Fig. 2b shows a case (for patient #18, case #29) of pre-seizure event that was detected 227 seconds before the actual seizure onset.
Moreover, Fig. 3 depicts the differences between the brain networks during normal brain activity (Fig. 3a) and seizure occurrence (Fig. 3b), inferred with ICON for the case presented in Fig. 2b. We can see that during a seizure, the strength of interactions between vertices, represented by weighted edges in Fig. 3b significantly increases. The summary of the results for all of the 24 patients is given in Table 2. Overall, a sensitivity of 78.79% was achieved, with a false detection rate of 0.02 FPs/h. However, we can achieve a sensitivity of 89.36% with a false detection rate of 0.01 FPs/h for only 17 patients considered for the FT method. The average seizure detection time was 4 seconds from the indicated seizure onset time and possible pre-seizure phenomena were observed for some patients (#11-#13, #15, and #18-#24).

Discussion
In this work, we provided a data-driven, dynmaic learning approach for seizure detection that was tested on the widely used benchmark CHB-MIT scalp EEG database 43 . The ICON method, which constructs a dynamic brain network to infer its connectivity from noisy measurement data, establishes a dynamic graph-theoretic approach to inferring the network dynamic topology and hence detecting the seizure onset and offset times by capturing the network's abnormalities. This work connects the brain network properties, such as connectivity and synchronizability 37 , with the occurrence of seizures. The hypothesis that there is a certain level of synchrony in brain signals during seizure led to a technique for inferring the epileptic brain network using mean phase coherence 47 .
This approach was used to analyze changes to graph centralities to shed more light on the role of constituents of evolving epileptic networks that recurrently transit into and out of seizures 48 . Furthermore, fundamental questions such as which nodes are connected by a predictive edge and which network modifications constitute a pre-seizure state were explored 48 . However, with ICON we mainly explore the seizure predictive capability of the second eigenvalue which is a measure of network algebraic connectivity. This method requires no prior information for data processing, and hence it can be applied to the entire dataset.
As seen in Table 2, the ICON method is sensitive to any network's abnormalities, which in turn guarantees the small latency for detecting seizure onset and offset times. Furthermore, since it can capture abnormality in the network dynamic topology, which also includes the one caused by neuronal disorders that might trigger the occurrence of seizures, this method can capture and reveal the patient-specific pre-seizure phenomena (See Fig. 2b). The challenge, however, is then how to relate these abnormalities to seizures while differentiating them from noise and artifacts, for example. This task will be difficult to achieve with the thresholding technique used in this paper, but using an appropriate machine learning model (e.g., a neural network) with multiple features extracted from the inferred network has the potential to uncover those hidden patterns that were difficult to detect, and can improve the detector accuracy or even predict seizure occurrence.
In general, the ICON method provides an innovative angle, through constructing a complex dynamic brain network, to detect epileptic seizures. This bridges the occurrence of seizures with its cause and effects, as shown in Fig. 3, which in turn enables future research studies on treating neuronal disorders like epileptic seizures.
While our ICON method can detect different types of seizures from the recovered network properties and hence, was tested on the entire dataset, the FT method cannot detect seizures that are characterized by AD. Therefore, patients #06, #12, #13, #14, #15, and #16 were excluded from the evaluation presented in Table 1. Although this method is not capable of detecting most seizure onsets for these patients, the false alarm rate is still very low. The graph-theoretic approach adopted here performed comparably well to some algorithms proposed in the literature (see Table 3). It is important to note that we achieved good performances without using sophisticated machine learning techniques and without additional features. There are other methods that are also suitable for real-time applications, e.g., the approach based on the phase-slope index of directed influence www.nature.com/scientificreports www.nature.com/scientificreports/ applied to multichannel electrocorticogram (ECoG) data 49 and which uses a threshold to assess seizures. In the paper by Rana et al., they used a more elaborate method for choosing the threshold based on the moving average of recent activity to accommodate variability between patients and slow changes within each patient over time 49 . Furthermore, by an adequate choice of parameters, they were able to achieve high sensitivity. In our work, on the other hand, we chose a constant threshold for each patient, and we are confident that with our approach, if we use a more elaborate technique for selecting a threshold or if we extract more features from the inferred synchronization network and combined with a good machine learning technique, we can improve the performance of the seizure detector. Recently, a paper has shown that using multiple features with a machine learning model (e.g., a neural network) improves the accuracy of the seizure detector 50 .
While the FT method can be used in automated real-time seizure control systems, the ICON method can help physicians better study pre-seizure phenomena. Indeed, this approach was able to identify brain activities that preceded seizure onsets as marked by specialists in the dataset. On average significant changes in the brain activity were observed 1.68 seconds before the actual seizure onset. The main advantage of the dynamic  www.nature.com/scientificreports www.nature.com/scientificreports/ graph-theoretic approach is its ability to reject noise and artifacts (see Fig. 1d), which allowed us to use a minimal number of features and yet, still achieve good performance of the detector. In particular, in Fig. 1d we show various features that have been obtained with other approaches such as the sum of spectral powers and the correlation method. The main observation is that our method generated a feature (λ 2 ) that is much smoother than the sum of spectral powers and λ n obtained with the correlation method. We also compared λ 2 to the eigenratio R = λ 2 /λ n , however, R was not as smooth as λ 2 as one can see in Fig. 1d, and for some other patients in the same dataset it was worse than shown in Fig. 1d, hence we do not recommend using it as a feature.

Methods
Data acquisition and processing. CHB-MIT scalp EEG dataset 43 is used for the validation of our proposed method, which consists of multi-variate EEG recordings from 23 epileptic patients at Boston Children's Hospital (5 males of age from 3 to 22; 17 females of age from 1.5 to 19). The data collection followed a protocol approved by the Committee on Clinical Investigations at the Beth Israel Deaconess Medical Center (BIDMC), Boston, Massachusetts, USA, and the Massachusetts Institute of Technology (MIT), Cambridge, Massachusetts, USA 51 . All procedures were performed in accordance with the relevant regulations and guidelines, and all informed consents were obtained before admission to BIDMC General Clinical Research Center took place. The EEG signals were recorded during many days using the International 10-20 system of 23 electrodes with a sampling rate of 256 Hz, and 182 epileptic epochs were marked by domain experts in a subset of 129 sessions of the entire dataset. Note that we only used the sessions with seizure onset epochs. After baseline removal and normalization, the EEG signals were re-referenced to average and a band-pass filter of 1-50 Hz was applied to reduce noise in the signals. These pre-processing steps were performed using functions in EEGLAB toolbox 52 . For real-time implementation, a non-overlapping sliding window of 1 second length was applied to each session of EEG recordings to split the signals into segments.
inferring connections of a brain network (icon method). In this paper, we provide a universal framework to diagnose seizure, by estimating a dynamic brain network and then revealing its connectivity from the noisy EEG data using our developed ICON (inferring connections of networks) method, which we have previously used to determine the topology of networks of oscillators, arising in electrochemistry, neuronal networks, and groups of mice 34 . Hence, by dynamically monitoring the connectivity of the inferred brain network, we capture abnormal activities such as seizures.
We study the evolution of a dynamic brain network consisting of N vertices. The brain activities are recorded using EEG sensors with N channels, and the dynamics of the brain regions near each channel is modeled by a state variable x i , where =  i N 1, , . The brain activities recorded at each vertex hence follow the dynamic law governed by the brain region's self-dynamics and the interaction among vertices, i n denotes the temporal state of brain region i at time t, the function f represents the inherent dynamics for the brain region i, and K ij , i, =  j N 1, , , is the coupling function between brain regions i and j (K ij can be different from K ji ). The dynamics, i.e., f and K ij , and thus the topology of such a directed brain network is problematic to infer from data because of the inherent nonlinearity within and between brain regions as well as the noisy measurement of the EEG data. However, the dynamic interactions between different regions will provide us insight to explain the brain activities through its first-order derivatives, which is more sensitive to capture the abnormality in this complex brain network and hence result in a lower false detection rate of identifying seizure onsets and/or offsets.
The central idea of our approach is to approximate the self and coupling dynamics of each region, f and K ij , respectively, by using complete orthogonal bases Q { } k k and P { } k k , in particular, the Fourier basis. Then, the dynamical law in Eq. (1) can be expressed as , that constitutes the Jacobian matrix J of the coupling term. Following this strategy based on the orthonormal basis representation and given the EEG data of each agent i in the network, the topology estimation can be formulated as a simple linear inverse problem for each agent i, given by where y (i) is the data vector and ∈ × + 2  is a matrix composed of the orthonormal bases, in which M is the number of data points for each sliding window of channel i, and r (r = 5) is the number of expansion terms in the truncated series; and z (i) is the coefficient vector that is being determined. Most importantly, this formulation enables independent estimation of the time-varying interactions between brain regions in such a network, so that the numerical computations become efficient (Fig. 4).
With the dynamic network topology revealed by ICON, we can study the synchronization behavior of a given brain network using eigenvalue analysis of the Jacobian matrix www.nature.com/scientificreports www.nature.com/scientificreports/ connectivity increases when a seizure occurs, the abnormality in the Fiedler eigenvalue, which indicates the synchronizability, would indicate the seizure onset. We achieved the seizure detection by using sliding windows of 10s-20s to guarantee the reliability of the topology recovered by ICON and efficiency of the detection. In this way, ICON can capture the abnormality of the brain network due to seizure efficiently and accurately.
inference of brain network from eeG power spectrum (ft method). For a significant number of patients suffering from epilepsy, seizures appear in EEG data as spiking trains, with frequencies often in specific bands, e.g., 4-8 Hz (theta band). However, for some patients' seizures can appear as amplitude depressions in the brain activity, or sometimes as amplitude death followed by spiking trains. In this paper, we exploit the burst of energy in EEG channels appearing in a given frequency band to form a time-varying synchronization network from which we then extract a single feature (the second eigenvalue, λ 2 of the graph Laplacian matrix). This feature is then compared to a fixed threshold to determine whether a seizure event occurred. The threshold selection was determined on a case by case basis because of the variability of EEG recordings from one patient to another. Furthermore, visual inspection was used to identify the seizure onsets. Various more elaborate methods have been proposed in previous studies, which could be adopted instead of using a simple threshold and visual inpsection, e.g., the method that determines the threshold based on the moving average of previous data points 49 and neural networks 50 .
Some of the automated seizure detection algorithms currently use Fourier 36 or wavelet transform 35 and look at the power spectrum in a certain frequency band 41 of some predetermined channels. These methods are usually patient dependent in the sense that one needs to know the morphology of the patient's seizure in addition to the channels that capture the seizure activity. However, it has been shown that algorithms that take into account different seizure morphologies, including amplitude depression, can have good performances and be less patient specific 53 .
In this work, we have presented the ICON method that infers the synchronization network from dynamical EEG measurements. Note that the coupling strengths are computed as α = ∑ ∑ ( ) coupling strength, we indeed see that the coupling strength between brain regions (vertices) i and j essentially represents the power of the interaction force between them (power of their interaction function). This analogy to the interaction signal power between vertices inspired a computationally efficient network inference approach, and hence enables online seizure detection.
In the simplified brain network inference approach, we employ the Fourier transform that is applied to the EEG data of each channel i. We choose a sliding window of the desired length, e.g., 1 second (or half a second) and compute the power and phase spectrum of each channel using the FFT (Fast Fourier Transform). After obtaining the power magnitudes p i (power spectrum) of each channel, we isolate a frequency band, typically between 2-10 Hz (may vary depending on patients), then set the magnitude of each frequency component outside the band [2,10] Hz, to zero. In order to extract only one feature, and hence suppressing the need to know in advance www.nature.com/scientificreports www.nature.com/scientificreports/ what channels are measuring seizure activities, we form a synchronization network as follows. Let i  be the sum of the squared power spectrum of the i th channel in a given time interval − − t t k k 1 (where k is the current epoch , , where N is the number of discrete frequencies in the desired band, e.g., [2,10] Hz. Hence, the ij th entry of the adjacency matrix A t ( ) k at time t k is obtained as Then the Laplacian matrix is computed as where D t ( ) k is the degree matrix. We can then extract a classification feature such as λ 2 at each time step t k (see Fig. 1d). Note that more features can be extracted from the network and used as inputs to a machine learning model, which we believe will improve the performance of the seizure detector. Note that for the FT method, the results were obtained with the EEG data normalized in the range [−1, 1], and the parameter γ = 2 was used. This parameter can be tuned to increase the robustness of λ 2 to noise and artifacts, however, a large value, i.e., γ > 3 might reduce seizure detection sensitivity because the value of λ 2 might be too small, hence causes some numerical issues.