Seizure tracking of epileptic EEGs using a model-driven approach

In 2000, a neurophysiological finding revealed that two types of ${\rm GABA}_A$ synaptic responses, that is, a fast one near the soma and a slow one in the dendrites, may play crucial roles in the formation of gamma rhythms in pyramidal cells [21]. Two such types of responses come from two types of GABAergic inhibition interneurons, which are dendritic inhibition interneurons (denoted by ${\rm GABA}_{A,slow}$) and somatic inhibition interneurons (denoted by ${\rm GABA}_{A,fast}$) respectively. Originated from this finding, Wendling et al [20] proposed a NMM to study the generation of high-frequency EEG activities in 2002, where the inhibition interneuron population was replaced by considering two populations (that is, ${\rm GABA}_{A,fast}$ and ${\rm GABA}_{A,slow}$) in the existing model proposed by Jansen and Rit [22]. For convenience, we name this model the 'Wendling model' in what follows. Figure 1 illustrates the topological structure of the Wendling model. Four populations in the model are included: the pyramidal population (PY), the excitatory interneuron population (E), the slow dendritic inhibitory interneuron population (I1) and the fast somatic inhibitory interneuron population (I2). Among four populations, PY receives a postsynaptic firing rate from E, I1 and I2, as well as send postsynaptic potential (PSP) to E, I1 and I2. Meanwhile, population I2 also receives a postsynaptic firing rate from I1 due to the physiological connection between them.

Zoom In Zoom Out Reset image size

In each population, the dynamics evolution is formulated by two operators. The first one is a convolution operator where the firing rate $m(t)$ is mapped into the PSP $v(t)$, that is,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v(t) = h(t) \otimes m(t). \label{WEQ} \nonumber \end{align} \tag{ 1 }$$

In equation (1), the response function $h(t)$ corresponds to three cases: excitatory h_e, slow inhibitory h_i, fast inhibitory h_f, which are denoted by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_e(t)= Aate^{-at}, t\geqslant 0, \nonumber \end{align} \tag{ 2 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_i(t)= Bbte^{-bt}, t\geqslant 0, \nonumber \end{align} \tag{ 3 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_f(t)= Ggte^{-gt}, t\geqslant 0. \label{h(t)} \nonumber \end{align} \tag{ 4 }$$

Here, A, B and G denote the synaptic gains, a, b and g are the lumped representation of time constants.

Conversely, the second operator transforms the PSP $v(t)$ into the firing rate $m(t)$, which is always selected to be a sigmoid function. It is generally described by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S(v(t))=\frac{2e_0}{1+e^{r(v_0-v(t))}} \label{sig} \nonumber \end{align} \tag{ 5 }$$

where e₀ is the maximum firing rate, $v_0$ is the firing threshold and r is the steepness of the sigmoidal transformation.

Since it is difficult to calculate the operator '$\otimes$' in a direct way, the Laplace transformation is applied so that equation (1) could be re-expressed by a pair of first-order ordinary differential equations equivalently, that is,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \dot{v}(t) &= v_1(t), \nonumber \\ \dot{v_1}(t)&= \Gamma \gamma m(t)- 2\gamma v_1(t)-\gamma^2v(t). \label{sod} \nonumber \end{align} \tag{ 6 }$$

Here, $\Gamma \in \{A,B,G\}$ and $\gamma \in \{a, b, g\}$. Figure 2(a) illustrates the computational block in each population.

Zoom In Zoom Out Reset image size

The Wendling model broke new ground in the field of NMM due to the successful design by considering two types of GABAergic inhibition interneurons, that is, ${\rm GABA}_{A,fast}$ interneurons (I2) and ${\rm GABA}_{A,slow}$ interneurons (I1). However, there is a lack concerning the response alteration of population PY when it interacts with I1 and I2 respectively. Meanwhile, with the recognition that time delay in neuronal signal transmission is a key factor to cause seizures [23], there is a lack in concerning the response alteration of populations E and I1, which are influenced by the entry of time delays in connections. For the sake of keeping modeling more flexibile in producing epileptiform activities, the computational blocks in populations should be re-formulated mathematically under new considerations. In addition, the neurons in the same type of population may have different responses for the same stimulus [19], which means various sub-populations with their own kinetics should be included in one population. And it will be more effective in capturing the diversity of neuronal networks. Motivated by such recognition, a new time-delay Wendling model with sub-populations (TD-W-SP model) will be proposed in this paper.

We first consider the case of pyramidal population, which can be regarded as a central hub in the construction of NMMs. In the original Wendling model, it accounts for the same response function for population PY when it interacts with populations I1 and I2 respectively, which means PY sends PSP into I1 and I2 in the same way. However, two situations should be addressed on PY in light of considering the different characteristics of ${\rm GABA}_{A,slow}$ and ${\rm GABA}_{A,fast}$. Specifically, we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_e^{\rm PY}(t) = \left\{\begin{array}{@{}ll@{}} Aate^{-at}, & {\rm if~PY~\rightarrow~E/I1}; \nonumber \\ A'a'te^{-a't}, & {\rm if~PY~\rightarrow~I2}. \end{array} \right. \label{MPY} \nonumber \end{align} \tag{ 7 }$$

Here, $h_e^{\rm PY}(t)$ represents the response function of PY in the TD-W-SP model with different parameters $(A,a)$ and $(A',a')$ under two situations. Figure 2(b) illustrates the computational block of PY with new formulation.

Next, we consider the case of interneuron population. In [23], Geng et al added an additional time-delay in the inhibitory feedback loop, and to simulate several different types of EEG activity including alpha wave and ictal EEG. This work succeeded in showing that the time delay in neuronal signal transmission is a key factor that causes seizures. Originating from this work, in order to further account for the impaired balance between excitatory and inhibitory populations, the time delay $\tau_a$ and $\tau_b$ are added both in the excitatory and inhibitory feedback loop. Based on these points, we formulate the computational blocks with time delay in E and I1 as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_e^E(t)= \delta(t-\tau_a)Aate^{-at}, \label{TDE} \nonumber \end{align} \tag{ 8 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_i(t)= \delta(t-\tau_b)Bbte^{-bt}, \label{TDI} \nonumber \end{align} \tag{ 9 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle t \geqslant 0 , \tau\geqslant 0 \nonumber \end{align} \tag{ 10 }$$

where $\tau_a$ and $\tau_b$ denote the time delays existing in the connections between PY and E, as well as PY and I1 respectively, $h_e^E(t)$ represents the excitatory response of E in the TD-W-SP model. Figure 2(c) illustrates the time-delay computational block.

Finally, according to [4], David et al revealed that the neurons in the same type of population may have different responses for the same stimulus, and verified that the Jansen & Rit model could capture more complex dynamics of a neural system if each population is extended by N parallel subpopulations with different kinetics. Motivated by such recognition, in this work, we divided each population into N sub-populations deploying in parallel with different kinetics, whose structure is shown in figure 2(d). Assume that $v_1, v_2, \cdots, v_N$ are the outputs and $w_1, w_2, \cdots, w_N$ the weights of sub-populations, then the final output of the population is formulated to be the weighted sum of them, that is,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v= w_1 v_1 + w_2 v_2 + \cdots + w_N v_N. \label{SUB} \nonumber \end{align} \tag{ 11 }$$

By combining the above modifications (say, equations (7)–(11)) with the Wendling model, the mathematical formulation of the TD-W-SP model in the case of N  =  2 can be represented by equation (12) (see page 5). Figure 3 illustrates the neuronal population organization of the TD-W-SP model. The final model output is the weighted sum of membrane potential of two excitatory sub-populations and four inhibitory sub-populations $wv_2+(1-w)v_8-wv_3-(1-w)v_9-wv_6-(1-w)v_{12}$.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{MathModel} \ddot{v}_1(t)&= AaS(wv_2(t-\tau_a)+(1-w)v_8(t-\tau_a)-wv_3(t-\tau_b)-(1-w)v_9(t-\tau_b)-wv_6(t)\nonumber \\ &\quad-(1-w)v_{12}(t))-2a\dot{v_1}(t)-a^2v_1(t);\nonumber \\ \ddot{v}_2(t)&= AaC_2S(C1(wv_1(t)+(1-w)v_7(t)))-2a\dot{v_2}(t)-a^2v_2(t)+Aap(t);\nonumber \\ \ddot{v}_3(t)&= BbC_4S(C_3(wv_1(t)+(1-w)v_7(t)))-2b\dot{v_3}(t)-b^2v_3(t);\nonumber \\ \ddot{v}_4(t)&= A'a'S(wv_2(t-\tau_a)+(1-w)v_8(t-\tau_a)-wv_3(t-\tau_b)-(1-w)v_9(t-\tau_b)-wv_6(t)\nonumber \\ &\quad-(1-w)v_{12}(t))-2a'\dot{v_4}(t)-a'^2v_4(t);\nonumber \\ \ddot{v}_5(t)&= BbC_6S(C_3(wv_1(t)+(1-w)v_7(t)))-2b\dot{v_5}(t)-b^2v_5(t);\nonumber \\ \ddot{v}_6(t)&= GgC_7S(wv_4(t)+(1-w)v_{10}(t)-wv_5(t)-(1-w)v_{11}(t))-2g\dot{v_6}(t)-g^2v_6(t);\nonumber \\ \ddot{v}_7(t)&= A_1a_1S(wv_2(t-\tau_a)+(1-w)v_8(t-\tau_a)-wv_3(t-\tau_b)-(1-w)v_9(t-\tau_b)-wv_6(t)\nonumber \\ &\quad-(1-w)v_{12}(t))-2a_1\dot{v_7}(t)-a_1^2v_7(t);\nonumber \\ \ddot{v}_8(t)&= A_1a_1C_2S(C1(wv_1(t)+(1-w)v_7(t)))-2a_1\dot{v_8}(t)-a_1^2v_8(t)+A_1a_1p(t);\nonumber \\ \ddot{v}_9(t)&= B_1b_1C_4S(C_3(wv_1(t)+(1-w)v_7(t)))-2b_1\dot{v_9}(t)-b_1^2v_9(t);\nonumber \\ \ddot{v}_{10}(t)&= A_1^{\prime}a_1^{\prime}S(wv_2(t-\tau_a)+(1-w)v_8(t-\tau_a)-wv_3(t-\tau_b)-(1-w)v_9(t-\tau_b)-wv_6(t)\nonumber \\ &\quad -(1-w)v_{12}(t))-2a_1^{\prime}\dot{v}_{10}(t)-a_1^{'2}v_{10}(t);\nonumber \\ \ddot{v}_{11}(t)&= B_1b_1C_6S(C_3(wv_1(t)+(1-w)v_7(t)))-2b_1\dot{v}_{11}(t)-b_1^2v_{11}(t);\nonumber \\ \ddot{v}_{12}(t)&= G_1g_1C_7S(wv_4(t)+(1-w)v_{10}(t)-wv_5(t)-(1-w)v_{11}(t))-2g_1\dot{v}_{12}(t)-g_1^2v_{12}(t). \nonumber \end{align} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

In this part, we focus on analyzing the effectiveness of the proposed TD-W-SP model in simulating seizure-like EEGs. Figures 7(a)–(c) illustrate the outputs of three different models, which are the original Wendling model, the Wendling model with the improved response function $h_e^{\rm PY}$ (abbreviated by 'Wendling  +  IRF-PY'), and the Wendling model with $h_e^{\rm PY}$ as well as the time delay $\tau_a, \tau_b$ (abbreviated by 'Wendling  +  IRF-PY  +  TD'). Figures 7(d)–(f) illustrates three real depth-EEG signals recorded in the human hippocampus (during SEEG exploration and using intracerebral multiple lead depth-electrodes) in [20]. We can observe from figure 7 that: (1) the output of the Wendling model is similar to normal background activity when all parameters are set to be their nominal values shown in table 2 (comparing subfigures (a) and (d)); (2) the rhythmic discharges appear when the new response function $h_e^{\rm PY}$ is added. Here, all parameters are nominal values except that $A'=4.25$ (comparing subfigures (b) and (d)). It resembles the slow quasi-sinusoidal activity from the morphology and rhythmicity point of view, which may be encountered during seizures [20]; (3) the seizure-like EEG signal occurs when delay factors $\tau_a, \tau_b$ are additionally included when all parameters are nominal values except that $A'=4.25, \tau_a=0.1~{\rm ms}, \tau_b=2~{\rm ms}$. It is analogous to the sustained discharge of spikes, which may be encountered at seizure onset or during seizures [20].

Zoom In Zoom Out Reset image size

Table 2. The nominal value of parameters in the TD-W-SP model.

Parameter	Nominal value	Unit
$A,A'$	3.25	mV
B	22	mV
G	10	mV
$a,a'$	100	s−1
b	50	s−1
g	500	s−1
C	135	/
$v_0$	6	mV
e0	2.5	s−1
$v_0$	0.56	mV−1

In order to further verify that the model with sub-populations is good at capturing the diversity of neuronal networks, we compare the results obtained by the original Wendling model and the Wendling model with sub-populations (abbreviated by 'Wendling  +  SP'). As shown in [20], the original Wendling model could produce two types of waveforms with a different combination of parameters $(B, G) \in [0,50] \times [0,30]$ and A  =  3, which are simulated normal background activity and slow rhythmic activity (i.e. 'type 1' and 'type 2' in figure 8(a)). However, we can observe that the 'Wendling  +  SP' could produce four types of waveforms under the same assumptions of parameters (i.e. all sub-figures 'type 1- type 4' in figure 8(a)), where type 3 corresponds to the simulated sustained discharge of spikes and type 4 corresponds to slow quasi-sinusoidal activity. These results validate that the Wendling model with sub-populations could not only produce more types of waveforms, but also the epileptic waveforms are included in them. Figure 8(b) illustrates the parameter distribution of $(B, G)$ corresponding to four types of EEGs.

Zoom In Zoom Out Reset image size

In this part, we focus on verifying the feasibility and efficiency of MTM-F-GA on artificial EEG data. Here, the time delay $\tau$ is taken as an example to show the results. A sequence of $\tau$ is firstly generated based on uniform distribution, which is shown as the red line in figure 9. Given the generated sequence of $\tau$ and other parameters with nominal values, we can then calculate a sequence of outputs of the TD-W-SP model by solving equation (12). Each output is regarded as an artificial EEG here. Next, we apply the method MTM-F-GA into the artificial EEGs, the optimal values of $\tau$ can be obtained, which are shown as the blue dotted line in figure 9. Obviously, there are only tiny differences between the real values and optimized values of $\tau$, which shows the good performance of the method MTM-F-GA. Note that similar results can be obtained for other parameters.

Zoom In Zoom Out Reset image size

In this part, we will verify (1) why it is reasonable to consider parameters $A, B, \tau$, and (2) how parameters $A, B, \tau$ work, in the aspect of exploring the evolution of 'Preictal-ictal-Postictal' in epileptic EEGs.

We first validate the capability of the proposed model-driven approach to recover and track parameters of interest ($A, B, \tau_a, \tau_b$) using the simulated data. Firstly, a 30 s-long simulated EEG segment is generated by the proposed model TD-W-SP (see figure 10(a)), where the normal-like EEG (first 10 s), slow rhythmic activity (middle 10 s) as well as seizure-like EEG (last 10 s) are obtained respectively with different parameter settings. Note that the parameter values for each case are known as prior and served as a ground truth, which are illustrated as the black solid lines in figures 10(b)–(e). Furthermore, the proposed model is trained by the simulated EEG to obtain the estimated parameter values, whose results are illustrated as the blue dotted line in figures 10(b)–(e). We can observe from figure 10 that the proposed model-driven approach could better recover and track the parameters on simulated data.

Zoom In Zoom Out Reset image size

Next, we verify that parameters $A, B, \tau$ are helpful in exploring the evolution of epileptic seizures. Given an EEG signal including three stages of a seizure (30 s pre-ictal, 20 s ictal and 30 s post-ictal), which are labeled by experts (see figure 11(a)), we segment it into a number of EEG segments by 2 s sliding windows with 1 s overlap. Then for each EEG segment, we can obtain the optimal values of parameters $A, B, \tau, a$ by equations (13) and (14) respectively. The variation trend of all parameters are illustrated in figures 11(b)–(e). As seen from figure 11(b), we can observe that a rapid increase of A around 15 s happens, and 25 s later, the values of A tend to be parallel after a rapid decrease. Similar results can be obtained in the cases of B and $\tau$ if we change the increase and decrease mutually (see figures 11(c) and (d)). However, such an apparent trend has not been consistently observed in the case of a, as shown in figure 11(e). These results verify that the variation trends of $A, B, \tau$ are very helpful in guiding us to explore the evolution of 'Preictal-ictal-Postictal' in epileptic EEGs.

Zoom In Zoom Out Reset image size

Finally, we verify the effectiveness of proposed I, TP-PI and TP-IP in the field of tracking different periods of a seizure. In order to calculate the threshold, we first select one EEG file without any seizure for each patient, and then artificially choose 200 s-long continuous EEG segment without a single spike or other epileptiform waves. Such obtained 200 s-long EEG segments are then defined as the 'normal EEG signal' for each patient. Then 100 normal (i.e. K_nor  =  100) EEG segments are obtained by a 2 s sliding window without any overlap. In our simulation, k1  =  2. Figure 12(a) illustrates an 80 s-long EEG segment including 30 s pre-ictal, 20 s ictal and 30 s post-ictal EEGs from 'patient 1' and the sequence of its indexes. The two thresholds are 0.4 and 0.79 respectively. Then the two turning points are 2982 s and 3020 s, which are 14 s earlier than the labeled onset time (2996 s) and 4 s later than the labeled offset time (3016 s). Figure 12(b) shows a 112 s-long EEG segment (30 s pre-ictal, 52 s ictal and 30 s post-ictal EEGs) from 'patient 3' and corresponding sequence of indexes. The TP-PI and TP-IP are 352 s and 415 s, which are 10 s earlier than the labeled onset time (362 s) and 1 s later than the labeled offset time (414 s). Also, we can observe from figure 12, the trend of indexes clearly reflect the cycle of seizures. On the other hand, we can regard the estimated TP-PI as an alarming time for early seizure detection. It is known that the more early the predicted onset time, the more timely the pre-control could be carried out. However, the more early the predicted offset time, the more difficult the information of seizures could be captured. For each seizure, we apply the pre-estimated duration (PRD) and the post-estimated duration (POD) to evaluate the performance of the proposed method. Due to the limited space, we only list the results of one seizure for each patient (see table 3). Table 3 shows the ictal duration whose left end is the onset time labeled by experts (LON) and the right end is the offset time labeled by experts (LOF), the ${\rm TP}_{\rm PI}$ and ${\rm TP}_{\rm IP}$ estimated by equations (18) and (19), the PRD and the POD respectively. We can observe from table 3 that the estimated ${\rm TP}_{\rm PI}$ generally present around 10 s earlier than the labeled onset times, and the estimated ${\rm TP}_{\rm IP}$ are 1 s–4 s later than the labeled offset times. For each patient, we define the 'tracking accuracy (TA)' to evaluate the performance of the proposed method, denoted by

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle {\rm TA} &= \frac{\#{\{{\rm seizures~satisfying}\ {\rm PRD} > 0 \& {\rm POD}>0\}}}{\#\{{\rm all~seizures~of~each~patient}\}},\nonumber \\ {\rm PRD} & = {\rm LON} - {\rm TP}_{\rm PI},\nonumber \\ {\rm POD} & = {\rm TP}_{\rm IP} - {\rm LOF}. \nonumber \end{align*}$$

Zoom In Zoom Out Reset image size

Table 3. Performances of the model-driven seizure tracking method.

Patient	labeled ictal duration (s)	${\rm TP}_{\rm PI}$ (s)	${\rm TP}_{\rm IP}$ (s)	PRD (s)	POD (s)	TA
1	(2996, 3016)	2982	3020	14	4	0.86
2	(3369, 3378)	3367	3381	3	3	1
3	(362, 414)	352	415	10	1	0.71
4	(7804, 7853)	7796	7857	8	4	1
5	(1086, 1196)	1077	1198	11	2	0.8
6	(2670, 2841)	2657	2842	13	1	1
7	(12 231, 12 295)	12 228	12 298	10	3	0.75
8	(6313, 6348)	6304	6349	9	1	0.86

Here, $\#{\{{\rm seizures~satisfying}\ {\rm PRD} > 0 \& {\rm POD}>0\}}$ is the number of seizures satisfying ${\rm PRD} > 0$ and ${\rm POD}>0$, which means TP-PI is earlier than LON and TP-IP is later than LOF (that is, the estimated ictal duration totally covers the labeled ictal duration). The corresponding results for each patient are listed in table 3.