An Explainable Statistical Method for Seizure Prediction Using Brain Functional Connectivity from EEG

Background Epilepsy is a group of chronic neurological disorders characterized by recurrent and abrupt seizures. The accurate prediction of seizures can reduce the burdens of this disorder. Now, existing studies use brain network features to classify patients' preictal or interictal states, enabling seizure prediction. However, most predicting methods are based on deep learning techniques, which have weak interpretability and high computational complexity. To address these issues, in this study, we proposed a novel two-stage statistical method that is interpretable and easy to compute. Methods We used two datasets to evaluate the performance of the proposed method, including the well-known public dataset CHB-MIT. In the first stage, we estimated the dynamic brain functional connectivity network for each epoch. Then, in the second stage, we used the derived network predictor for seizure prediction. Results We illustrated the results of our method in seizure prediction in two datasets separately. For the FH-PKU dataset, our approach achieved an AUC value of 0.963, a prediction sensitivity of 93.1%, and a false discovery rate of 7.7%. For the CHB-MIT dataset, our approach achieved an AUC value of 0.940, a prediction sensitivity of 93.0%, and a false discovery rate of 11.1%, outperforming existing state-of-the-art methods. Significance. This study proposed an explainable statistical method, which can estimate the brain network using the scalp EEG method and use the net-work predictor to predict epileptic seizures. Availability and Implementation. R Source code is available at https://github.com/HaoChen1994/Seizure-Prediction.


Introduction
Epilepsy is a group of chronic neurological disorders characterized by the abnormal and excessive fring of brain neurons, called epileptic seizures [1]. According to the the spontaneous electrical activity generated by brain neurons with high time resolution over a while. Scientists have found that it can be categorized into four diferent waveforms for scalp EEG records of patients with epilepsy. In the view of brain functional connectivity, these four diferent waveforms can be represented by four brain functional connectivity structures [7], corresponding to four diferent states of epilepsy seizures: (1) preictal state, which is the state before a seizure occurs; (2) ictal state, which is the onset state of seizure; (3) postictal state, which refers to the immediate state after a seizure; (4) interictal state, that is the state between postictal state and preictal state [8]. Figure 1 shows the sketch of the four states. Predicting seizures can be realized by detecting the preictal state, which can be achieved through discovering the changes in brain connectivity networks from interictal state to preictal state [9][10][11][12]. However, it is clinically difcult to identify the preictal state by visual inspection of scalp EEG signals to observe changes in the structure of brain connectivity networks. Terefore, powerful and explainable statistical methods are needed to determine the preictal state from scalp EEG recordings for seizure prediction.
Nowadays, brain functional connectivity modeling approaches have been proven to be a crucial tool in the neuroscience research feld [12]. Te brain can be seen as a complex network in which each brain region communicates and cooperates to carry out diferent functions [13]. However, the dysfunction in certain areas would interfere with the processing of upcoming information, consequently leading to network disorders and changes in a person's behavior [14]. Current research has shown that epilepsy is a specifc disease related to the brain network abnormalities, and they also suggest that the brain functional connectivity of a particular patient would abnormal dynamic changes during seizures, and the forms of brain functional connectivity are diferent among four diferent states of epilepsy seizures [15]. Hence, it is reasonable and adequate to employ brain functional connectivity modeling approaches to predict seizures. Some previous studies have applied brain functional connectivity modeling methods to study epilepsy disease. For example, Williamson et al. [16] constructed multiple spatiotemporal correlation structure features from EEG data to classify the patients' preictal or interictal states. A potential limitation of this study is that it only used a cortical network rather than a whole-brain network and cannot extract all essential features, making it difcult to achieve excellent predictive performance. Varotto et al. [17] proposed a method that employed a partially directed coherence method to depict the brain functional connectivity network. However, this study did not use the brain network features to predict seizures. Furthermore, a more recent work [18] has proposed an automatic seizure prediction method based on a graph convolutional network. Tis method could achieve a good seizure prediction performance by exploring the critical brain network features. To the best of our knowledge, this method is an excellent approach for seizure prediction. However, there are two fatal issues with this method. First, since this method is based on the deep learning technique, the interpretability is relatively weak. Second, the algorithm of this method is too complicated to be applied by clinicians.
To address the issues mentioned above, we referred to the existing novel statistical analysis framework called simultaneous diferential network analysis and classifcation for matrix-variate data (SDNCMV) [19] and based on the characteristics of scalp EEG data and epilepsy disease, we proposed an explainable statistical model for patient-specifc seizure prediction. Before introducing the proposed method, we frst briefy described the SDNCMV approach. Tis method was a two-stage data-driven approach that deals with fMRI data. Te frst stage estimated each subject's brain functional connectivity network and converted this network data into vector data for prediction in the next step. In the second stage, an ensemble prediction procedure was used to conduct the prediction results. In our study, we focused on using scalp EEG data to address patient-specifc prediction problems. Since EEG data and fMRI data are in the same data format as matrix data [20,21], we can refer to SDNCMV. However, there are still some diferences between the scenarios of this study and those of Chen et al. [19] that used fMRI data to classify Alzheimer's disease. Te fMRI data is a matrix data for each subject, while the EEG data is a matrix data for each epoch [22], which is artifcially generated. More specifcally, since the scalp EEG data used in this study have higher temporal resolution than fMRI data and the data epochs are small, we cannot use the SDNCMV method directly. We should assume that the brain functional connectivity for each subject is time-varying and then modify the frst stage of SDNCMV to estimate the dynamic brain functional connectivity to obtain better performance. More details of our method will be introduced in the following section.
Te signifcance of this study is that we proposed an explainable statistical method, which could be used to predict seizures based on the brain functional network. In addition, the proposed method is also computationally effcient, and the results can be easily interpreted. Terefore, this method can be better applied to the clinical and is more conducive to helping more patients with epilepsy.

Scalp EEG Data.
Te scalp EEG data used in this study are obtained from Children's Hospital Boston, the Massachusetts Institute of Technology (CHB-MIT) database [23] and Peking University First Hospital (FH-PKU) database.
Te CHB-MIT database is available with open access at https://physionet.org/physiobank/database/chbmit/. Te EEG recordings were collected from 23 children with intractable seizures, of which fve males with age from 3 to 22, 17 females with age from 1.5 to 19, and one child with missing gender and age data. Tese recordings were grouped into 24 cases since the EEG data of patient ID chb21 was obtained one-half year after chb01 from the same child. Te sampling frequency for this database was 256 Hz, and the international 10-20 EEG electrode positions and nomenclature system was used for these recordings. Te diference between two adjacent electrodes obtained the signal measurement for each electrode in this scalp EEG data. In most of the 24 cases, there are 21 unique signals, while a few cases contain less or more. Hence, only the recordings that include these 21 unique signals were selected in this study to keep the results consistent. Te FH-PKU database is a private database containing 17 patients or cases. Tis scalp EEG data were recorded at a sampling rate of 500 Hz and used 19 signals in the international 10-20 system. In this database, a diferent method was used from the CHB-MIT database to measure the signal of each electrode. Tis method used the signal diference between the electrode and the fxed reference electrode. Furthermore, to ensure the accuracy of the results, the physiological state of the patients in diferent epochs was roughly the same.
In this study, we focus on patient-specifc seizure prediction performance. For each patient in these two databases, the ictal state, the period when the patient experienced seizure onset, is easily detected from raw signals by doctors. Although the preictal state is challenging to identify and there is no gold standard, based on the ictal state, the preictal state can be defned by ourselves, which is the 30-minute window before a seizure occurs and is seen as the case state in this study. However, the interictal state, which is seen as the control state, is more difcult to defne; it is hard to recognize the period of the postictal state. Hence, to eliminate the noise efect of the postictal state, recordings within 2 hours after the end of the seizure are removed. Te period from this time to the next preictal state is defned as the interictal state. Moreover, if the time between two seizure periods is less than 2 hours, only the frst one is selected for this study. We consider all epochs as samples in this study and divide the continuous EEG data within a preictal state and interictal state into nonoverlapping 60-second epochs. To reduce computational complexity, we average the data obtained every second for each signal. Ten, the data of each epoch is in matrix form with 21 or 19 columns and 60 rows for the patients from CHB-MIT or FH-PKU database.

Dynamic Brain Functional Connectivity Estimation.
Tis section introduces the procedure to estimate the individual-specifc dynamic brain functional connectivity measures. For each individual, we denote X c ∈ R p×q and Y ϕ ∈ R p×q as the raw scalp EEG data matrix of c-th epoch for preictal state and ϕ-th epoch for the interictal state, respectively, where p represents the number of electrodes and q represents the number of time points. Based on the assumption that not every region of interest in our brain is connected, we estimate the dynamic brain functional connectivity within c-th epoch for the preictal state and ϕ-th epoch for the interictal state via sparse precision matrices for X c and Y ϕ , which estimate the strength measures of brain connectivity via partial correlations. Here we mainly focus on the procedure of how to address the raw scalp EEG data matrix in the preictal state X c , while Y ϕ can be dealt with similarly.
Before introducing the detailed procedure, we follow the classical matrix normal distribution framework to defne the distribution of X c . Assume X c follows a matrix normal distribution for each c, denoted as represent the covariance matrices of p electrodes locations and q time points for c-th epoch, respectively. Ten, for each time If the brain functional network is stable within each epoch, there are lots of existing approaches to estimate the sparse precision matrix Ω c X S � (Σ c X S ) − 1 in the high-dimensional setting, such as Graphical Lasso [24] and CLIME [25]. However, in the current study of epilepsy, it is more reasonable to assume that the brain functional network is changing over time. Hence, we need to estimate the dynamic sparse precision matrix Ω c X S (t) for each time point based on time-varying covariance matrices, which can be achieved by is a weighted covariance matrix, and we adopt a symmetric non-negative kernel function K(·) to generate the over time weights as w it � K(|i − t|/h n ). It is easy to fnd that this objective function is based on the Graphical Lasso and given the estimated Σ c X S (t), we can use the same algorithm to solve this optimization problem. Please refer to Graphical Lasso [24] for details. In practice, to obtain a better prediction performance, we substitute Σ c X S (t) by q t�1 Σ c X S (t)/q, which is the average of Σ c X S (t) over q time points and then a unifed sparse precision matrix estimation within each epoch, instead of multiple diferent sparse precision matrices Ω c X S , which can be used to estimate the brain functional connectivity measures. Although, here the sparse precision matrix is same for each time points within a specifc epoch, it is generated via a time varying covariance matrix, so the brain functional network can be considered dynamic. In addition, we use Gaussian kernel function and set h n � n 1/3 to calculate the weights w it in the real application, where n is the sample size.
To sum, in this study, we adopt two symmetric matrices W c to measure the dynamic brain functional connectivity strengths for c-th epoch in preictal state and ϕ-th epoch in interictal state, respectively. Here, we vectorize them via extracting the upper triangular elements by row for each matrix and connecting them together, and defne these vectors as V c for each epoch in diferent states, in which each element represent an edge in brain functional network and the dimension d is equal to p(p − 1)/2. Assume there are n 1 epochs in preictal state, n 2 epochs in interictal state and totally n � n 1 + n 2 epochs. Hence, we can defne the predictor matrix used in predictive model as V ∈ R n×p(p− 1)/2 , where V can be expressed as the stack of matrix V X S ∈ R n 1 ×p(p− 1)/2 and matrix V Y S ∈ R n 2 ×p(p− 1)/2 . In the following of this study, we are using V to serve as "Network Predictor Matrix." Given the network predictor matrix, due to the complexity and high dimensionality of the data, we adopt a penalized and ensembled logistic regression method to predict seizures. For details, we use Lasso penalty to deal with highdimension problem and bootstrap procedure to ensemble this penalized logistic regression. Let Z as the binary response variable and its observations are Z 1 , . . . , Z n , in which Z k � 1(k � 1, . . . , n 1 ) means corresponding observations are from preictal state and Z k � 0(k � n 1 + 1, . . . , n) means corresponding observations are from the interictal state. We denote P as the probability of Z � 1 and B(b � 1, . . . , B) is the bootstrap times. Now, we introduce this predictive method briefy, and for more details, please refer to Chen et al. [19]. We randomly sample n 1 epochs in the preictal state and sample n 2 epochs in the interictal state, respectively, with replacement. Ten, we repeat the resampling B times, and for each time, we employ the high-dimensional logistic regression model. We defne the β (b) as the corresponding regression coefcients vector estimated by b-th model. If there is the coeffcient in β (b) which is not equal to 0, it indicates that the corresponding edge in the brain network is meaningful for distinguishing between the preictal state and the interictal state. Finally, after whole resampling procedure, we defne the estimated coefcients vector as β i ≠ 0) to denote the weight for corresponding edge in the brain network. Te greater the weight, the more important this edge is. Tis demonstrates the interpretability of our method.

Performance Evaluation Measures.
Te evaluation measures that we adopt for the performance of seizure prediction are Sensitivity Rate (SENS), False Discovery Rate (FDR), and Area Under Curve (AUC). Te SENS measures the proportion of epochs from the preictal state with a positive result, and FDR is defned as the proportion of all epochs predicted from the preictal state, which is not. Since the values of these two criteria change with the cutof value, we select the cutof with the highest prediction accuracy. To avoid diferent cutofs afecting performance, we also present the AUC values, which is a cutof-independent measure and would be the most comprehensive measure.

Results
Tis section illustrates the results of the proposed method in seizure prediction by applying it to the CHB-MIT and FH-PKU databases. Table 1 presents the seizure prediction results for 24 patients in the CHB-MIT database. To prove that features derived via a dynamic brain functional network contain more information and improve prediction accuracy, we convert the raw scalp EEG data matrix for each epoch to a vector and stack these n vectors into a data matrix of dimension n × pq. Ten, we feed this data matrix into our ensemble prediction model for comparison. Furthermore, considering that the raw data may contain a small amount of information, we also combine the predictor matrix derived from brain network data (BN Data) and the predictor matrix derived from raw scalp EEG data (Raw Data) as a predictor to observe its prediction performance. Te results in Table 1 show that satisfactory prediction results can be obtained by using network features, while the prediction results obtained via using raw data as input are no diferent from random guessing. In addition, using the combination of brain network data and raw EEG data can also achieve satisfactory prediction. However, it is still worse than using brain network data only. Tis is because the raw data cannot provide any valuable information for prediction, so increasing these redundant variables from raw data makes prediction performance worse.
Te seizure prediction results for 17 patients from the FH-PKU database are presented in Table 2. As done in Table 1, we compare the prediction performance of our method using three kinds of input features. In Table 2, we show that network features can also achieve accurate predictions. Although the prediction result obtained using the raw data is higher than random guessing, it is still unsatisfactory. Furthermore, unlike the CHB-MIT data, in this dataset, since the raw data can provide some valuable information for prediction, it is found that the best prediction performance can be obtained using the combined data.

Discussion
In this study, we have proposed an explainable statistical method to predict epileptic seizures, which is helpful in raising the alarm before seizures. More concretely, our method uses scalp EEG data to construct a dynamic brain functional connectivity network via a time-varying precision matrix estimation approach. Ten, we treat these brain functional connectivity measures as predictor variables for the ensembled prediction model. Finally, through the proposed method in this study, we can obtain accurate prediction results by using these electrodes with overactive electrical discharges as predictors. In the following, we would like to discuss the fndings of this study, and at the end, we will present some future research directions.

Relationship to Other Studies.
Our study is not the frst to use scalp EEG data to extract brain connectivity signatures for seizure prediction. Tere have been lots of state-of-theart approaches, such as Gemein et al. [26], Tsiouris et al. [27], and Truong et al. [28]. So far as we know, the method called STS-HGCN-AL proposed by Yang et al. [18] can achieve a better seizure prediction performance among these approaches. We speculate that if our method outperforms STS-HGCN-AL, our method will outperform all existing methods. Hence, in this study, we only compare the prediction performance of our method with the method STS-HGCN-AL in the public database CHB-MIT database. It should be noted that it is difcult to draw a direct comparison due to diferent data preprocessing, such as how to choose the length of the epoch, when the preictal state starts, etc. In addition, the method STS-HGCN-AL has more strict requirements for the raw data, so this method can address part of patients' data within the CHB-MIT database, while our method can deal with all patients' data in terms of seizure prediction. We choose the AUC value as a measure of prediction performance to compare these two methods for the part of patients' data. From the results in Table 1, it is not difcult to see that our method can predict more patients and has a higher AUC value of 94% in these patient data, which can suggest that our method performs better than the existing method STS-HGCN-AL in terms of predicting epileptic seizures.
Furthermore, since our method uses brain functional connectivity features to predict seizures, in this study, in addition to the prediction performance, we also briefy discuss the performance of critical feature selection. We  Computational Intelligence and Neuroscience 5 randomly select two patients from the FH-PKU database, the patient 210486 and patient 210494. Based on the order of the derived weights ψ for each edge in the brain network, we show the top 10 critical brain functional connectivity features identifed by our model in Table 3 and Figures 2 and 3. Ten, based on the weights for each edge in Table 3, we may consider the features with larger weights as identifying potential connectome biomarkers. For the identifcation of potential connectome biomarkers, there are many existing methods, such as Song et al. [29], Lu et al. [30], and Ding et al. [31]. However, there is no gold standard for the public datasets identifying potential connectome biomarkers. We cannot prove the efectiveness of our method for this issue and thereby cannot compare our method with existing methods. We just put the identifcation results here without evaluating the performance, and for this issue, we will leave it as a direction for future research.

Limitations and Future
Directions. Although our method achieves satisfactory results in terms of seizure prediction, it still leaves much space for improvement to obtain even more realistic models. Te limitations of our method proposed in this study mainly concentrated on three aspects. At frst, for the sake of simplicity of our model, we have assumed that the raw scalp EEG data of each patient comes from a normal distribution. However, in the real world, we cannot know the actual distribution of the raw data. Hence, we are currently applying some statistical methods to relax the normal assumption for this problem.
Secondly, in our model, we have used the same time window to defne the preictal state of each patient, while in the real world, each patient has its heterogeneity and the time window of the preictal state is diverse. Hence, we plan to focus on how to estimate an optimal time window of the preictal state in the future. In the end, the study only briefy discussed how to identify potential connectome biomarkers using our method but did not test whether these biomarkers actually afected epilepsy. Terefore, in the future, we will conduct hypothesis tests in this feld to demonstrate the efciency and accuracy of our method in identifying potential connectome biomarkers.

Potential Applications of the Method in Treatment of
Epilepsy. Epileptic seizures are sudden and have no apparent signs. Te prediction of epileptic seizures can signifcantly enhance the efect of epilepsy treatment, improve the quality of life of patients with epilepsy, and reduce the mortality due to epileptic seizures, so the accurate prediction of epileptic seizures in the clinical application has a vital signifcance. Our experimental results and the comparison with previous work demonstrate that the proposed method is efcient. Tis gives the patient enough time to take action to cope with the seizure and reduce anxiety and trauma. Patients with epilepsy after regular ASMs treatment, there is still one-third of patients with epilepsy that cannot be controlled. Uncontrolled seizures have severe impacts on patients' cognition, memory, quality of life, social psychology, and the growth and development of children.
In recent years, imaging, electroencephalography, genetics, and other diagnostic techniques have been continuously improved, and the efcacy and safety of surgical resection have been recognized. For patients with drugresistant epilepsy with a clear epileptogenic focus and a low surgical risk, surgical resection should be considered as soon as possible. Te accuracy of the connectome biomarker identifcation may help to determine the epileptic region before epilepsy surgery, which is the key to ensuring the success of epilepsy surgery. Hence, fnding the right target remains the essential prerequisite for our new drug development, and accurate connectome biomarker identifcation can provide potential therapeutic targets. In this study, we have briefy discussed that our   Computational Intelligence and Neuroscience method may enable potential connectome biomarker identifcation, which is essential for physicians to conduct the preoperative evaluation and develop new drugs or treatments for epilepsy, but this needs to be further validated in the future.

Consent
Te parents of the patients signed written informed consent and agreed with their children's participation in this study and allowing the use of the relevant data and information for scientifc research.

Data Availability
Te original contributions presented in this study are included in the article; further inquiries can be directed to the corresponding authors.

Ethical Approval
Ethical approval for this study was obtained from the Ethics Committee of the Peking University First Hospital.

Disclosure
Hao Chen and Taoyun Ji are the co-frst authors.

Conflicts of Interest
Te authors declare that the research was conducted without any commercial or fnancial relationships that could be construed as a potential conficts of interest.