Using Constrained Square-Root Cubature Kalman Filter for Quantifying the Severity of Epileptic Activities in Mice

(1) Background: Quantification of severity of epileptic activities, especially during electrical stimulation, is an unmet need for seizure control and evaluation of therapeutic efficacy. In this study, a parameter ratio derived from constrained square-root cubature Kalman filter (CSCKF) was formulated to quantify the excitability of local neural network and compared with three commonly used indicators, namely, band power, Teager energy operator, and sample entropy, to objectively determine their effectiveness in quantifying the severity of epileptiform discharges in mice. (2) Methods: A set of one normal and four types of epileptic EEGs was generated by a mathematical model. EEG data of epileptiform discharges during two types of electrical stimulation were recorded in 20 mice. Then, EEG segments of 5 s in length before, during and after the real and sham stimulation were collected. Both simulated and experimental data were used to compare the consistency and differences among the performance indicators. (3) Results: For the experimental data, the results of the four indicators were inconsistent during both types of electrical stimulation, although there was a trend that seizure severity changed with the indicators. For the simulated data, when the simulated EEG segments were used, the results of all four indicators were similar; however, this trend did not match the trend of excitability of the model network. In the model output which retained the DC component, except for the CSCKF parameter ratio, the results of the other three indicators were almost identical to those using the simulated EEG. For CSCKF, the parameter ratio faithfully reflected the excitability of the neural network. (4) Conclusion: For common EEG, CSCKF did not outperform other commonly used performance indicators. However, for EEG with a preserved DC component, CSCKF had the potential to quantify the excitability of the neural network and the associated severity of epileptiform discharges.


Introduction
Epilepsy is a common disorder of the central nervous system caused by abnormal or excessively synchronized firing of neurons. Sixty percent of adults with epilepsy are classified as having focal epilepsy, with temporal lobe epilepsy accounting for most cases [1,2]. Since the foci of temporal lobe epilepsy are commonly found in the hippocampus, many studies on epilepsy have focused on this region [3]. Oral antiepileptic drugs are currently the most common epilepsy treatment, however about 30% of patients cannot be effectively treated with drugs. Many alternative treatments have been investigated, including electrical stimulation [4][5][6][7]. Many parameters and waveforms can be adjusted when designing electrical stimulation protocols. The frequency range was roughly divided into low Figure 1 depicts a physiologically based neural mass model of electrical activities of hippocampus [34]. The model includes a pyramidal cell population, fast and slow inhibitory interneuron populations, and an excitatory interneuron population. The relationship between the presynaptic firing rate and postsynaptic membrane potential of each neuron population was modeled using a second-order linear transfer function. The postsynaptic membrane potential was then transformed into the firing rate of the postsynaptic neuron population using a sigmoid function (S), where e s , h s , and v s are constants, and their values are listed in Table A1 of the Appendix B.
The corresponding state-space equations are given below. .
where C i , i = 1 to 7, and C S , a, b, g, and G are constants, and their values are listed in Table A1 of Appendix B. A and B are parameters that are adjusted to generate normal and epileptic EEGs. The right-most transfer function, s/(s + 0.5), is a high pass filter which emulates the signal processing of an EEG recorder. The output, y, is the simulated EEG.
Where es, hs, and vs are constants, and their values are listed in Table A1 of the Appendix B. The corresponding state-space equations are given below.
where Ci, i = 1 to 7, and CS, a, b, g, and G are constants, and their values are listed in Table  A1 of Appendix B. A and B are parameters that are adjusted to generate normal and epileptic EEGs. The right-most transfer function, s/(s + 0.5), is a high pass filter which emulates the signal processing of an EEG recorder. The output, y, is the simulated EEG.

Formulation of CSCKF and Bifurcation Analysis
In this section, CSC−F was used to estimate the states of the mathematical seizure model, as shown in Figure 1, as well as the time-varying parameter B(t). First, the state space model (Equation (2)) was expanded, 4 of 17 where w c (t) and ν c (t) are the state and output noise, respectively, in which Gaussian distribution was assumed, with zero mean and covariance matrix Q and R, respectively. x c is the augmented state variable array, where x 1 to x 11 are the original state variables and x 12 corresponds to B(t). In other words, this approach assumes that B(t) is a state variable. On the other hand, A(t) and G(t) are time-invariant constants, i.e., A(t) = 5 and G(t) = 10. In this way, f c (x c ,t) and h c (x c ,t) were redefined as Because the experimental data were all digitalized signals, the next step was to digitize the model, i.e., .
x c (t(k)) where the sampling rate is 1/T, k is the indexed time point, and x c (t(k)) and y c (t(k)) are abbreviated as x k and y k , respectively.
The Kalman filter was calculated using a recursive algorithm of two alternating steps.
Step 1 was the time update, i.e., updating the state estimates to the present time (k) based on the data of previous time (k − 1). S k−1|k−1 was obtained from step 2 of the previous time point (k − 1). The states at the cubature points (i = 1, 2 . . . m; m = 2n x ) of the previous and present times were estimated, where [ζ i ] represents the cubature points. Using the definition of cubature points, the estimates of present states based on the data up to the previous time point were obtained, and the square-root factor (S k|k−1 ) of the predicted error covariance was updated to the present time, where Tria is a general triangularization function to produce a lower triangular matrix, and where Tria is a general triangularization function to produce a lower triangular matrix, and ▌ (End of step 1) Step 2 was the measurement update, i.e., updating the data to be used when estimating the present states to the present time. First, Sk|k−1 was obtained from the results of step 1, and the states (Xi,k|k−1) and system output (yi,k|k−1) at the cubature points (i = 1, 2 … m; m = 2nx) of the present time were estimated, Then, using the definition of cubature points, the estimate of system output at the present time point was obtained, (13) and the square-root factor of the predicted error covariance was updated to the present time, where Next, cross-covariance (Pxy,k|k−1) was calculated, (End of step 1) Step 2 was the measurement update, i.e., updating the data to be used when estimating the present states to the present time. First, S k|k−1 was obtained from the results of step 1, and the states (X i,k|k−1 ) and system output (y i,k|k−1 ) at the cubature points (i = 1, 2 . . . m; m = 2n x ) of the present time were estimated, Then, using the definition of cubature points, the estimate of system output at the present time point was obtained,ŷ (13) and the square-root factor of the predicted error covariance was updated to the present time, where Next, cross-covariance (P xy,k|k−1 ) was calculated, and Kalman gain (W k ) and states were updated to the current time point (k) based on the data of the present time, Finally, the square-root factor (S k|k ) of the predicted error covariance was updated based on the data of the present time, where [ζi] represents the cubature points. Using the definition of cubature points, the estimates of present states based on the data up to the previous time point were obtained, (9) and the square-root factor (Sk|k−1) of the predicted error covariance was updated to the present time, where Tria is a general triangularization function to produce a lower triangular matrix, and ▌ (End of step 1) Step 2 was the measurement update, i.e., updating the data to be used when estimating the present states to the present time. First, Sk|k−1 was obtained from the results of step 1, and the states (Xi,k|k−1) and system output (yi,k|k−1) at the cubature points (i = 1, 2 … m; m = 2nx) of the present time were estimated, Then, using the definition of cubature points, the estimate of system output at the present time point was obtained, The ratio I B/A = x 12 /A, representing the estimated B(t)/A(t), was used as the indicator of excitability of the model in the later analyses. Potential types of dynamic behavior and bifurcation in the epilepsy model were analyzed using MatCont software [35].

Animal Preparation and Experiment Setup
The C57BL/6 mice were housed under standard conditions in the animal facility of the College of Medicine, National Cheng Kung University (NCKU) and kept on a 12 h light/dark cycle with food and water ad libitum. The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Institutional Review Board of National Cheng Kung University (IACUC Approval No: 108223).
A set of three electrodes with two electrodes in spiral form were implanted into the CA3 region of the hippocampus (AP − 1.34 mm, ML ± 1.25 mm and DV − 1.8 mm from the bregma) using a stereotaxic technique under isoflurane (Attane, Step, Taipei, Taiwan) anesthesia. Two stainless steel screws were fixed on the skull above the cerebellum and olfactory bulb as the ground and reference electrodes, respectively. The spiral electrodes were used for electrical stimulation and the other implanted electrode and the screw electrodes were used for EEG recording. The experiments were performed 2 weeks after electrode insertion when the animals had fully recovered from surgery. Acute seizures were induced by intraperitoneal (ip) injection of a bolus of kainic acid (ab120100, abcam, Boston, MA, USA) (10 mg/kg), followed by boluses of 5 mg/kg per 0.5 h until the behavior of mice was rated as 5 on the Racine scale. The injection rate was then reduced to boluses of 2.5 mg/kg/h. After the experiment, diazepam (Dupin, China Chemical & Pharmaceutical Co., Taipei, Taiwan) (ip: 10 mg/kg) was given to stop the seizures.
EEG signals were first amplified 1000× and hardware filtered (0.1 Hz-1000 Hz bandpass and 60 Hz notch filters) by a biopotential amplifier (Differential AC Amplifier model 1700, A-M System, Sequim, WA, USA), and then sampled at 3840 Hz into a personal computer using a control system development tool (DS1104, dSPACE, Paderborn, Germany). The behavior of the mice was videotaped throughout the experiment to assist later analysis.
The DS1104 also output control signals to drive an electrical stimulator (A395 Linear Stimulus Isolator, World Precision Instruments, Sarasota, FL, USA) to produce electrical stimulation in constant current mode. For the RNS, the control signal was band-limited (101-640 Hz) white noise, and the amplitude was adjusted so that the extreme amplitude of actual current output was ±500 µA. For the RPS, the control signal was a sequence (5 s) of 128 Hz biphasic and charge balanced square-wave pulses, with the positive phase 458.5 µA in height and 781.25 µs in duration, and the negative phase −50.94 µA in height and 7031.25 µs in duration. Sham stimulation had identical waveform parameters as the real stimulation, except the height was only 1% that of the real stimulation.
The electrical stimulation was given manually when type III epileptic waveforms, i.e., 4-12 Hz spikes lasting longer than 5 s, were detected by visual inspection. The rationale for choosing this type of seizure was (1) that the animal had a Racine score of 3 and showed no obvious convulsive movements, thus preventing injury to the animal and electrode setup, and (2) that it was relatively objective to judge whether epileptic activities were suppressed.
In total, 20 C57BL/6 mice (age: 46.2 ± 5.6 weeks and body weight: 33.3 ± 4.5 g) were used in this study, of which 6 received RNS, 9 received RPS, and the remaining 5 were the controls and did not receive electrical stimulation. The stimulation shifted between real and sham stimulation pseudo-randomly. Consecutive stimulations were separated by at least 1 min, and repeated until no more consistent type III epileptic waveforms were noted.
EEG was recorded throughout the whole experiment. A 209-order equiripple finite impulse response low-pass filter with a pass band below 50 Hz, a stop band above 100 Hz, pass band ripple less than 0.1 dB, and stop band attenuation of 60 dB was used to remove stimulation artifacts.

Formulation of Other Quantifying Indicators of Epileptic Activities
Another three indicators were formulated for comparison. The first was band power (P b ) in the frequency range of 4-50 Hz. The second was Teager energy operator (E T ) [36].
The third was sample entropy (SEn) [26,37], a variant of approximate entropy with two advantages: data length independence and relatively trouble-free implementation. Sample entropy was calculated as follows [26]. Initially, a set of vectors, X m (i), was formed from an EEG segment of length N, , was defined as the Chebyshev distance and the number of vector pairs of d[X m (i), X m (j)] < r·SD and d[X m+1 (i), X m+1 (j)] < r·SD, where r is a constant and SD is the standard deviation of the EEG segment about the mean, denoted as B m (r) and A m+1 (r), respectively. SEn was then defined as In this study, m, r, and N were set as 4, 0.1, and 2400, respectively. The time-window width of the data segment to calculate all three indicators was 5 s.
All the calculations and simulations were carried out using the commercial software package, Matlab and Simulink (MathWorks, Natick, MA, USA).

Statistical Analyses
Group results were summarized as mean ± one standard deviation. Two-way mixed design ANOVA with treatment, i.e., electrical stimulation and time (before, during and after stimulation) as the main factors, and post-hoc analysis by non-parametric Wilcoxon signed-rank test were used to test the efficacy of electrical stimulation. The significance level was set at p < 0.05. Multiple comparisons using non-parametric Mann-Whitney U tests with Bonferroni correction were used to test the difference of performance indicator between types of EEG. The significance level was set at α < 0.01.

Performance of Quantitative Indicators in Simulated EEG
The mathematical model described in the Methods section was used to generate three sets of one normal and four types of epileptic activities as shown in Figure A1 in Appendix A. Each combination of A and B was used to generate 50 EEG segments of 5 s. The excitability of the neural network model decreased from Type I to Type IV, i.e., the severity of epileptic discharges decreased from Type I to Type IV.
The results of the four indicators are shown in Figure 2. The results of I B/A largely deviated from those expected, i.e., I B/A would increase from Type I (around 3) to Normal (around 10). For Type I epileptiform discharge, some estimated I B/A values were around 3 as expected, but many values were at the other extreme (around 15). The estimates of Types II, III, and IV, in the range of 1 to 4, were lower than those of Type I. The estimates of normal EEGs matched the expected values. P b showed a trend compatible with the subjective judgement, i.e., the power increased from Type I to Type II, and then decreased from Type II to Type IV. In addition, the corresponding P b values also increased with A, the parameter of the model that made a major contribution to the amplitude of the spikes. The wave amplitudes of all types increased with A but at different rates, so that the relative amplitudes of Type I and Type IV could reverse due to different A values. The results of E T and SEn were similar to those of P b . In short, all four indicators showed a trend compatible with subjective human judgement. However, the trend did not match the trend of excitability of the neural network, i.e., B/A.
Of note, the trends of the four indicators were quite similar. Except for Type I EEG, the severity of seizures decreased from Type 2 EEGs to normal EEGs. In addition, the variance was the largest for Type I EEGs in all four performance indicators, and especially for I B/A .

Performance of Quantitative Indicators in Simulated Output of the Seizure Model
The right-most transfer function after the model output, z(t) ( Figure A2), was used to emulate the high-pass filtering properties of commonly used EEG machines (Figure 1). The split phenomenon in estimating parameters of Type 1 EEGs by CSCKF was suspected to be due to the filter. Therefore, the performance of the four quantitative indicators was re-investigated using z(t) as the output from the model (Figure 3). Except for I B/A and P b , the results of the other three indicators were almost identical to those using y(t) as the output. For P b , all values increased, but the relative trend was maintained. For I B/A , the split phenomenon disappeared, and all estimates were concentrated around the expected value of 3 (Figure 4 right). In addition, the trend of estimates was consistent with the setting, i.e., graded increase of I B/A from Type 1 EEGs to normal EEGs. The variance of estimates for each type of EEG also decreased.
as expected, but many values were at the other extreme (around 15). The estimates of Types II, III, and IV, in the range of 1 to 4, were lower than those of Type I. The estimates of normal EEGs matched the expected values.
Pb showed a trend compatible with the subjective judgement, i.e., the power increased from Type I to Type II, and then decreased from Type II to Type IV. In addition, the corresponding Pb values also increased with A, the parameter of the model that made a major contribution to the amplitude of the spikes. The wave amplitudes of all types increased with A but at different rates, so that the relative amplitudes of Type I and Type IV could reverse due to different A values. The results of ET and SEn were similar to those of Pb. In short, all four indicators showed a trend compatible with subjective human judgement. However, the trend did not match the trend of excitability of the neural network, i.e., B/A.

Performance of I A/B in the Experimental Data without Electrical Stimulation
The estimated I B/A values of experimental EEG data without electrical stimulation are shown in Figure 5. It is clear that most of the normal EEGs were estimated to have a larger I B/A value of around 10, while most of the epileptiform EEGs were estimated to have a lower I B/A value of around 3. However, some epileptiform EEGs had higher I B/A values than those of normal EEGs. In addition, the normal and epileptiform EEGs were segregated into two groups, and there were only a small number of intermediate points between these two peaks, indicating a trend of dichotomy. In other words, I B/A could distinguish between normal and epileptiform discharges, but could not quantify their severity. This phenomenon, i.e., split of estimates to two ends in Type 1 EEGs, was also similar to those seen in simulated EEGs (Figures 2 and 4).
the results of the other three indicators were almost identical to those using y(t) as the output. For Pb, all values increased, but the relative trend was maintained. For IB/A, the split phenomenon disappeared, and all estimates were concentrated around the expected value of 3 (Figure 4 right). In addition, the trend of estimates was consistent with the setting, i.e., graded increase of IB/A from Type 1 EEGs to normal EEGs. The variance of estimates for each type of EEG also decreased.    Figures 2 and 3, respectively. In the case of z(t), the trend of estimates was consistent with the setting, i.e., graded increase of IB/A from Type 1 EEGs to normal EEGs. In addition, three bands, corresponding to three values of A ( Figure A1), can be seen. The boxplot for each group has a standard form, displaying the inter-quartile range (box), the range of inliers (error bar), and the median (orange line).

Performance of IA/B in the Experimental Data without Electrical Stimulation
The estimated IB/A values of experimental EEG data without electrical stimulation are shown in Figure 5. It is clear that most of the normal EEGs were estimated to have a larger  Figures 2 and 3, respectively. In the case of z(t), the trend of estimates was consistent with the setting, i.e., graded increase of I B/A from Type 1 EEGs to normal EEGs. In addition, three bands, corresponding to three values of A ( Figure A1), can be seen. The boxplot for each group has a standard form, displaying the inter-quartile range (box), the range of inliers (error bar), and the median (orange line). gated into two groups, and there were only a small number of intermediate points tween these two peaks, indicating a trend of dichotomy. In other words, IB/A could dis guish between normal and epileptiform discharges, but could not quantify their sever This phenomenon, i.e., split of estimates to two ends in Type 1 EEGs, was also simila those seen in simulated EEGs (Figures 2 and 4).

Performance of Four Quantitative Indicators in RNS Experiments
Thirty triplets, before, during, and after stimulation, of EEG segments (5 s) from mice receiving RNS were used in the following analysis ( Figure 6). The group resu (mean ± standard deviation) of IB/A before, during, and after real stimulation were 5.8 4.46, 5.80 ± 4.63, and 2.04 ± 0.92, respectively. Two-way mixed design ANOVA show that both time (F2,96 = 10.80, p < 0.001) and treatment (F1,48 = 13.84, p < 0.001) were signific factors, and that there was a significant interaction between the effects of the two fact (F2,96 = 4.34, p = 0.016).

Performance of Four Quantitative Indicators in RNS Experiments
Thirty triplets, before, during, and after stimulation, of EEG segments (5 s) from six mice receiving RNS were used in the following analysis ( Figure 6). The group results (mean ± standard deviation) of I B/A before, during, and after real stimulation were 5.80 ± 4.46, 5.80 ± 4.63, and 2.04 ± 0.92, respectively. Two-way mixed design ANOVA showed that both time (F 2,96 = 10.80, p < 0.001) and treatment (F 1,48 = 13.84, p < 0.001) were significant factors, and that there was a significant interaction between the effects of the two factors (F 2,96 = 4.34, p = 0.016). The group results (mean ± standard deviation) of Pb before, during, and after real stimulation were −3.02 ± 5.09, −4.78 ± 4.79, and −2.52 ± 4.39, respectively. Two-way mixed design ANOVA showed that time was a significant factor (F2,116 = 5.488, p = 0.005), but treatment was not (F1,58 = 0.043, p = 0.837), and there was a significant interaction between the effects of the two factors (F2,116 = 4.523, p = 0.013). The group results of ET before, during, and after real stimulation were −2.99 ± 0.53, −2.76 ± 0.21, and −2.61 ± 0.24, respectively. Two-way mixed design ANOVA revealed that time was a significant factor (F2,106 = 10.581, p < 0.0001), but treatment was not (F1,53 = 2.83, p = 0.098), and there was a significant interaction between the effects of the two factors (F2,106 = 9.482, p < 0.001). The group results of SEn before, during, and after real stimulation were 0.33 ± 0.12, 0.62 ± 0.21, and 0.23 ± 0.12, respectively. Two-way mixed design ANOVA showed that both time (F2,112 = 49.51, p < 0.001) and treatment (F1,56 = 6.12, p = 0.016) were significant factors, and that there was a significant interaction between the effects of the two factors (F2,112 = 56.55, p < 0.001). Figure 6. Results of the four indicators on EEG segments collected from mice before, during, and after random noise stimulation. IB/A: estimated B/A from CSCKF, Pb: band power in the frequency range of 4-50 Hz, ET: Teager energy operator, and SEn: sample entropy. The colored bar shows the mean value of each group, and the error bar shows the corresponding range of ±1 standard deviation. The number within the parentheses is the number of EEG segments, where the outliers have been excluded by the 1.5 inter-quartile range rule. **: Significant difference (p-value < 0.01) was tested by using non-parametric Wilcoxon signed-rank tests.

Performance of the Four Quantitative Indicators in RPS Experiments
Twenty-seven triplets before, during, and after stimulation of EEG segments (5 s) from nine mice receiving RPS were used in the following analysis (Figure 7). The group results (mean ± standard deviation) of IB/A before, during, and after real stimulation were Figure 6. Results of the four indicators on EEG segments collected from mice before, during, and after random noise stimulation. I B/A : estimated B/A from CSCKF, P b : band power in the frequency range of 4-50 Hz, E T : Teager energy operator, and SEn: sample entropy. The colored bar shows the mean value of each group, and the error bar shows the corresponding range of ±1 standard deviation. The number within the parentheses is the number of EEG segments, where the outliers have been excluded by the 1.5 inter-quartile range rule. **: Significant difference (p-value < 0.01) was tested by using non-parametric Wilcoxon signed-rank tests.
The group results (mean ± standard deviation) of P b before, during, and after real stimulation were −3.02 ± 5.09, −4.78 ± 4.79, and −2.52 ± 4.39, respectively. Two-way mixed design ANOVA showed that time was a significant factor (F 2,116 = 5.488, p = 0.005), but treatment was not (F 1,58 = 0.043, p = 0.837), and there was a significant interaction between the effects of the two factors (F 2,116 = 4.523, p = 0.013). The group results of E T before, during, and after real stimulation were −2.99 ± 0.53, −2.76 ± 0.21, and −2.61 ± 0.24, respectively. Two-way mixed design ANOVA revealed that time was a significant factor (F 2,106 = 10.581, p < 0.0001), but treatment was not (F 1,53 = 2.83, p = 0.098), and there was a significant interaction between the effects of the two factors (F 2,106 = 9.482, p < 0.001). The group results of SEn before, during, and after real stimulation were 0.33 ± 0.12, 0.62 ± 0.21, and 0.23 ± 0.12, respectively. Two-way mixed design ANOVA showed that both time (F 2,112 = 49.51, p < 0.001) and treatment (F 1,56 = 6.12, p = 0.016) were significant factors, and that there was a significant interaction between the effects of the two factors (F 2,112 = 56.55, p < 0.001).

Performance of the Four Quantitative Indicators in RPS Experiments
Twenty-seven triplets before, during, and after stimulation of EEG segments (5 s) from nine mice receiving RPS were used in the following analysis ( Figure 7). The group results (mean ± standard deviation) of I B/A before, during, and after real stimulation were 1.84 ± 0.58, 5.95 ± 4.29, and 3.82 ± 4.81, respectively. Two-way mixed design ANOVA showed that neither time (F 2,96 = 1.839, p = 0.165) nor treatment (F 1,48 = 3.220, p =0.079) was a significant factor, but that there was a significant interaction between the effects of the two factors (F 2,96 = 6.312, p = 0.003).  In summary, the results of the four indicators were inconsistent during both types of electrical stimulation. Because of the lack of a gold standard, it was unclear which indicator had the best performance. The group results (mean ± standard deviation) of P b before, during, and after real stimulation were −0.53 ± 5.59, −5.72 ± 4.39, and −4.59 ± 7.54, respectively. Two-way mixed design ANOVA showed that time was a significant factor (F 2,102 = 10.85, p < 0.001), but treatment was not (F 1,51 = 0.12, p = 0.75), and there was a significant interaction between the effects of the two factors (F 2,102 = 7.73, p < 0.001). The group results of E T before, during, and after real stimulation were −2.68 ± 0.54, −3.27 ± 0.36, and −3.07 ± 0.62, respectively. Two-way mixed design ANOVA showed that time was a significant factor (F 1,52 = 21.57, p < 0.001), while treatment was not (F 1,52 = 1.677, p = 0.201), and there was a significant interaction between the effects of the two factors (F 2,104 = 14.66, p < 0.001). The group results of Sen before, during, and after real stimulation were 0.29 ± 0.69, 0.31 ± 0.07, and 0.18 ± 0.08, respectively. Two-way mixed design ANOVA showed that time was a significant factor (F 2,92 = 17.90, p < 0.001), while treatment was not (F 1,46 = 0.331, p = 0.568), and there was a significant interaction between the effects of the two factors (F 2,92 = 5.730, p = 0.005).
In summary, the results of the four indicators were inconsistent during both types of electrical stimulation. Because of the lack of a gold standard, it was unclear which indicator had the best performance.

Bifurcation Analysis
Because the results of the performance indicators were inconsistent, we investigated this issue using simulated data. Three sets of one normal and four types of epileptic activities were generated using a seizure model to evaluate the performance indicators, and the results of bifurcation analysis offered an explanation for the discrepancy in indicator performance when using different outputs (Figure 8, right plot). Within the range between the saddle-node (SN, at B/A = 7.824) and saddle-saddle (SS, at B/A = 12.897) bifurcation points, the dynamics of the epilepsy model had one stable and two unstable fixed points, while it had only either one stable or unstable fixed point outside this range. The dynamics underwent supercritical Andronov-Hopf (AH) bifurcation (negative Lyapunov number) at B/A = 2.478 as I B/A increased from below, and then a limit cycle persisted until B/A reached the SN bifurcation point. As shown, the Type II epileptiform activity generated by the epilepsy model corresponded to the dynamics of limit cycles with small amplitudes, and the Type III epileptiform activity was generated by large-amplitude limit cycle dynamics. The irregular Type IV and clustered Type V epileptiform activities appeared when B/A was slightly larger than the SN bifurcation point. The irregular or clustered large spikes in Type IV and Type V epileptiform activities originated from the limit cycle behavior due to escape from the attraction domain of the stable fixed point (lower solid line) in the presence of noise in the model input to the unstable fixed point region (upper dashed lines). Normal activity was generated when B/A was sufficiently far away from the SN bifurcation point, so that the input noise was hardly able to lift z(t) to the unstable region.
The right plot of Figure 8 shows the results when the model output z(t) was used for analysis, and the left plot shows the results when the simulated EEG y(t) was used for analysis. Because of the high-pass filtering property of EEG machines, the DC component was removed from z(t), after which the stable points were all zero. As explained above, Type I EEGs were generated by the smallest B/A ratio, i.e., highest excitability, in four types of epileptiform discharge. The oscillation of Type I EEGs was very fast and small, and may have been obscured by background noise. If the DC component was removed by the high-pass filter of common EEG machines, it became difficult to differentiate the Type I EEGs from normal EEGs, as seen in the left plot of Figure 8. This is a problem of observability.
Type I EEGs were generated by the smallest B/A ratio, i.e., highest excitability, in four types of epileptiform discharge. The oscillation of Type I EEGs was very fast and small, and may have been obscured by background noise. If the DC component was removed by the high-pass filter of common EEG machines, it became difficult to differentiate the Type I EEGs from normal EEGs, as seen in the left plot of Figure 8. This is a problem of observability.  Figure A1).

Discussion
In this study, the performance indicator derived from CSCKF was compared with three conventional performance indicators to evaluate the severity of epileptiform discharges using both simulated and experimental data. For the experimental data, the results of the four indicators were inconsistent during both types of electrical stimulation, although a trend of lower epileptiform discharges during electrical stimulation was observed. Due to the lack of a gold standard, it is unclear which indicator had the best performance. When the simulated EEG segments were used, the results of all four indicators were similar; however, this trend did not match the trend of the parameters of the model. When the model output was used, except for IB/A, the results of the other three indicators  Figure A1).

Discussion
In this study, the performance indicator derived from CSCKF was compared with three conventional performance indicators to evaluate the severity of epileptiform discharges using both simulated and experimental data. For the experimental data, the results of the four indicators were inconsistent during both types of electrical stimulation, although a trend of lower epileptiform discharges during electrical stimulation was observed. Due to the lack of a gold standard, it is unclear which indicator had the best performance. When the simulated EEG segments were used, the results of all four indicators were similar; however, this trend did not match the trend of the parameters of the model. When the model output was used, except for I B/A , the results of the other three indicators were almost identical. For I B/A , the estimates faithfully reflected the severity of epileptiform discharges.
The detection and quantification of seizure severity using electrophysiological measures is an important topic, and many signal processing techniques have been developed [38]. The main difficulty is the non-linearity and non-stationarity of epileptiform discharges. A simple and robust method to overcome the non-linearity and non-stationarity is therefore needed. Another issue is the lack of a gold standard. Conventionally, the severity of epileptiform discharge is assessed by experts through visual inspection. Most of the studies intuitively used spike counts per unit time (frequency equivalent) [24] or amplitude [25] as the quantifier without explicit evidence of its fidelity. When frequency is relatively stable, amplitude is proportional to band power (P b ). Teager energy operator (E T ) takes both frequency and amplitude into consideration. Entropy (SEn) is a nonlinear indicator that measures randomness, based on the fact that regularity increases during seizure attack [39]. However, there may not be consensus among experts, when frequency and amplitude change in opposite direction, as seen with the Type I and Type II seizure activities in Figure A1. From the formulae, P b is the band power of a signal between 4 and 50 Hz. The more EEG spikes are present and the greater the amplitude of each spike, the higher the P b value. However, if the frequency of the EEG spike is out of the defined range of P b , then this description may not be true. E T is another method of calculating signal energy, and it can be used to quantify instantaneous changes in signal amplitude and frequency. The indicator E T used in this study is the average value of E T over a period of time, and it was also positively correlated with the amplitude and frequency of the EEG spikes. Sample entropy is a nonlinear metric used to quantify the regularity or complexity of a signal. In general, fewer variations in the spike frequency and amplitude will lead to lower sample entropy of the EEG signal. In addition, for EEG signals with similar spike frequencies, the general amplitude of the waveform has little effect on their sample entropy values. As seen in Figure A1, although there were clear differences in the general amplitudes between the three Type II signals, their SEn values were nearly identical. Therefore, special care should be taken when using sample entropy as an indicator of epileptiform discharge severity. In this paper we used the ratio of inhibitory weight to excitatory weight in the local neural network as the indicator, which is used as an indicator in a study focused on postictal generalized EEG suppression [40]. The advantage is that it is more theoretically based, while the disadvantage is that it is difficult to estimate.
RNS suppressed epileptiform discharges only during stimulation, while the suppressive effects of RPS persisted after the stimulation was stopped, as reveal by P B and E T . E T showed largest difference between RNS and RPS. However, if the post stimulation suppression was used as the indicator, RNS was superior to RPS. We think the contradiction may reflect the different suppressive mechanisms underlying these two types of stimulation. As is revealed by Figure 9, the amplitude of electrical activity decreased while the mean frequency unchanged, implying RNS partially reduced synchronization but the seizure loop persisted. On the other hand, RPS directly disrupted seizure loop and the main spike frequency disappeared. RPS possibly inhibits a spike train by maintaining a post-spike refractory period. The mechanism of RNS is still under active research. As mentioned in the Introduction, the effects of RNS may depend on the amplitude of RNS. At a larger amplitude, RNS may take over the neural circuit, exciting the neuron population randomly and interrupting the regularity of spikes. . The group amplitude and mean frequency of EEGs before and during the two stimulation protocols, respectively. Significant difference (p-value *: < 0.05, **: < 0.01 and ***: < 0.001) was tested by using non-parametric Wilcoxon signed-rank tests.
The results of bifurcation analysis imply that it may be feasible to estimate the excitability of a neural network and quantify the severity of epileptiform discharges if DC EEG machines are used to record neural activities. However, the problem of more artefacts when using DC EEG machines has to be overcome. Further studies are needed to investigate this issue and validate our findings.

Conclusions
Based on the artificially generated and experimentally derived EEG data, the four performance indicators showed similar but varying trends, partially conforming to the severity of seizure activity. However, when the DC component of EEG was preserved as in the simulated data, only CSCKF could quantify the severity of seizure activity.  . The group amplitude and mean frequency of EEGs before and during the two stimulation protocols, respectively. Significant difference (p-value *: < 0.05, **: < 0.01 and ***: < 0.001) was tested by using non-parametric Wilcoxon signed-rank tests.
The results of bifurcation analysis imply that it may be feasible to estimate the excitability of a neural network and quantify the severity of epileptiform discharges if DC EEG machines are used to record neural activities. However, the problem of more artefacts when using DC EEG machines has to be overcome. Further studies are needed to investigate this issue and validate our findings.

Conclusions
Based on the artificially generated and experimentally derived EEG data, the four performance indicators showed similar but varying trends, partially conforming to the severity of seizure activity. However, when the DC component of EEG was preserved as in the simulated data, only CSCKF could quantify the severity of seizure activity. Funding: This work was partly funded by a grant (MOST 110-2314-B-006 -080) from the Ministry of Science and Technology, Taiwan. The funding source was not involved in conceptualization, study design, data collection, analysis or interpretation, manuscript writing or journal submission.

Institutional Review Board Statement:
The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Institutional Review Board of National Cheng Kung University (IACUC Approval No: 108223, date: 3 May 2019).

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to its size.

Acknowledgments:
The data that support the findings of this study are available on reasonable request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest. The work described has not been published previously.  Appendix B Table A1. Parameters for the mathematical seizure model of hippocampus in mice.  Appendix B Table A1. Parameters for the mathematical seizure model of hippocampus in mice.