Proposing a two-level stochastic model for epileptic seizure genesis

Abstract

By assuming the brain as a multi-stable system, different scenarios have been introduced for transition from normal to epileptic state. But, the path through which this transition occurs is under debate. In this paper a stochastic model for seizure genesis is presented that is consistent with all scenarios: a two-level spontaneous seizure generation model is proposed in which, in its first level the behavior of physiological parameters is modeled with a stochastic process. The focus is on some physiological parameters that are essential in simulating different activities of ElectroEncephaloGram (EEG), i.e., excitatory and inhibitory synaptic gains of neuronal populations. There are many depth-EEG models in which excitatory and inhibitory synaptic gains are the adjustable parameters. Using one of these models at the second level, our proposed seizure generator is complete. The suggested stochastic model of first level is a hidden Markov process whose transition matrices are obtained through analyzing the real parameter sequences of a seizure onset area. These real parameter sequences are estimated from real depth-EEG signals via applying a parameter identification algorithm. In this paper both short-term and long-term validations of the proposed model are done. The long-term synthetic depth-EEG signals simulated by this model can be taken as a suitable tool for comparing different seizure prediction algorithms.

Appendices

Appendix A. Details of the depth-EEG model

The structure of depth-EEG model is shown in Fig. A1. As shown in this figure, in the depth-EEG model each synaptic process (excitatory, fast and slow inhibitory processes) is a combination of a linear system (pulse-to-wave conversion at a synaptic level) and a nonlinear function (wave-to-pulse operator at the soma of neurons, i.e., $S(v)=2{e}_0/\left(1+{e}^{r\left({v}_0-v\right)}\right)$), where v is the average potential of the pre-synaptic cells and S(v) is the mean firing rate of the post synaptic cells. This process is expressed by the following differential equations:

Fig. A1

Organization of the depth-EEG model including four neuronal subsets: pyramidal cells, excitatory interneurons, dendritic projecting interneurons with slow synaptic kinetics and somatic-projecting interneurons (grey rectangle) with faster synaptic kinetics, x(t) denotes for the step function. (from (Wendling et al. 2002))

$$ \begin{array}{l}{\dot{y}}_0\left(\mathrm{t}\right)={y}_5\left(\mathrm{t}\right)\\ {}{\dot{y}}_5\left(\mathrm{t}\right)= Aa S\left[{y}_1\left(\mathrm{t}\right)-{y}_2\left(\mathrm{t}\right)-{y}_3\left(\mathrm{t}\right)\right]-2\tau {y}_5\left(\mathrm{t}\right)-{\tau}^2{y}_0\left(\mathrm{t}\right)\\ {}{\dot{y}}_1\left(\mathrm{t}\right)={y}_6\left(\mathrm{t}\right)\\ {}{\dot{y}}_6\left(\mathrm{t}\right)= Aa\left\{x\left(\mathrm{t}\right)+{C}_2S\left[{C}_1{y}_0\left(\mathrm{t}\right)\right]\right\}-2a{y}_6\left(\mathrm{t}\right)-{a}^2{y}_1\left(\mathrm{t}\right)\\ {}{\dot{y}}_2\left(\mathrm{t}\right)={y}_7\left(\mathrm{t}\right)\\ {}{\dot{y}}_7\left(\mathrm{t}\right)= Bb{C}_4S\left[{C}_3{y}_0\left(\mathrm{t}\right)\right]-2b{y}_7\left(\mathrm{t}\right)-{b}^2{y}_2(t)\\ {}{\dot{y}}_3\left(\mathrm{t}\right)={y}_8\left(\mathrm{t}\right)\\ {}{\dot{y}}_8= Gg{C}_7S\left[{C}_5{y}_0\left(\mathrm{t}\right)-{C}_6{y}_4\left(\mathrm{t}\right)\right]-2g{y}_3\left(\mathrm{t}\right)-{g}^2{y}_3\left(\mathrm{t}\right)\\ {}{\dot{y}}_4\left(\mathrm{t}\right)={y}_9\left(\mathrm{t}\right)\\ {}{\dot{y}}_9= Bb S\left[{C}_3{y}_0\left(\mathrm{t}\right)\right]-2b{y}_4\left(\mathrm{t}\right)-{b}^2{y}_4\left(\mathrm{t}\right)\\ {}{y}_{out}\left(\mathrm{t}\right)={y}_1\left(\mathrm{t}\right)-{y}_2\left(\mathrm{t}\right)-{y}_3\left(\mathrm{t}\right)\end{array} $$

(A1)

Interpretation of linear system parameters (C ₁ to C ₇ , a, b, g, A, B, G, τ) and nonlinear function parameters (e ₀ , v ₀ and r) and their standard values are presented in Table A1 based on (Wendling et al. 2002). The adjustable variables of the depth-EEG model are A, B, G, and τ, while other parameters are taken to be constant.

Table A1 Depth-EEG Model parameters, interpretation and standard values used to produce background EEG activity (Wendling et al. 2002)

In (A1) the signals named y ₀ to y₉ are the state signals underlying each subset of the hippocampus. According to Fig. A1, y₀ to y₄ are post-synaptic potentials of each subset of the model. y _out is the model output and x is the model input corresponding to a white noise signal for the depth-EEG model of a single area.

Differential equations (A1) should be solved by a stochastic numerical algorithm; but, like (David and Friston, 2003) for the sake of simplicity, we opted for solving the equations by Euler algorithm.

Appendix B. The parameter tracking algorithm

To track the variation of parameters during the long-term depth-EEG signals, the optimization problem (3) should be solved consecutively for short windows of the signal.

For the first window of the signal there is no a priori knowledge about the amount of the parameter vector. In other words, the feasible parameter space of the first window of signal is as follows (this space is achieved through discretization):

$$ {\varOmega}_1=\left\{\begin{array}{l}\left.\left[ i\varDelta, j\varDelta, k\varDelta, l\varDelta \right]\right|\\ {}\kern3.25em i=1:R,j=1:R,\\ {}\kern3.25em k=1:R,l=1:R\end{array}\right\} $$

(B1)

But, after parameter identification of the first window, for other windows the full search is performed in the vicinity of the parameters obtained for their preceding windows. This ensures the extraction procedure of the parameter sequences underlying depth-EEG signals are performed efficiently in time. This vicinity is determined as follows: if for the tth window of signal the parameter θ _t  = [A _t , B _t , G _t , τ _t ] = [IΔ, JΔ, kΔ, LΔ] has been identified, the feasible parameter space for the t + 1th window is:

$$ {\varOmega}_{t+1}=\left\{\begin{array}{l}\left.\left[ i\varDelta, j\varDelta, k\varDelta, l\varDelta \right]\right|\\ {}\kern3.25em i=I-1:I+1,j=J-1:J+1,\\ {}\kern3.25em k=K-1:K+1,l=L-1:L+1\\ {}\kern2.75em and\ {F}_{onset}\left(\left[ i\varDelta, j\varDelta, k\varDelta, l\varDelta \right],{W}^{*}\right)< th\end{array}\right\} $$

(B2)

where W ^* is the optimum weight vector obtained for the considered parameter vector. We set th = 10. However, when the type of EEG activity changes (e.g., transition from normal state to the sustained discharge of spikes), a sudden jump in the parameter sequence may occur, such that the feasible parameter space Ω _t + 1 becomes a null set. In these situations the search space of parameter vectors must be reset again, i.e. Ω _t + 1 = Ω ₁.

Although according to this procedure, which assumes the continuity of the parameters, the optimum parameter sequence may not be obtained, but a smoother sub-optimal parameter sequence has been concluded in much shorter time.