A Feature Extraction Method for Seizure Detection Based on Multi-Site Synchronous Changes and Edge Detection Algorithm

Automatic detection of epileptic seizures is important in epilepsy control and treatment, and specific feature extraction assists in accurate detection. We developed a feature extraction method for seizure detection based on multi-site synchronous changes and an edge detection algorithm. We investigated five chronic temporal lobe epilepsy rats with 8- and 12-channel detection sites in the hippocampus and limbic system. Multi-site synchronous changes were selected as a specific feature and implemented as a seizure detection method. For preprocessing, we used magnitude-squared coherence maps and Canny edge detection algorithm to find the frequency band with the most significant change in synchronization and the important channel pairs. In detection, we used the maximal cross-correlation coefficient as an indicator of synchronization and the correlation coefficient curves’ average value and standard deviation as two detection features. The method achieved high performance, with an average 96.60% detection rate, 2.63/h false alarm rate, and 1.25 s detection delay. The experimental results show that synchronization is an appropriate feature for seizure detection. The magnitude-squared coherence map can assist in selecting a specific frequency band and channel pairs to enhance the detection result. We found that individuals have a specific frequency band that reflects the most significant synchronization changes, and our method can individually adjust parameters and has good detection performance.


Introduction
Epilepsy seizures are brain disorders predominantly characterized by recurrent and unpredictable interruptions in normal brain function [1], which affect over 70 million people worldwide, with an incidence of 2.4 million per year [2][3][4]. Seizure detection is an important method to study the physiological mechanism, development process and treatment of epilepsy and has been widely used in clinical practice and scientific research. However, the most common detection method depends on the experience of clinical doctors and technicians, which places a burden on personnel and resources, and detection efficiency is subjective. Therefore, many automatic detection methods have been developed to help doctors save time and improve efficiency. For example, seizure signals can be automatically detected by finding the changes in the amplitude of the Electroencephalogram (EEG) signal. But there are still some problems with automatic seizure detection methods: the performance of detection methods in terms of accuracy, false alarm rate, and time delay needs to be improved; the computational complexity of some algorithms is high and difficult to be truly implemented; the detection features sometimes are not well adapted to different epilepsy seizure types. Therefore, further research is needed to improve the effectiveness of epilepsy detection.
In seizure detection algorithms, proper signal feature selection is the key factor in improving the efficiency of the detection result. Many features of single-channel seizure detection methods have been widely used, such as line length [5], amplitude [6], spectral features [7], time-frequency domain features [8][9][10], and features relating to information theory, such as approximate entropy [11]. Some features can achieve excellent performance but consume many computing resources, such as machine learning methods [12,13]. Some are preferable for online detection, with high computational efficiency [14], but the detection effect is not satisfactory, such as in linear detection methods. Synchronization is also a common feature in a multi-site method of seizure detection. The most used method for synchronization calculation is phase synchronization (PC) [15], which measures the interactions between rhythmic signals by detecting the instantaneous phase, and is unaffected by signal amplitude [16]. Depending on the position of the implanted electrode, the calculated synchronization results may increase [16][17][18] or decrease [19,20]. Further, synchronization changes between channels in different frequency bands have been found in recent studies [18,21,22]. In the PC method, it is difficult to observe synchronization changes in the full frequency band.
Another method to observe the synchronization change between signals in the frequency domain is called magnitude-squared coherence (MSC): a normalized cross-spectral density function that measures the strength of association and relative linearity between two stationary stochastic processes on a scale from zero to one [23]. In frequency-domain objective response detection, researchers compared MSC to PC using simulations with specific signal-to-noise ratios (SNRs) and a varying number of subaverages, and MSCs were superior to PC [24]. The MSC method is widely used in the fields of signal detection, time delay estimation, and SNR estimation [25][26][27][28], such as in human auditory models [29], electrocardiograms (ECGs), EEGs, magnetoencephalography (MEG) analysis [30][31][32] and seismic waves [33]. However, the current utilization of MSCs is limited in the frequency domain. MSCs only perform the calculation on a selected frequency, making it difficult to obtain a full view of the signals.
In this study, based on the existing MSC methods, we extended the MSC method to the time-frequency domain by drawing an MSC map to detect the synchronization change of two local field potential (LFP) signal channels. We used Canny edge detection algorithm to obtain information regarding the frequency band and when synchronization increases. We subsequently selected the frequency band where the synchronization change was the most apparent and calculated the maximal cross-correlation coefficient curve of the selected channel pairs after filtering. We used the average value (AVG) of correlation coefficient curves as the feature to detect seizure onset and added standard deviation (STD) as a parameter to reduce false alarms. This study focused on synchronization changes in the limbic system in a chronic spontaneous temporal lobe epilepsy (TLE) animal model. We found that individuals have a specific frequency band that reflects the most significant synchronization changes. Compared with other seizure detection features, the results showed that our features can individually adjust parameters and has good detection performance.

Animal Model and Data Acquisition
All surgical and experimental procedures were approved by the Animal Care and Use Committee of Zhejiang University. Pilocarpine-induced chronic TLE rat model was adopted in this study since this model closely resembles TLE in human. The procedure of inducing a chronic spontaneous TLE model is similar to our previous work [34]. Six adult Sprague-Dawley rats weighing about 250 g were injected with lithium chloride (12.7 mg/100 g i.p.) intraperitoneally. After 24 h, rats were treated with pilocarpine (3.5 mg/100 g i.p.) 30 min after being treated with atropine sulfate and supplemented by an extra half dose of pilocarpine after 30 min if status epilepsy (SE) was not observed. Once rats had SE status lasting over 90 min, a dose of diazepam (2 mg/100 g i.p.) was given to the rats to terminate the SE.
One week after drug induction, the rats in which spontaneous seizures were observed were selected for electrode implantation. All the rats were implanted with electrodes in the bilateral subiculum (SUB, AP: −6.0, ML: 3. The rats recovered for a week after surgery before being connected for long-term signal recording. One rat died two weeks after recording. The remaining five rats were included in the experiments. A homemade recording system was used according to a previous study [4]. This system contains a front-end recording (RHD2132 Intan Tech, LLC, Los Angeles, CA, USA) to transform neural electrical signals into digital signals. A field-programmable gate array (FPGA) module (XEM6010-LX45, Opal Kelly, Inc., Portland, OR, USA) was used for signal processing and system control. Each channel's LFP signal was recorded with a sample rate of 1 kHz and 16-bit resolution. The FPGA implemented a second-order IIR band stop filter of 50 Hz to filter out the power-supply noise. The data were then sent to an onboard secure digital card to be saved for subsequent analyses.

Construction of MSC Maps
As shown in Figure 1, the signals were calculated using the preprocessing and detection steps for seizure detection.
In the preprocessing step, the MSC map between every two channels was initially calculated to obtain the change in synchronization of different brain areas. In contrast to the phase coherence method, MSC is a power spectral density function that contains both amplitude and frequency information. Therefore, the MSC map of brain signals may reflect more seizure information [35]. The calculation formula for the MSC is as follows: where P xx ( f ) and P yy ( f ) are the power spectral densities of signals x and y, P xy ( f ) is the cross-spectral density of signals x and y. The C xy ( f ) indicates the coherence between x and y at each frequency. The result is between 0 and 1. Our study estimated the MSC function using Welch's overlapped average periodogram [36]. The Hamming window was used to divide the signal into segments with a 50% overlap on each segment. To construct a time-frequency scale MSC map to illustrate the change in coherence on the time scale, we split the signals by a 400 ms window (overlap 200 ms, minimum resolvable frequency f min = 2.5 Hz) and analyzed MSC results for every segment. Then, we spliced each result in the direction of the time series and obtained an MSC map with time on the horizontal axis and frequency on the vertical axis, as shown in Figure 1. The process of multi-site seizure detection method. The processing is divided into two steps: preprocessing and detection. In preprocessing step, 12-channel signal fragments are first calculated between every two channels to get MSC maps. The edge detection is then performed by the Canny's edge detection algorithm, and the result is shown as the red line in the figure. Select the original signals of the channel pair with detected edges for subsequent detection. Select the intersection of detection results as the frequency band. In detection step, the selected channels are filtered with the selected frequency band, and the seizure detection is completed through the steps of correlation coefficient calculation, feature extraction and threshold iteration.

Canny Edge Detection Algorithm
The Canny edge detection algorithm was applied to determine the MSC map's earliest and strongest coherence change area. Canny edge detection is a method of smoothing and derivative, and optimization approximation is obtained based on the measurement of SNR and location prodseuct [37]. The basic steps for the Canny edge detection algorithm are as follows [38,39]: (1) The original image is smoothed. The Gaussian filter is usually used, and in our study, the σ of the Gaussian filter was set to √ 2. The 2-D Gaussian filter function is as follows: (2) Calculate the magnitude and direction of the gradient using the finite-difference of the first-order partial derivative. The x direction first-order partial derivative: The y direction first-order partial derivative: The magnitude of the gradient: The direction of the gradient: (3) Non-maximum suppression: each pixel, M[x, y], is compared with its two neighboring pixels in the direction of the gradient θ[x, y]. If pixel M[x, y] is smaller by comparison, then M[x, y] = 0.
(4) A dual-threshold algorithm is used to detect and connect the edges. We set the high threshold value to 0.90 to remove most of the noise. The low threshold value is set to 0.4 times the high threshold value to add a connection between the edges.
As shown in Figure 1, the most dramatic correlation changes (red line) in the MSC map were detected, demonstrating the frequency range with the strongest correlation between the two channels. Other three different edge detection algorithms were also compared in this study, including Roberts, Sobel, and Prewitt. The comparison results are shown in Appendix A.
Under a uniform threshold, not all edges can be detected by the Canny edge detection algorithm in MSC maps, such as channels 11 and 12 in Figure 1. It indicated that some channel pairs do not increase visibly in correlation with the other. To screen out proper channels for detection, we selected a 5-s segment in each MSC map that contained the seizure onsets. We quantified the results by summing all pixel values detected by the Canny edge detection algorithm. We ranked the quantized results from largest to smallest and selected the first ten channel pairs as the frequency band and seizure detection channels. Physically, the edges were the frequency bands of the increment of the correlation. We intersected the edge of the selected channel pairs to find the common frequency band for later detection.

Calculation of Maximal Cross-Correlation Coefficient
The calculations of MSC and Canny edge detection algorithm were complicated, and detection cannot be performed promptly. Therefore, the maximal cross-correlation coefficient was further applied in the detection step as a synchronization indicator to reduce computing resources and improve detection efficiency.
The selected signals were passed with zero-phase digital filters, and the frequency range was decided upon in the preprocessing step. As in the construction of MSC map, we split the signals with a 400 ms window and an overlap of 200 ms. We subsequently calculated the maximal cross-correlation coefficient in each window and arranged the results in the time series to construct the cross-correlation curves. The maximal crosscorrelation coefficient was calculated for every channel pair as follows: where x(m) and y(m) are the discrete signal sequences, m is the number of each function value appearing in the sequence, and n is the time lags between x(m) and y(m). N is the length of signal sequence. R xy (n) can be normalized between 0 and 1 using the following formula: R xy, normalized (n) is the sequence associated with the time lags, where the maximal value is the maximal cross-correlation coefficient.
According to the Wiener-Khinchin theorem, the cross-correlation calculation results in the time domain can be transformed into the product form in the frequency domain: P xy ( f ) is the result of cross-spectral density of signal x and y. Therefore, the crosscorrelation results can represent the MSC results.
After calculating the maximal cross-correlation coefficient curves of selected channel pairs, we extracted the AVG of these maximal cross-correlation coefficient curves as a detection feature.The AVG value represents the overall enhancement of correlation, but is susceptible to common-mode noise. Therefore, it is necessary to add ten channel pairs with insignificant correlation enhancement, and calculate STD with the previously selected channel pairs. It is an important parameter to reduce false alarms.

Threshold Iteration
This study aimed to find a stable way to identify the starting time point of seizure activity. Two parameters, AVG and STD, were introduced to characterize the threshold of the seizure start time point. We used the initial five seizure activities of each rat to calculate the initial threshold using the following formula: where u is the average value of the AVG and STD curves, and sd is the standard deviation of the AVG and STD curves. We verified the reliability of the thresholds every three days, and the data from every third day was used for the iterative calculation of thresholds. If the initial threshold detected all seizure activities, and the false alarm rate was less than ten times per hour, the thresholds were reserved for later detection. If one or more seizure activities were not detected, or if the false alarm rate was more than ten times per hour, the threshold was recalculated based on the current AVG and STD data. The updated thresholds were subsequently used to detect seizures within the last data and ensure that all seizure activities could be detected. If the thresholds were valid, they were reserved for later detection (Figure 2).

Comparison
One of the multi-channel features, mean phase coherence (MPC) [15]was used to compare the effect with our method using the same data sets. MPC is a frequency-domain analysis method to quantify the degree of phase synchronization for two time series, x and y, defined as [20]: where 1/∆t is the sample rate of the discrete time series of length N. R is restricted to the interval [0, 1]. (When two signals are synchronized, R reaches the value of 1 and inversely reaches the value 0.) The φ a and φ b are defined as the instantaneous phases for signals x and y: wherex is the Hilbert transform of the signal x.
In addition, the maximal cross-correlation coefficient (CC) feature [40,41] without frequency band selection and channel selection is also added to the comparison. We also selected 7 common single channel features for comparison, including amplitude (AMP), line length (LL), variance (VAR), maximum slope (MSP), total power (TP), power spectral ratio (PSR) and approximate entropy (ApEn). The detection thresholds in the nine comparison features were treated the same way as ours. The calculation formulas refer to Appendix B.

Data Analysis
The data and statistical analyses were performed using MATLAB R2020b software (The Mathworks, Natick, MA, USA). Continuous variables with normal distribution are presented as mean ± standard deviation; non-normal variables are reported as medians. All recorded seizure activities were reviewed and labeled by experienced neurologists. The seizure detection efficiency is evaluated by detection rate, detection delay, and false alarm rate. Detection rate is calculated as the rate of automatic detection number and total seizure number. Detection delay is calculated as the difference between automatically detected time and manually marked time. False alarm rate is calculated as the number of false positives per hour.

Statistics of Seizure Activities
The number of intercepted seizure segments of the five rats was shown in Table 1. Experienced neurologists reviewed seizure activities from recorded brain signals combined with videos reviewed. A total of 206 seizure segments were identified.

Analysis of MSC Maps
Despite the same drug-induced seizure procedure in all rats, the results of MSC maps showed that each rat had its seizure activity mode. Two typical results are shown in Figure 3A,B. Canny edge detection showed a strong correlation change located in the high-frequency band of the MSC map of R-DG and R-CA3 channels from Rat-1 ( Figure 3A). With R-CA1 and L-DG channels from Rat-2, Canny edge detection result was located in the low-frequency band ( Figure 3B). We found the exact time point of the edge detection result, which in this case, appeared immediately before the seizure-like waves in the LFP signal. We further filtered the above LFP signals with the frequency band detected by the Canny edge detection algorithm to illustrate the signal correlation coefficient value between the two channels. The LFP signals from Rat-1 were filtered by frequency band of 180-500 Hz, and signals from Rat-2 were filtered by frequency band of 1-13 Hz. Both rats illustrated a marked change in CC values. By taking the detected edge as the dividing line, the CC value of Rat-1 changed from 0.57 before the edge was detected, to 0.95 after it was detected; Rat-2 changed from 0.45 to 0.85, respectively. For a single rat, not all channels showed an increased synchronization. We defined the quantified cross-correlation strength (QCCS) as the sum of all pixel values within 5 s before and after the edge detected in the MSC map: where t and f represent the two cumulative dimensions of time and frequency in the MSC map, the time dimension starts from point p start , 2.5 s before the edge, to 2.5 s after the edge. The frequency dimension ranges from 0 Hz to 500 Hz. The calculation result of QCCS is the cumulative sum of all pixel values V pixel in the region. Figure 4A illustrates the QCCS matrix of all MSC maps during one seizure onset of Rat-1, which had a strong correlation change in the high-frequency band. The coherence enhancement was mainly in the right and left hippocampal regions of Rat-1, indicating that seizure activity may originate from the hippocampal area and be enhanced by oscillating inside the two hippocampi. Figure 4B shows the result of the seizure onset of Rat-2, which had a strong correlation change at low frequencies. The coherence enhancement of Rat-2 mainly appeared between the right and left hippocampus and with the left ANT. Therefore, the seizure activity on Rat-2 may spread wider from one hippocampus to the other and other brain areas.
After the channel selection, we further determined the frequency band by intersecting the edge results of the selected channels. We calculated the frequency band every six days to ensure stability. The final frequency band selection is shown in Figure 4C, and the details are provided in Table 2. From the results of frequency band selection, rats 1, 4, and 5 had a strong correlation change in the high-frequency band and Rat-2 and Rat-3 in the low-frequency range.  For the selected channels, the stability is also evaluated every six days. The QCCS matrix is used for the evaluation, and the results are shown in Figure 5. As time changes, the synchronization between some channels changes. However, the intensity of synchronization and the distribution of related sites are relatively fixed. The channel pairs with high synchronization intensity of Rat-1, Rat-4 and Rat-5 are mainly distributed in the left and right hippocampus, and the channel pairs with high synchronization intensity of Rat-2 and Rat-3 are mainly distributed between the left and right hippocampus, as shown in the red box in Figure 5. The channels are selected from the red box to ensure the stability of long-time detection while ensuring the high discrimination ability in the detection process.

Detection Result and Validity
A typical example of the detection results from Rat-1 is shown in Figure 6A. The start of seizure activity was accompanied by a sharp increase in AVG and STD. Seizure onset was identified by both the AVG and STD thresholds. Generally, the thresholds of each rat were stable after five iterations. The threshold results for all five rats are shown in Table 3. The first five seizure activities calculated the initial thresholds. The thresholds of Rat-4 are stable at initialization, and the rest of the rats have minor threshold changes.  The frequency band selection was sensitive to seizure detection. As shown in Figure 6B, we compared the AVG and STD results in three different frequency bands.The curve clearly changed at seizure onset in the high-frequency band (280-420 Hz), which is selected by MSC maps and Canny edge detection algorithm. However, the curve changes slightly at seizure onset in the full-band and low-frequency band . The frequency band selection would affect the detection accuracy, false alarm rate, delay and other results. The channel selection was also sensitive to seizure detection. As shown in Figure 6C, we compared the AVG and STD results of the selected channels, all channels and unselected channels. Using the selected channels, the curve has more significant discrimination, which is more conducive to detection.

Comparison Results
We used this method to detect all seizures in this study, and the results are shown in Table 4. The detection rate of all rats is 96.60%. The average detection delay time after seizure onset is 1.25 ± 0.18 s. The false alarm rate of all rats is 2.63/h. Compared with the CC feature without frequency band and channel selection, our method can greatly improve the detection effect, and the effect is the best among all comparison features. In addition, multi-channel features such as CC and MPC detection performance are better than single-channel features. It should be noted that every feature has the possibility of optimization. This comparison only reflects this dataset's results, which may be limited and incomprehensive.
We also compared our method with other previous works, as is shown in Table 5. The comparison methods contain the features we use for comparison. It can be seen that the increase in the number of features used can lead to an increase in detection efficiency. In-sample parameter optimization can also improve detection results. Our method uses only one single feature, but by optimizing the selection of feature, we can also achieve the desired detection results. Schelter et al. [44] iEEG phase NO 20 70 (Sensitivity) synchronization Alaei et al. [45] scalp EEG mean-phase YES 88 100 (Sensitivity) coherence

Discussion
Our study developed a method to find a suitable feature for seizure detection in TLE model rats. We used the MSC map and Canny edge detection to verify the synchronization change phenomenon in the seizure rats and find the strongest change frequency band. We used the maximal cross-correlation coefficient to indicate the synchronization change in a computationally efficient and appropriate manner. We extracted two features of the correlation coefficient curves and iterated the thresholds to obtain the available feature thresholds. The feature extraction method for seizure detection proposed in this paper has the following advantages: high detection rate, low detection delay, low false alarm rate. In the preprocessing step, the feature extraction method has individually adaptive but has a certain computational complexity. Therefore, in the detection step, the feature extraction results can be transformed into filtering and cross correlation calculation to reduce the computational complexity.
Because it is difficult to build a Pilocarpine-induced chronic TLE rat model, we have limited data in this study. However, we can still get some interesting phenomena. One of them is the selection of the specific frequency band and channels. Although many studies have found synchronization changes in different frequency bands [21], some studies have investigated synchronization changes by constructing networks [46,47]. However, the reason for the synchronization changes still not fully understood. In this work, we found that the epileptic transmission structure may cause the frequency band difference; for the low-frequency band between the left and right hemispheres, the high-frequency band is shown in the hemispheres. This phenomenon needs to be confirmed by more experimental samples, and more experiments and methods are needed to explain the biological mechanism. This phenomenon can help differentiate epileptic types.
Furthermore, detecting lesion sites and key nodes is also an important task in epilepsy research, especially in clinical treatment, where doctors use CT, MRI, PET [48], and other medical imaging methods to locate the suspected lesions. Electrophysiological signals are the most direct and final means of confirmation for some patients with occult lesions. Localizing epileptic lesions is difficult and depends on the location of the implanted electrodes. The graph theory [49] and causal analysis [50] methods are widely used in network analysis and attempt to analyze epileptogenic areas or key sites with a limited number of electrodes. However, the computation becomes more complex with an increase in the number of electrodes. Our method may be helpful to preliminarily determine the correlation between channels and rank the importance of channel pairs. It provides a reference for researchers, and they can select nodes of interest for analysis according to their judgment and reduce computational complexity.
Based on the limitations of the current experiment, future research into the following is warranted: (1) Experimental verification of different species of organisms. The epilepsy characteristics of different species have similarities and differences. The synchronization of epilepsy in different individuals requires further investigation. Ultimately, the seizure detection method must to serve as an early warning and treatment for epilepsy. Our future studies will utilize the MSC map method to detect seizure activities in the epileptic brain of humans, confirming whether the method is effective.
(2) Identification of the seizure sites. Our study focused on finding the onset of seizure activity, and some important channel pairs can be selected. However, our study also had a limitation. The ten important channel pairs are changed during the progression of epilepsy, but we only used the initially selected channel pairs to detect the seizure to reduce computational complexity. The results of this detection are not superior. The location of lesions has also not been studied in detail. We can determine when the seizure begins by which sites have a strong correlation, which may be related to the location of the lesion, but further study is warranted.
(3) Dynamic analysis process. The MSC map can show the synchronization changes over the entire period, including the preictal, ictal, postictal, and interictal. It can help us to analyze the dynamic changes during seizures. Especially in recent years, the emergence of technicals such as sub-scalp EEG has made it possible to monitor patients' seizures for a very long time [51]. The method is suitable for long-term seizure recording and can help us determine what happens in the epileptic brain.

Conclusions
This study developed a feature extraction method for seizure detection based on multi-site synchronous changes in TLE. We investigated five chronic TLE rats with 8-and 12-channel detection sites in the hippocampus and limbic system. We used the MSC map and Canny edge detection algorithm to determine the synchronization changes in the specific frequency band and channel. We developed this phenomenon into a detection method and tested the seizure data of five rats. The method achieved high performance, with an average 96.60% detection rate, 2.63/h false alarm rate, and 1.25 s detection delay. The experimental results show that synchronization is an appropriate feature for seizure detection. The MSC map and edge detection algorithm can assist in selecting a specific frequency band to enhance the detection result. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available in the article.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Comparison of Detection Effects of Different Edge Detection Algorithms
In our study, we compared four different edge detection algorithms, including Canny, Roberts, Sobel, and Prewitt. The canny edge detection algorithm was finally selected to compare detection effects comprehensively. The matrix representation of the other three algorithms is shown in Figure A1A. The edge detection effect is shown in Figure A1B.

Appendix B. The Calculation Formula of Comparison Features
In our study, seven single-channel features were added for comparison. The calculation formula for these features is as follows: Amplitude: Line length: Variance: Maximum slope: Total power: Power spectral ratio: Approximate entropy: In the above formulas, x(i) represents a single channel signal,X represents the average of x(i). In the ApEn formula, the calculation formula of φ m is as follows: For the reconstructed m-dimensional vector X(1), X(2),... X(m), where X(i) = [x(i), x(i + 1), ...x(i + m − 1)], C m i (r) counts the number of vectors satisfying ||X(i), X(j)|| ≤ r (1 ≤ i ≤ N − m + 1). Generally, the value of parameter m is set as 2, and the value of parameter r is set as 0.2 times the standard deviation of data.