A model-based method for computation of correlation dimension, Lyapunov exponents and synchronization from depth-EEG signals

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

journal homepage: www.intl.elsevierhealth.com/journals/cmpb

A model-based method for computation of correlation dimension, Lyapunov exponents and synchronization from depth-EEG signals F. Shayegh a,∗ , S. Sadri a , R. Amirfattahi a , K. Ansari-Asl b a

Digital signal Processing Lab, Department of Electrical and Computer Engineering, Isfahan University of Technology, 84156-83111 Isfahan, Iran b Electrical Department, Faculty of Engineering, Shahid Chamran University of Ahvaz, Ahvaz, Iran

a r t i c l e

i n f o

a b s t r a c t

Article history:

In order to predict epileptic seizures many precursory features, extracted from the EEG sig-

Received 9 April 2013

nals, have been introduced. Before checking out the performance of features in detection

Received in revised form

of pre-seizure state, it is required to see whether these features are accurately extracted.

14 August 2013

Evaluation of feature estimation methods has been less considered, mainly due to the lack

Accepted 28 August 2013

of a ground truth for the real EEG signals’ features. In this paper, some simulated longterm depth-EEG signals, with known state spaces, are generated via a realistic neural mass

Keywords:

model with physiological parameters. Thanks to the known ground truth of these synthetic

Depth-EEG generator

signals, they are suitable for evaluating different algorithms used to extract the features.

Synchronization

It is shown that conventional methods of estimating correlation dimension, the largest

Largest Lyapunov exponent

Lyapunov exponent, and phase coherence have non-negligible errors. Then, a parameter

Correlation dimension

identification-based method is introduced for estimating the features, which leads to bet-

Seizure

ter estimation results for synthetic signals. It is shown that the neural mass model is able

Accurate feature extraction

to reproduce real depth-EEG signals accurately; thus, assuming this model underlying real depth-EEG signals, can improve the accuracy of features’ estimation. © 2013 Elsevier Ireland Ltd. All rights reserved.

1.

Introduction

Epileptic patients constitute about 1% of the world’s population [10]. Life quality of epileptic patients depends on the use of anti-epileptic drugs. Long-term usage of these drugs may create unavoidable side effects. Since seizures occur without any obvious herald, the possibility of predicting seizures will open new horizons of remedies to prevent, or at least control, seizures. Reduction of drug dosage, confining the use of drugs to emergencies, stimulation of the vague nerve, and

other techniques will improve the quality of patients’ lives by reducing their intense feelings of helplessness. Different algorithms have been used to assess whether seizures can be anticipated or not. These algorithms involve extracting some characteristic features from signals, consisting of both uni-variate features, like correlation dimension [9,24] and Lyapunov exponents [15], and bi-variate features, like synchronization [32] and dynamical entrainment [16–18]. Three main features that have dominated the seizure prediction field during the last years are correlation dimension, largest Lyapunov exponent, and synchronization.

∗

Corresponding author. Tel.: +98 9133173220. E-mail addresses: [email protected], [email protected] [email protected] (R. Amirfattahi), [email protected] (K. Ansari-Asl). 0169-2607/$ – see front matter © 2013 Elsevier Ireland Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cmpb.2013.08.014

(F.

Shayegh),

[email protected]

(S.

Sadri),

324

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

Elger and Lehnertz [9,24] applied the correlation dimension, as a measure of dynamical complexity, to discriminate between inter-ictal and pre-ictal states. But, later, by some studies [2,14] the ability of the correlation dimension to predict seizures was questioned. For example, data by which Martinerie and his colleagues [28] distinguished pre-ictal and inter-ictal states, via correlation density, were re-evaluated by [30], who showed that this measure reflects only the variance of EEG signals. In fact, it was suspected if the presence of non-linearity in a signal can justify the use of non-linear and complicated measures to characterize dynamical changes (see also [36]). It is important to prove that whether these complicated methods indeed outperform simpler linear measures or not (see also [4]). However, correlation dimension is still considered for seizure prediction, e.g., in combination to other features [1,31]. On the other hand, since seizure generation is attributed to synchronous firing of neurons, something like regularity of neural networks; so, by taking the complex brain as a chaotic system, one would think of seizure as a state of decrement in chaotic behavior [15] estimated the largest Lyapunov exponent as an indicator of this phenomenon. They showed a decrease in chaoticity of the brain some minutes before an epileptic seizure. Then, studies by Lai et al. [22,23] raised doubts about the suitability of the Lyapunov exponent for seizure prediction. Furthermore, Mormann et al. [32] considered the variation of synchronization between different brain areas before and during seizure. They measured phase synchronization between signals recorded from two channels, and observed such patterns of phase synchronization before seizure that were not found in seizure-free recordings. Le Van Quyen et al. [25,26] confirmed these findings through the study on eight patients suffering from neocortical epilepsy. Then, in the controlled studies comprising defined groups of patients with pre-ictal and inter-ictal control recordings, phase synchronization [33,34] was shown to be suitable for distinguishing between inter-ictal and pre-ictal data. In 2003, Chavez et al. [5] claimed that pre-ictal changes in phase synchronization occur predominantly in the beta band. Also, pre-ictal synchronization changes were found to be locally, restricted to specific channels, rather than a global phenomenon. Although using the concept of seizure time surrogates [3] bi- and multivariate measures were reported to have better performance than uni-variate measures [17,27,35], power of synchronization measures in seizure prediction is still under debate [11]. Briefly, none of these features can predict seizures accurately enough to be applied in the clinical care. This failure may be a result of the complex nature of seizures, which prohibits determination of a unique pattern for pre-seizure state in different seizures, even in one single patient. One the other hand, it is possible that the aforementioned features are not computed accurately, i.e., the features may be approximated in such a way that prevents the prediction to be accomplished. In the absence of a ground truth for the real signals, performances of the features’ computation methods are under question. In this paper, at first, a multi-channel depth-EEG model is explained, and then via a parameter identification procedure, it is shown that this model is able to reproduce real

depth-EEG signals. So, according to the synthetic signals produced by this spontaneous seizure generation model [40,41], we showed that the conventional methods lead to inaccurate estimations of features. From this point of view, it is not of main importance which seizure generation scenario underlies the synthetic signals [7], or what kind of seizure is simulated by them. The important thing is that some depth-EEG-like signals are at hand, by which it is possible to show that conventional methods can not compute their features accurately that may be the first reason of preventing seizures to be predicted. Then, it is shown that thanks to the proposed algorithm that is based on parameter identification procedure, the features of synthetic multi-channel signals are obtained accurately. Although this good result is due to the exact matching of the synthetic signals and the assumed model; but, because of the ability of model to produce real depth-EEG signals, the proposed algorithm can be generalized to real depth-EEG signals. Nevertheless, it is noteworthy that more feasible the assumed model, more accurate the features. The conventional synchronization measures have been taken into consideration in [41], in which it has been shown that in contrast to the conventional synchronization measures, the connectivity values identified based on the realistic depth-EEG model, would effectively quantify the synchronization between different areas of brain. But it is not discussed there [41] if the proposed algorithm still works in conditions in which the connectivity values change during time or not. Here, the concept of time varying connectivity values is covered. Thus, the structure of this paper is as follows. At first, the model generating depth-EEG signals is illustrated, and similarity if its output signals to the real depth-EEG signals is validated. Then the conventional algorithms of feature extraction, as well as the proposed one based on knowledge about the signal generation mechanism, are described. These algorithms are applied to the synthetic signals, and it is shown that the proposed algorithm can extract features of synthetic signals more accurate than the conventional techniques.

2.

Materials

2.1.

The seizure generator model

Although all factors playing role in ictogenesis procedure are not exactly known, unbalanced excitation and inhibition activity of populations of neurons are supposed to be the cause of seizure initiation [43]. In Shayegh et al. [40] neuronal excitation and inhibition variation are assumed to be stochastic processes, i.e., it is supposed that both excitatory and inhibitory gains of neural synapses change randomly. Thereupon, in this paper a two level stochastic model is assumed for generating depth-EEG signals in which the first layer simulates the stochastic variation of excitation and inhibition parameters. In the second layer there is a seizure genesis neural mass model whose parameters are excitation and inhibition properties of populations of neurons. This two-layer model is well described in [42] and is repeated in the Appendix A, and here we suffice to a short review of it. In the first layer of model four physiological parameters, including excitation gain, excitation time constant, fast

325

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

2.2.

Synthetic signals

Although number of channels does not influence our analysis technique, in this paper we considered five-channel signals. To model five nearby areas of hippocampus among which two are onset areas and others are drawn into seizure due to the effect of coupling, intrinsic parameters of onset areas (channels 1 and 4) are synthesized by running CHMM [42], and those of non-onset areas (channels 2, 3 and 5) are synthesized by random-walking. To have long-term EEG signals, for each channel about 100 h of these parameter sequences (A, B, G, ) were generated. A short (about 100 s) segment of these physiological sequences is represented in Fig. 1.

A

20 10

B

0 50 30 10 50 G

inhibition gain and slow inhibition gain, are considered. For the onset areas these random processes should behave such that the value of physiological parameters values can change easily. According to some physiological laws, Shayegh et al. [40] modeled these random processes by 4-dimensional Markov chain, but after being able to extract these physiological parameters via identifying the parameters of a neural mass model [41], Hidden Markov Model (HMM) was preferred to illustrate the excitation and inhibition behavior. Thus, the onset areas’ physiological parameter sequences are produced by running a 5-state continuous Hidden Markov Model. Since the parameter sequences are real-valued time series, a continuous hidden Markov model (CHMM) is considered in which the state-dependent density function is modeled by the mixture of Gaussian functions. But, to model the behavior of non-onset areas’ parameters that is fluctuation about some standard values, random walk process is more attractive. In the second layer, each area of the brain is modeled by a local computational model whose the underlying parameters are excitation gain, excitation time constant, fast inhibition gain and slow inhibition gain and is able to reproduce normal and epileptic activities [46]. The variables of this model (excitatory gain, fast and slow inhibitory gains and pyramidal excitatory constant time) are called A, B, G, and  respectively that are gathered in a parameter vector  = [A, B, G, ]. The equations of the model in the state space form are presented in (A.2). In this model [46] the effects of the sub-cortical and other cortical areas are simulated by a white Gaussian noise at the model input: p(t). This model is shown in Fig. A1. By connecting some of these local (single-channel) models together a multi-channel depth-EEG model is made. In order to extend the single-channel model to a multichannel one, each area of the hippocampal cortex is separately modeled by the single-channel model with the parameter  j whose input is no longer a pure white noise, but the combination of other channels’ output signals is added to it [8]. In order to simplicity, delays are neglected, i.e., the effects of each pair of areas are assumed to be instantaneous. This assumption is acceptable for multi-channel signals recorded from nearby areas. A schematic display of the multi-channel model is drawn at Fig. A2. The scales are denoted by coupling coefficients between each pair of EEG signals that are displayed by cij in Fig. A2. The parameter cij is the weight by which jth channel’s model output influences the input of the ith channel.

0 100 50 0

10

20

30

40

50 60 Time (sec)

70

80

90

100

Fig. 1 – A part of physiological parameters simulated via the CHMM for the onset areas and through the random walk process for non-onset areas.

Synthetic depth-EEG signals could be generated via driving the multi-channel depth-EEG model described in Section 2.1 by these intrinsic parameter sequences. But, to produce multi-channel depth-EEG signals, the coupling values must also properly be modeled. Here, for the sake of simplicity, all coupling values, except the coupling coefficients of the first onset area toward other four areas, are forced to be constant. / 1 are constant, and those In other words, the values of cij , j = of ci1 , i = 2,. . .,5 are allowed to change. The assumed scenario of this variation is shown in Fig. 2. Eventually, the intrinsic random parameter sequences and simulated coupling coefficients are used to drive the multichannel model in such a way that five-channel depth-EEG signals with 100 h in length would be at hand. The long-term dynamics of these synthesized depth-EEG signals are similar to what is observed in real [42]: very similar to what actually occurs on the depth-EEG signals of a patient in different physical and mental activities, as well as his/her seizure times, the background activity of the resultant signals is seizure free (i.e., interictal signals), and stochastically different activities (especially seizure-like activities) appear on these background EEG signals. A short segment of multi-channel signals produced by the seizure generator are shown in Fig. 3. Appearance of seizure-like signals on the onset channels that simulates the electrical seizure onset is shown in Fig. 3. It is assumed that electrical onset of seizures lead to a clinical manifestation only when seizure-like signals that spontaneously appeare on the onset channels (1 or 4) can propagate to other channels. As it is shown in Fig. 3, the seizure propagation occurs either when seizure appear on both onset channels simultaneously (rectangles with dot-dash line), or when only one channel initiate a seizure electrically, and at the same instant the areas are coupled together more strongly, i.e., in our case, when a synchronous increment of ci1 , i = 2,. . .,5 values

Fig. 2 – A scenario for the variation of coupling values ci1 , i = 2,. . .,5.

326

Ch1

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

5 0 -5

Ch2

1 0

Ch4

Ch3

-1 .5 0 -.5 5 0 -5

Coupling Value Ch 5

2 0 -2 10 5 0

0

10

20

30

40

50

60

70

80

90

100

Time(sec)

Fig. 3 – A segment of multi-channel signals produced using physiological parameters shown in Fig. 1. The scenario of the variation of coupling sequences ci1 is assumed to be the same for all values of i = 2,. . .5, and is shown at the last row of figure. Seizure-like signals can be appeared on channels one or four spontaneously, but seizures would propagate to other areas either when seizure appear on both onset channels (one and four) simultaneously (rectangles with dot-dash line), or when seizure appear on channel one while simultaneously the value of ci1 , i = 2,. . .,5 increase (rectangles with solid line). Otherwise, clinical manifestation of seizures would not be observed. In addition, there are some conditions in which the values ci1 , i = 2,. . .,5 increase but there is no seizure (rectangles with dash line). In the synthetic depth-EEG signals seizure occurrence is not a direct effect of synchronization increment, but conditionally is related to it.

and spontaneous seizure appearance on channel one is happened (rectangles with solid line); otherwise seizures would not be observed clinically. It is noteworthy that there are some situations in which the values ci1 , i = 2,. . .,5 increase, but no seizure is appeared (rectangles with dash line). Therefore, the proposed model of depth-EEG signals can simulate the states of epilepsy with two separate foci, in which seizures do not occur only due to an increase of neural synchronization, but is conditionally related to it. Accordingly, throughout the 100 h of simulated depth-EEG signals, 24 seizures are occurred. Once again, it should be mentioned that we do not use the synthetic EEG signals to evaluate how accurately predictable are seizures. Also, we do not aim to discuss whether the conventional features are able to predict seizures or not. Our claim is just that conventional methods of quantifying features are not reliable. In this manner the advantage of synthetic signals is that we know the generative equation set, as well as the parameter sequences underlying each channel of these signals, and also the coupling coefficients between channels. In other words, the phase space of the system generating the synthetic signals is exactly determined by the simulated parameter sequences, in such a way that true values of the phase space relevant features (e.g., correlation dimension and largest Lyapunov exponent) are known: to compute these features, Procaccia and Grassberger technique [12] can be applied to the true phase space. Thus, we can use these signals in order to evaluate the accuracy of these features’ quantifying algorithms. Also, synchronization quantification algorithms can be assessed according to the known coupling coefficients.

2.3.

Real signals

In this work the real signals under investigation are depthEEG signals of three temporal lobe epileptic patients from the

FSPEEG database ([39]; http://epilepsy.uni-freiburg.de). In the Freiburg database, the inter-ictal and ictal activities are available as separate records. Ictal activity files contain the signals from at least 50 min before seizure onset up to approximately 5 min after seizure termination. The onset and termination times of each seizure are also indicated in the database. Sampling frequency of the signals is 256 Hz. The number of recorded depth-EEG signal channels varies for different patients, but there are three patients for whom five depth-EEG signal channels are available. Five channels including two in-focus and three out-of-focus channels are considered.

2.4.

Validation of synthetic signals

To be able to use synthetic signals as the surrogate of real depth-EEG signals, in order to evaluate seizure predictor features, the similarity of synthetic and real signals must be assured. To confirm the ability of the model to simulate the realistic depth-EEG signals, in this paper, we have sought for the intrinsic physiological parameters of each channel, and the coupling coefficients corresponding to real depth-EEG signals, i.e., the parameters by which the model can reproduce the real signals. So, our plan is to do a parameter identification procedure, such that the synthetic signals corresponding to real depthEEG signals are the output of multi-channel model driven by the identified parameters.

2.4.1.

Parameter identification procedure

In order to identify the parameters underlying a segment of signal a suitable cost function is introduced and minimized. As it is demonstrated in [41,42], the single-channel model described by Eq. (A2) is a one-to-one function, such that assuming each channel of the observed signal as the output

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

of a single-channel model, the relevant input of each channel at each parameter vector  could be uniquely obtained. In the case of the multi-channel model, the optimum parameter values of each block must hold the best matching of that block’s input signal to a linear combination of the other areas’ signals. Thus, the cost function of the jth block is defined as the following:

Jj (j , Cj ) =

N  k=1

(x

k

(O1:N ; ) − obs,j j

M 

2

cji × Okobs,i )

(1)

i=1,i = / j

where xk (O1:N ;  ) is the estimated input signal related to the obs,j j , for parameter  j . N observed signal at the jth channel, O1:N obs,j is the number of signal samples, cij is the coupling value of channel i to channel j, and Cj is a vector containing cij ’s for all channels. The superscript k denotes for the kth sample of signal. Thus, assuming the observed signals as the outputs of the known multi-channel model, the connectivity values and intrinsic parameters of each channel can be estimated by minimizing the function (1). Now, it is possible to reproduce the real depth-EEG signals according to the identified parameters. But, to be able to drive the multi-channel depth-EEG model by these identified parameters some white Gaussian noise are required. It is evident that the value of cost function is not zero at the optimum parameters, i.e., there is always a residual term. Therefore, to obtain some synthetic signals very similar to the real ones, this white noise signal is retrieved from the residue term. At each channel, the whitened residue is taken as the white noise.

2.4.2. Comparison of real depth-EEG signals and their corresponding synthetic ones The parameters underlying the real depth-EEG signals were obtained through sequential solving the optimization problem for 1-s windows of signals. Driving the model by the resultant parameter sequences and the white noise obtained by whitening the residue term of each window gives us the synthetic signals corresponding to real depth-EEG signals. Two exemplar segments of real depth-EEG signals and their corresponding synthetic signals obtained according to the abovementioned procedure are shown in Fig. 4. Similarity of the signals of Fig. 4 is a small evidence of the fact that our seizure generator model is able to synthesize any real depth-EEG signal. In fact, although coupling is not wellmodeled in our proposed multi-channel model, it does not lose generality. This result can be generalized to any real signal. The visual similarity seems to be sufficient to convince anyone about validity of the model in simulating depth-EEG signals. But, to be more precious, the cross spectral difference between real and synthetic signals has be taken as a criterion. For all available real signals, the small average value of cross spectral difference about 0.03 confirms the validity of model. Hereupon, we safely can work on synthetic signals, while we are confident of the possibility of generalization the synthetic signals’ results to real depth-EEG signals.

3.

327

Methods

As it is mentioned in the Section 1, there are different reports about the efficiency of different features in predicting seizures. In other words, still there are many ambiguities about predictability of seizures based on EEG signals using these features [11]. It is not evident whether these features are predictors of seizures, and why there are both optimistic reports and pessimistic ones. From one point of view, one can think that these contradictory results may be the result of deficiency of feature estimation techniques, rather than the nature of these features. In other words, it is thought that all of the estimation approaches result in only some approximated values of the features, and not their precise values, in such a way that different approximation techniques could lead to different reports about the effectiveness of features as seizure precursory. To decide about the performances of these features in predicting seizures, there is no alternative except accurate computation of these features. Only in this way the results will be user-independent and the comparison of features’ performances would be fair. In this paper, we investigate about the accuracy of conventional methods in estimating features by using the realistic long-term simulated depth-EEG signals. It worths mentioning again that the advantage of synthetic signals is that we know their generative equation set, as well as the parameter sequences underlying each channel of these signals, and also the coupling coefficients between channels. In other words, the phase space of system generating the synthetic signals is attainable, such that true values of the phase space relevant features (e.g., correlation dimension and largest Lyapunov exponent) are known, i.e., to compute these features, Procaccia and Grassberger technique [12] can be applied to the original phase space without any requirement to reconstruct it. Also, considering the true coupling coefficients as original measures of effective connectivity, synchronization quantification algorithms can be properly evaluated. Thus, via these synthetic signals, evaluation of the feature extraction algorithms is feasible. Signal of the onset channel is used for evaluation of uni-variate features, and multi-channel signals are used for assessment of the bi-variate features. The conventional methods of estimating these features and the proposed one are illustrated in this section; then, in the next section, based on these ground truth values about the long-term simulated depth-EEG signals, the feature approximation methods would be assessed. As an example of the synthetic seizure-related depth-EEG signals used in the evaluation study, a short segment of the first onset channel signal containing a seizure, as well as a short segment of the five-channel signals are shown in Fig. 5.

3.1. Conventional methods for the estimation of features 3.1.1.

Correlation dimension

Correlation dimension is defined for a dynamic system in its phase space [19]. In this space the sequence of state points make a route that remains in a limited region or diverges, depending on the system. Usually points are attracted to a

328

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

Ch5

Ch4

Ch3

Ch2

Ch1

(a)

0

1

2

3

4

5

6

7

Time (sec)

Ch5

Ch4

Ch3

Ch2

Ch1

(b)

0

1

2

3

4

Time (sec)

5

6

7

Fig. 4 – Two exemplar segments of real depth-EEG signals (solid line) and their corresponding synthetic signals (dashed line) obtained according to the abovementioned parameter identification procedure. subset of space called attractor that can be a fixed point, a limit cycle or even a strange attractor with fractal structure. Correlation dimension corresponds to dimension of this attractor. The correlation dimension is equal to “0” for a fixed point and “1” for a limit cycle [29]. The result of [12] technique to compute correlation dimension from the state space of the system, is sufficiently accurate, even using limited (even noisy) points of the phase space. In this technique correlation dimension is obtained via the following equation:

  2 (ε − xi − xj ) N(N − 1) N

C(ε, N) =

N

i=1 J=i+1

(2)

∂ ln C(ε, N) d(ε, N) = ∂ ln ε where xi and xj are state vectors at discrete times i and j, respectively.  is the Heaviside function, C(∝ , N) is the

correlation sum such that its derivative d(␧, N) is the correlation dimension. However, for an observed signal, without explicit knowledge about its descriptor equations, still Grassberger and Procaccia technique can be used to approximate the correlation dimension value. But, this algorithm should be applied on reconstructed phase space, whose reconstruction is done by using the Takens method [45], i.e., time-delayed embedding approach. In Taken’s method the state vectors are constructed from several samples of the signal. In this order two parameters must be set: number of states (dimension), and delay between samples (delay). The closer is the reconstructed space to the true phase space, more accurate the estimation of correlation dimension will be. Thus, selection of the optimal values of (delay and dimension) these parameters are crucial. Theiler [44] proposed a slight modification of the standard algorithm in which the estimate of dimension is improved. In different papers, the estimation of these values has been done in different ways [20,21,37].

(a)

400

Samplar Output of the Seizure Denerator

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

300

329

200 100 0 -100 -200 0

1

2

3

4

5

6

7

8

9

10

11

12

Time(sec)

Ch5

Ch4

Ch3

Ch2

Ch1

(b)

0

0.5

1

1.5 Time (sec)

2

2.5

3

Fig. 5 – (a) A part of the signal of the channel one that includes a simulated seizure initiation and termination. Signal of the channel one is used to evaluate the performance of nonlinear features estimations. (b) A part of the multi-channel, showing a propagated seizure. Multi-channel signals are used to evaluate the performance of computing the synchronization measures.

3.1.2.

Largest Lyapunov exponent

The rate of divergence of the trajectories in the phase space is measured by the largest Lyapunov exponent, which is positive for a chaotic attractor [19]. The nearby points of the chaotic phase space diverge unceasingly. In contrast, two nearby points will stay close together forever, if they are on or near a limit cycle attractor, for which the largest characteristic exponent is zero. Two such points in the vicinity of a fixed point attractor will invariably get closer, i.e., the largest characteristic exponent is negative. In other words, the instability of a dynamic system is quantified by the largest Lyapunov exponent. In order to calculate the Lyapunov exponent, every point set of the phase space (sn0 ) must be assessed to compute the rate of trajectory divergence. The Eq. (3) fulfills this procedure:

⎛

1 1  ı(n) = ln ⎝  U(sn0 ) N N

n0 =1



⎞   sn0 +n − sn0 ⎠

(3)

sn0 ∈ U(sn0 )

where sn0 and sn0 +n are the reconstructed points of the embedded phase space at times n0 and n0 + n, respectively. U(s  n0 ) is the set of neighbor points of the sn0 in phase space, and U(sn0 ) is the number of members of this set. Lyapunov exponents are the exponents by which ı(n) grow on each dimension of the

phase space. Largest Lyapunov exponent is a good measure of the complexity changes’ rate [47]. Again, without knowing the true states of the system, computing the Lyapunov exponent, just based on the observed signal, requires reconstruction of the phase space. So, suitable delay and dimension parameters must be applied; otherwise, value of largest Lyapunov exponent would not be reliable.

3.1.3.

Synchronization

Synchronization is defined as the appearance of certain relations between the phases and frequencies of some interacting objects [38]. These relations may be of the kind of detuning caused by coupling. Therefore, in order to measure the value of synchronization in a particular experiment, we must either have access to some parameters of the oscillators (or to the external force) that govern the frequency detuning, or be able to vary the coupling values. Actually, we have to look what happens to the frequencies and/or phase difference under variation of one of these parameters, i.e., an “active” experiment should be performed. When one cannot change the parameters of the systems and/or amount of couplings, and just observes the signals under free-running conditions, like the situation for brain signals, synchronization cannot be detected from the data

330

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

without additional assumptions. In fact, since synchronization is not a state, but a process of adjustment of phases and frequencies, its presence or absence cannot be established from a single observation. Nevertheless, different endeavors are done to measure synchronization of signals as interdependence between two (or more) signals. These techniques explicitly uses the assumption that the data are generated by several (at least two) interacting self-sustained systems. If the probability of having the same phase/frequency contents for two signals is too low per se, these measures may provide useful information about the interrelation of the systems underlying signals [38]. This is usually the case for non-stationary signals, but is not always true. Since we cannot check whether this similarity happens occasionally or it is a manifestation of interaction, we just obtain some approximate information about the interrelation of the systems in a passive manner. Phase coherence, as formulated in Eq. (4), is one of the passive synchronization approximators: 2

2 Rn,m = E[cos (ϕn,m )] + E[sin (ϕn,m )]







2

(4)



where ϕn,me  = nϕx (t) − mϕy (t) , and ϕx (t) and ϕy (t) are the instantaneous phase of the signals x and y that can be obtained either using the analytic approach according to the Hilbert transform, or using a definition based on the Wavelet transform [36]; and E(·) denotes mathematical expectation.

state space underlying the system of each channel of observed signal could be reconstructed, in such a way that the estimated features would be more accurate than algorithms in which the phase space is reconstructed using Taken’s method. In other words, the parameter identification procedure described in Section 2.4.1 is the core of the proposed feature extraction algorithm. Thus, at first the parameter identification procedure should be applied to the consequent windows of depth-EEG signal. In the parameter identification procedure, the parameters are supposed to have constant values during each window of the signal, therefore, short segments of the signal must be considered. However, the segments must be large enough to reveal the frequency contents of the signal. So, the overlapped segments of the signals are taken into account. Accordingly, parameters are identified from overlapped windows of the signal with 1 s length (771 samples) and 66% overlap. The estimated parameter sequences of the signal shown in Fig. 5(a) are given in Fig. 6. It is evident from Fig. 6 that although the small fluctuation of the parameter sequence can not properly be estimated, difference between true and estimated intrinsic parameters is negligible. According to the identified parameters, the model equation set of each channel implies a state space very close to the true state space. At the second step it is proposed to obtain the uni-variate features of each channel from this state-space via the Procaccia and Grassberger technique.

4.

3.2. Proposing a better approach to estimate the features Given some depth-EEG signals, it is expected that knowledge about the system underlying them would help us to estimate the values of different features more accurately. According to the ability of the multi-channel model to produce realistic depth-EEG signals, the proposed feature estimation approach is based on estimating the parameters of multichannel model, i.e., intrinsic parameters of each channel, as well as the sequences of coupling coefficients. Then, according to these estimated parameters and the equation set (A2), the

Results

The above-described algorithms (conventional and the proposed one) are applied to approximate the features of the synthetic EEG signals of 100 h in length. As it is implied in Section 3.2, features are extracted from overlapped windows of the signal with 1 s length (771 samples) and 66% overlap. For the signals of the first onset area the correlation dimension and largest Lyapunov exponents are calculated, while the synchronization is considered for the five channels altogether. To be more conscious, the obtained values of features are compared with their true values in the MSE manner according

A

100

True values Estimated

50

B

0 100 50

G

0 100 50 0 100 50 0 0

1

2

3

4

5

6 7 Time (sec)

8

9

10

11

12

Fig. 6 – The true and identified parameters of the synthetic signal shown in Fig. 5(a). The solid diagrams show the true values and dotted ones show the estimated values.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

to (5). It worths mentioning again that for synthetic signals true values of the features are at hand.



n

 ¯ 2 MSEj =  (fij − fˆ ij )

(5)

i=1

where n is the number of the signal windows under consideration, and i is the index of the of the under process window. j denotes the true value of feature vector, in which j is the feature number, and respectively correspond to the correlation dimension (j = 1), largest Lyapunov exponent (j = 2), and the coupling values (j = 3:22). Similarly fˆ is the estimated feature vector. To compensate the situations in which the true and estimated features do not have the same ranges, but are scaled versions of each other, the normalized values of all features are considered in the MSE calculation. For example, mean phase coherence values (in contrast to the true coupling coefficient values) are between “0” and “1”, thus normalization of coupling coefficients is required to make their subtraction possible. The over-bar implies the normalized values of features. Normalization is done according to the maximum values of each feature along time.

4.1.

Calculating the features by conventional methods

The conventional methods described in Section 3.1 are used to estimate correlation dimension and largest Lyapunov exponent of signals of the onset channel. The choice of the embedding dimension is done based on a reliable method that corrects the undesired systematic effects emerged due to temporal oversampling, autocorrelation, and changing time lag [21]. For determining the time delay, values of the mutual information between the signal and its delayed versions are calculated; then the index of first minimum of these values is chosen as the time delay. The mean phase coherence described in Section 3.1 is used to quantify the amount of synchronization between different channels of the multi-channel signal. The MSE values for all 100 h of synthetic signal, either during seizure or non-seizure states, are reported in Table 1 for the conventional methods. In this section fˆi1 and fˆi2 are correlation dimension and the largest Lyapunov exponent respectively, and fˆij , j = 3 : 22 are the mean coherence values, e.g., fˆi3 is equivalent to the mean phase coherence between the signal of the first channel and that of the second channel (i.e., R21 ), for the ith window. MSE values of correlation dimension and the largest Lyapunov exponent, and 20 synchronization (connectivity) values are given in the rows of Table 1. Accordingly, large values of MSEs show the inability of conventional feature estimation algorithms. To be more intuitive, the results obtained for the signals of Fig. 5 are shown in Figs. 7 and 8. The dashed diagrams in Figs. 7(a) and (b) and 8 stand for the estimated values of correlation dimension, the largest Lyapunov exponent and amount of synchronization between signals, respectively. In these figures, the solid diagrams show the ground truth of these features. As it is seen from Fig. 7(a) and (b), the estimated values of nonlinear features are far from their true values.

331

Table 1 – The mean square difference between true values of the features and features estimated by conventional methods, for 100 h of signal. The first two rows of the table are assigned to the correlation dimension, and the largest Lyapunov exponent, respectively. MSE values of 20 connectivity values are given in rows 3–22. It is noteworthy that both true and estimated features are normalized according to their maximum values over time. Features

MSE

Correlation dimension Largest Lyapunov exponent c21 (fi3 = c21 , fˆi3 = R21 ) c31 (fi4 = c31 , fˆi4 = R31 ) c41 (fi5 = c41 , fˆi5 = R41 ) c51 (fi6 = c51 , fˆi6 = R51 ) c12 (fi7 = c12 , fˆi7 = R12 ) c32 (fi8 = c32 , fˆi8 = R32 ) c42 (fi9 = c42 , fˆi9 = R42 ) c52 (fi,10 = c52 , fˆi,10 = R52 ) c13 (fi,11 = c13 , fˆi,11 = R13 ) c23 (fi,12 = c23 , fˆi,12 = R23 ) c43 (fi,13 = c43 , fˆi,13 = R43 ) c53 (fi,14 = c53 , fˆi,14 = R53 ) c14 (fi,15 = c14 , fˆi,15 = R14 ) c24 (fi,16 = c24 , fˆi,16 = R24 ) c34 (fi,17 = c34 , fˆi,17 = R34 ) c54 (fi,18 = c54 , fˆi,18 = R54 ) c15 (fi,19 = c15 , fˆi,19 = R15 ) c25 (fi,20 = c25 , fˆi,20 = R25 ) c35 (fi,21 = c35 , fˆi,21 = R35 ) c45 (fi,22 = c45 , fˆi,22 = R45 )

0.752 0.439 0.291 0.314 0.346 0.305 0.412 0.375 0.405 0.236 0.396 0.332 0.325 0.206 0.572 0.436 0.519 0.550 0.312 0.271 0.212 0.532

Even for some cases dynamics of the features are not followed by the estimated values properly. In Fig. 8(a) the true synchronization values of ci1 , i = 2,. . .,5 and ci4 , i = 1,2,3,5 are displayed. The estimated mean phase 2 ,i = coherence values corresponding to these channels, i.e., Ri1 2 , i = 1, 2, 3, 5, are shown in Fig. 8(b) and (c). The 2, ..., 5 and Ri4 time course of Fig. 8(a) is not appropriately preserved in the 2 , i = 2, ..., 5. Also, as mentioned in phase coherence values Ri1 Section 2.2, ci4 , i = 1, 2, 3, 5 is supposed to be constant over 2 , i = 1, 2, 3, 5 values show some fluctuation. time, but the Ri4 As reported in Table 1, for other pairs of channels, the mismatch between the true and approximated synchronization values is also present. In other words, the mean phase coherence is not a reliable measure for quantifying the amount of synchronization.

4.2. Calculating the features by the proposed algorithm As mentioned in Section 3.2, by optimizing the Eq. (1), sequentially for the short segments of signal with 1-s length, the ˆ ˆ ) of the signals of every channel, ˆ B, ˆ G, intrinsic parameters (A, ˆ j would be obtained. as well as the connectivity values C It is evident from Fig. 6 that these estimated parameters are not exactly the original ones, due to the errors of the identification algorithm, as well as the assumption of constant parameters during short segments of signal. By replacing the

332

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

Fig. 7 – Values of (a) correlation dimension and (b) largest Lyapunov exponent for the signal shown in Fig. 5(a). Features are obtained for subsequent segments of the signal with 1-s duration and 66% overlap. Result of the conventional embedding approach (dash line) and the proposed approach (dotted line) are compared to the true (original) values of these features (solid line).

identified physiological parameter sequences of the 100 h signals in the equation set (A2) the state space of the system underlying each channel of signals have been constructed. Then, from this state space, the correlation dimension and largest Lyapunov exponents are computed using Grassberger method [12]. Clearly, these reconstructed state spaces are very close to the true ones, but not exactly the same, because of the errors of the identified parameters. Due to this fact the features are not exactly the same as their true values. For all 100 h of signal, the MSE values between true and estimated features, either during seizure or non-seizure states, are reported in Table 2 for the proposed methods. Two first rows of the tables are assigned to the correlation dimension and the largest Lyapunov exponent, respectively. In the proposed method, for the multi-channel signals the identified coupling values can be considered as a measure of synchronization. In other words, fˆij , j = 3 : 22 are the pair-wise coupling coefficients, e.g., fˆi3 means the coupling coefficient between the signal of the first channel and that of the second channel (i.e., cˆ 21 ). The MSE values between 20 estimated connectivity values and their true values are given in Table 2, in rows 3–22. Although the MSE values reported in Table 2 are not equal to zero, due to the non-avoidable errors of the identification algorithm caused by assuming constant values during short intervals of signal, but these values are very smaller than the MSE values of Table 1 (estimated values using conventional methods). The obtained results shown in Table 2 imply that by using “Evolution” equations, which embody the system’s dynamics, the obtained measures outperform statistical techniques of detection of synchrony. In addition to the MSE values, the uni-variate features of the onset channel signals of Fig. 5(a), and the amount of synchronization between five-channel signals of Fig. 5(b) obtained by the proposed technique are compared to their true values in Figs. 7 and 8. Dotted diagrams of Fig. 7(a) and (b) displays the correlation dimension and largest Lyapunov exponent obtained by the proposed algorithm, and the dashed diagrams stand for the conventional methods. The superiority of the proposed method compared with the conventional method is evident

from these figures. Using the proposed algorithm, both correlation dimension and largest Lyapunov exponent values are closer to their true values than the results of conventional approximation methods. The coupling values from channel 1, cˆ i1 , i = 2, ..., 5, and channel 4, cˆ i4 , i = 1, 2, 3, 5 are shown in Fig. 8(d) and (e). According to these figures, it is seen that the identified connectivity values are very close to their true values.

Table 2 – The mean square difference of true values of the features and estimated features by the proposed method, for 100 h of signal. The first two rows of the table are assigned to the correlation dimension and the largest Lyapunov exponent, respectively. MSE values of 20 connectivity values are given in rows 3–22. It is noteworthy that both true and estimated features are normalized according to their maximum values over time. Features

MSE

Correlation dimension Largest Lyapunov exponent c21 (fi3 = c21 , fˆi3 = cˆ 21 ) c31 (fi4 = c31 , fˆi4 = cˆ 31 ) c41 (fi5 = c41 , fˆi5 = cˆ 41 ) c51 (fi6 = c51 , fˆi6 = cˆ 51 ) c12 (fi7 = c12 , fˆi7 = cˆ 12 ) c32 (fi8 = c32 , fˆi8 = cˆ 32 ) c42 (fi9 = c42 , fˆi9 = cˆ 42 ) c52 (fi,10 = c52 , fˆi,10 = cˆ 52 ) c13 (fi,11 = c13 , fˆi,11 = cˆ 13 ) c23 (fi,12 = c23 , fˆi,12 = cˆ 23 ) c43 (fi,13 = c43 , fˆi,13 = cˆ 43 ) c53 (fi,14 = c53 , fˆi,14 = cˆ 53 ) c14 (fi,15 = c14 , fˆi,15 = cˆ 14 ) c24 (fi,16 = c24 , fˆi,16 = cˆ 24 ) c34 (fi,17 = c34 , fˆi,17 = cˆ 34 ) c54 (fi,18 = c54 , fˆi,18 = cˆ 54 ) c15 (fi,19 = c15 , fˆi,19 = cˆ 15 ) c25 (fi,20 = c25 , fˆi,20 = cˆ 25 ) c35 (fi,21 = c35 , fˆi,21 = cˆ 35 ) c45 (fi,22 = c45 , fˆi,22 = cˆ 45 )

0.22 0.29 0.15 0.19 0.19 0.13 0.18 0.26 0.07 0.20 0.19 0.15 0.26 0.08 0.06 0.14 0.13 0.151 0.17 0.26 0.08 0.30

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

333

Fig. 8 – (a) The sequences of the coupling coefficients ci1 , i = 2,. . .,5 of the multi-channel signal shown in Fig. 5(b). (b) The 2 , i = 2, ..., 5, (c) The estimated phase coherence values of channel 4, estimated phase coherence values of channel 1, i.e., Ri1 2 , i = 1, 2, 3, 5, (d) The estimated values of synchronization measure c ˆ i1 , i = 2, ..., 5 by proposed approach, and (e) i.e., Ri4 cˆ i4 , i = 1, 2, 3, 5. Synchronization values estimated by the proposed approach are very close to the true coupling coefficients.

5.

Discussion and conclusion

In this paper, a model-based evaluation of the algorithms that estimate of some precursory features for epileptic seizure is presented. In fact, in order to compare different approaches used to predict seizure, a standard dataset is required. Although the Freiburg database1 [39] has been used by many researchers recently, since true values of features of these signals are unknown, the performance of feature extraction techniques could not be evaluated. This problem is eliminated

1 https://epilepsy.uni-freiburg.de/freiburg-seizure-predictionproject/eeg-database.

for synthetic signals, for which true values of features are known. Three features, i.e., correlation dimension, the largest Lyapunov exponent, and synchronization measures are considered in this paper. In this paper, via a physiological plausible multi-channel seizure generator that its ability to reproduce real depthEEG signal is shown, some synthetic multi-channel depth-EEG signals are generated. In other words, these synthetic depthEEG signals are similar to the real ones. According to these synthetic signals it is illustrated that when the feature extraction techniques just work based on the signals, without any attention to their origins, they can not successfully estimate features. In other words, results represented in Table 1 and Figs. 7 and 8 show that these techniques provide a rough approximation of true values of the features. For example,

334

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

phase coherence may show some decrease or increase (Fig. 8(c)) when the coupling between brain areas does not change at all. In other words, it falsely detects synchronization increase when signals become more similar to each other, even when the coupling increase does not occur. Finally, as shown in Table 2, by considering the origin of depth-EEG signals, i.e., using the realistic multi-channel model as the structure of the signal generator system, the accuracy of the estimated feature would be improved significantly.

Appendix A The model In this paper, a two level stochastic model is assumed for depth-EEG signals in which the first layer simulates the stochastic variation of excitation and inhibition parameters and in the second layer there is a computational model of seizure genesis whose parameters are excitation and inhibition properties of populations of neurons [42].

A.1. The first layer In the first layer of model four physiological parameters, including excitation gain, excitation time constant, fast inhibition gain and slow inhibition gain are considered. The onset

areas’ physiological parameter sequences are produced by running a 5-state Hidden Markov Model. HMM is a powerful statistical tool for modeling an observable sequence that can be generated by an underlying probabilistic process. In fact we assume that intrinsic physiological parameters of an onset channel are some probabilistic function of the states that are not observable. Since the parameter sequences are real-valued time series, a continuous hidden Markov model (CHMM) is considered in which the state-dependent density function is modeled by the mixture of Gaussian functions. For an N-state CHMM, in which all states have K mixture components, the characteristic vector is f = { N , PN×N , CN×K , N×K } , where is the initial-state probability vector, P is the transition matrix, C is the mixture-coefficient matrix, and  is composed of the Gaussian parameters ≥ik = {. . . ik , ≤ ik } for the kth mixture component of the ith state with ik and k as mean and covariance matrix, respectively. The value of CHMM parameters used to simulate the onset parameter sequences are given in Table A1. It is reasonable to assume that the parameter values of non-onset areas have small changes around their standard values. Thus, random walk process is suitable to model the non-onset areas’ parameter behavior.The sequence {Sn }∞ n=1 is called a random walk, if it was the summation of the form S0 + X1 + X2 + ... + Xn , where {Xk }∞ k=1 is a sequence of

Table A1 – Characteristic of the CHMM used to simulate the parameter sequences of the onset areas. CHMM has five states (N = 5), and the probability density function of each state contains three Gaussian function (K = 3).  is the initial-state probability vector, P is the transition matrix, C is the mixture-coefficient matrix, and  and indicates the mean and variance of the Gaussian functions at every state. CHMM parameter

P

C

A , ˙ A

B , ˙ B

G , ˙ G

 , ˙ 

0 ⎡

1

0

0

0



Value

⎤

0.76 0.0654 0.0866 0.0024 0.0857 0.8338 0.0030 0.0267 0.0222 ⎥ ⎢ 0.1143 −4 ⎢ 7.5278 × 10 2.003 × 10−56 0.8126 0.0223 0.1644 ⎥ ⎣ ⎦ 0.0706 0.0110 0.8796 0.0245 0.0143 0.0845 0.1084 0.0560 0.7355 ⎡ 0.0156 −4 ⎤ 9.2383 × 10 0.0019 0.9971 −4 0.1249 0.8746 ⎥ ⎢ 5.6427 × 10−4 ⎢ 6.8483 × 10 6.8483 × 10−4 0.9986 ⎥ ⎣ ⎦ 0.4409 0.3642 0.1949 −4 −4 ⎡ 2.5619 × 10 4.0479 × 10⎤ ⎡ 0.9993−6 ⎤ 96.8909 93.5964 1.0553 1 × 10 3.0542 0.0680 −6 ⎢ 92.9598 2.6454 1.8245 ⎥ ⎢ 1 × 10−6 2.3103−6 0.7892 ⎥ ⎢ 3.0919 29.1822 1.2400 ⎥,⎢ 1 × 10 1 × 10 0.2133 ⎥ ⎣ ⎦⎣ ⎦ 21.4415 1.6916 85.4345 31.8043 3.4993 8.7646 −6 −6 ⎡ 65.1664 65.1664 1.5519 ⎤ ⎡1 × 10 −6 1 × 10 0.4531 ⎤ 71.3227 77.8916 8.2208 1 × 10 0.0459 4.4592 −6 15.1439 ⎥ ⎢ 98.3448 28.3879 24.4952 ⎥ ⎢ 1 × 10−6 26.6670 ⎢ 74.2511 75.5335 2.0167 ⎥,⎢ 1 × 10 1 × 10−6 0.2133 ⎥ ⎣ ⎦⎣ ⎦ 82.6736 37.2641 6.3176 11.5213 22.8568 4.9197 −6 −6 7.0219 ⎡ 70.1789 70.1879 9.2160 ⎤ ⎡ 1 × 10−6 1 × 10 ⎤ 22.3531 0.4758 1 × 10 8.7196 6.5182 5.8417 −6 18.2984 ⎥ ⎢ 76.3469 45.6823 30.5070 ⎥ ⎢ 1 × 10−6 28.8496 ⎢ 8.4513 5.0384 1.2932 ⎥,⎢ 1 × 10 1 × 10−6 0.3603 ⎥ ⎣ ⎦⎣ ⎦ 47.0159 23.5231 5.8135 30.3430 22.6436 4.3407 −6 −6 2.1863 ⎡ 44.3004 44.3004 3.1558 ⎤ ⎡ 1 × 10−6 1 × 10 ⎤ 97.9264 10.3684 40.7030 1 × 10 3.5048 20.1739 −6 ⎢ 14.1291 2.1501 35.6002 ⎥ ⎢ 1 × 10−6 0.8971−6 18.0538 ⎥ ⎢ 21.8738 1.0744 28.9958 ⎥,⎢ 1 × 10 1 × 10 16.7564 ⎥ ⎣ ⎦⎣ ⎦ 25.2375 28.8457 53.2104 25.2375 28.8457 28.2104 −6 −6 1 × 10 10.3882 95.5205 95.5205 17.2711 1 × 10

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

335

Fig. A1 – Basic structure of EEG models including linear dynamic and nonlinear static transforms, respectively at the synapses and soma.

independent, identically distributed discrete random variables [13]. Since there are four physiological parameters in the modeling, the common range of the Xk ’s is R4 , i.e. the random vector Xk contains four variables x1 , x2 , x3 , x4 . It is assumed that the walk starts from the standard value of parameters at time n = 0, here S0 = [3.25,22,10,100] (Table A2). We consider the simplest non-trivial case of a random walk, namely the case where the common distribution function of the random variables Xk is given by:

fX (x) =

⎧  4 ⎪ ⎨ 1 x1 = ±1 & x2 = ±1 & x3 = ±1 & x4 = ±1 2

⎪ ⎩ 0 others

Fig. A2 – The local model, the relation of excitatory neurons (both pyramidal cells and interneurons), slow inhibitory interneurons and fast inhibitory interneurons (the gray rectangle) are shown; h(t), he (t), hg (t) are the impulse response of these systems, respectively, and A, B and G are their synaptic gains [46].

(A1)

A.2. The second layer In the second layer each area of the brain is modeled by a local computational model. By connecting some of these local models together a global (multi-channel) model is made.

A.2.1. Local (single channel) model n area of the hippocampus consisting of neuronal subpopulations, i.e. pyramidal neurons and interneurons, has been modeled in different manners. Two transforms, a linear dynamic and a nonlinear staticas shown in Fig. A1, respectively at the synapses and soma, are the essence of these models [6]. The sub-populations interact with each other through inhibitory or excitatory relations. The pyramidal

neurons are always modeled by an excitatory sub-population, but the interneurons have been modeled by different types of sub-populations. In the reality fast and slow inhibitory and also excitatory interneurons coexist [43]; however, models may consider just some of these sub-populations. By considering the three sub-populations of interneurons, a model can generate normal and epileptic activities, even the fast discharge activity [46]. The model described in [46] is a pure cortical model of four sub-populations: pyramidal cells, excitatory interneurons, fast and slow inhibitory interneurons. In this model, the effects of the sub-cortical and other cortical areas are simulated by a white Guassian noise at the model input p(t), as shown in Fig. 2. This model is thus a local model of hippocampal cortex, and it can represent only the depth-EEG signal of a single channel [46]. This model is shown in Fig. A2. The model is described mathematically by equation set (A2) in the state space. The physiological facts corre-

Table A2 – The interpretation and standard values of the local model parameters [46]. Parameter A B G  a b g C1 , C2 C3 , C4 C5 , C6 C7 v0 , e0 , r

Interpretation Average excitatory synaptic gain Average slow inhibitory synaptic gain Average fast inhibitory synaptic gain Dendritic average time constant in the feedback excitatory loop of pyramidal neurons Dendritic average time constant in the feedback excitatory loop of interneurons Dendritic average time constant in the slow feedback inhibitory loop Somatic average time constant in the fast feedback inhibitory loop Average number of synaptic contacts in the excitatory feedback loop Average number of synaptic contacts in the slow feedback inhibitory loop Average number of synaptic contacts in the fast feedback inhibitory loop Average number of synaptic contacts between slow and fast inhibitory interneurons Parameters of the nonlinear asymmetric sigmoid function

Standard value 3.25 (mV) 22 (mV) 10 (mV) 1/100 (s) 1/100 (s) 1/50 (s) 1/500 (s) C1 = 135, C2 = 108 C3 = C4 = 33.75 C5 = 40.5, C6 = 13.5 C7 = 108 v0 = 6 (mV), e0 = 2.5 (s−1 ), r = 0.56 (1/mV)

336

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

the simplicity, the delays of the information transmission between areas are neglected. This assumption is reasonably acceptable for adjacent electrodes. The scales are denoted by coupling coefficients between EEG signals, represented by cij in Fig. A3. The parameter cij is the weight by which jth area influences the input of the ith area.

references

Fig. A3 – The global model, the coupling between areas is modeled by adding the scaled outputs of different areas to the white noise at the input of one area. Each block is equivalent to a single-channel model shown in Fig. A2.

sponding to the model coefficients and their standard values are presented in Table A2. The parameters of this model are excitatory gain (A), fast and slow inhibitory gains (B and G), and pyramidal excitatory constant time () of the linear dynamic part of the subpopulations’ transfer function. These synaptic parameters are gathered in a parameter vector  = [A, B, G, ], whose behavior is modeled by the processes illustrated in Section A.1. y˙ 0 (t) = y5 (t) y˙ 5 (t) = AS[y1 (t) − y2 (t) − y3 (t)] − 2y5 (t) −  2 y0 (t) y˙ 1 (t) = y6 (t) y˙ 6 (t) = Aa{p(t) + C2 S[C1 y0 (t)]} − 2ay6 (t) − a2 y1 (t) y˙ 2 (t) = y7 (t) y˙ 7 (t) = BbC4 S[C3 y0 (t)] − 2by7 (t) − b2 y2 (t)

(A2)

y˙ 3 (t) = y8 (t) y˙ 8 = GgC7 S[C3 y0 (t) − C6 y4 (t)] − 2gy3 (t) − g2 y3 (t) y˙ 4 (t) = y9 (t) y˙ 9 = BbS[C3 y0 (t)] − 2by4 (t) − b2 y4 (t) O(t) = y1 (t) − y2 (t) − y3 (t)

A.2.2. Global (multi channel) model To generalize the single-channel model to a multi-channel one, each area of the hippocampus, is separately modeled as a single-channel model with parameter  j . The connection between these areas is modeled by adding delayed and scaled outputs of all other areas to the white noise at the input of each single-channel model [8]. Indeed, the input of the single-channel model is no longer a pure white noise, but the combination of EEG signals observed at other channels is added to it. In fact, the white noise signal corresponds to the effects induced by unrecorded areas. A schematic display of the multi-channel model is drawn at Fig. A3. For the sake of

[1] A. Aarabi, B. He, A rule-based seizure prediction method for focal neocortical epilepsy, Clinical Neurophysiology 123 (2012) 1111–1122. [2] R. Aschenbrenner-Scheibe, T. Maiwald, M. Winterhalder, H.U. Voss, J. Timmer, A. Schulze-Bonhage, How well can epileptic seizures be predicted? An evaluation of a nonlinear method, Brain 126 (2003) 2616–2626. [3] R.G. Andrzejak, F. Mormann, T. Kreuz, C. Rieke, A. Kraskov, C.E. Elger, K. Lehnertz, Testing the null hypothesis of the non-existence of the pre-seizure state, Physical Review E 67 (2003) 010901. [4] R.G. Andrzejak, F. Mormann, G. Widman, T. Kreuz, C.E. Elger, K. Lehnertz, Improved spatial characterization of the epileptic brain by focusing on nonlinearity, Epilepsy Research 69 (2006) 30–44. [5] M. Chavez, M. Le Van Quyen, M.V. Navarro, M.M. Baulac, J. Martinerie, Spatio-temporal dynamics prior to neocortical seizures: amplitude versus phase couplings, IEEE Transactions on Biomedical Engineering 50 (2003) 571–583. [6] F.H.L. da Silva, A. Hoeks, H. Smits, L.H. Zetterberg, Model of brain rhythmic activity: the alpha-rhythm of the thalamus, Kybernetik 15 (1974) 27–37. [7] L. da Silva, W. Blanes, S. Kalitzin, J. Parra, P. Suffczynski, D.N. Velis, Epilepsies as a dynamical disease of brain systems: basic models of the transition between normal and epileptic activity, Epilepsia 44 (2003) 72–83. [8] O. David, K.J. Friston, A neural mass model for MEG/EEG: coupling and neuronal dynamics, NeuroImage 20 (2003) 1743–1755. [9] C.E. Elger, K. Lehnertz, Seizure prediction by non-linear time series analysis of brain electrical activity, European Journal of Neuroscience 10 (1998) 786–789. [10] J. Engel, T.A. Pedley, Epilepsy: A Comprehensive Textbook, Lippincott Williams & Wilkins, 2007. [11] M.G. Frei, H.P. Zaveri, S. Arthurs, G.K. Bergey, C.C. Jouny, K. Lehnertz, J. Gotman, I. Osorio, T.I. Netoff, W.J. Freeman, J. Jefferys, G. Worrell, M. Le Van Quyen, S.J. Schiff, F. Mormann, Controversies in epilepsy; debates held during the fourth international workshop on seizure prediction, Epilepsy and Behavior 19 (2010) 4–16. [12] P. Grassberger, I. Procaccia, Estimation of the kolmogorov entropy from a chaotic signal, Physical Review A 28 (1983) 2591–2593. [13] C.M. Grinstead, J.L. Snell, Introduction to Probability, American Mathematical Society, 1997. [14] M.A.F. Harrison, M.G. Frei, I. Osorio, Accumulated energy revisited, Clinical Neurophysiology 116 (2005) 527–531. [15] L.D. Iasemidis, J.C. Sackellares, H.P. Zaveri, W.J. Williams, Phase space topography and the Lyapunov exponent of electrocorticograms in partial seizures, Brain Topography 2 (1990) 187–201. [16] L.D. Iasemidis, P. Pardalos, J.C. Sackellares, D.S. Shiau, Quadratic binary programming and dynamical system approach to determine the predictability of epileptic seizures, Journal of Combinatorial Optimization 5 (2001) 9–26.

c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 1 3 ( 2 0 1 4 ) 323–337

[17] L.D. Iasemidis, D.S. Shiau, P.M. Pardalos, W. Chaovalitwongse, K. Narayanan, A. Prasad, K. Tsakalis, P.R. Carney, J.C. Sackellares, Long-term prospective on-line real-time seizure prediction, Clinical Neurophysiology 116 (2005) 532–544. [18] L.D. Iasemidis, D.S. Shiau, W. Chaovalitwongse, J.C. Sackellares, P.M. Pardalos, J.C. Principe, et al., Adaptive epileptic seizure prediction system, IEEE Transactions on Biomedical Engineering 50 (2003) 616–627. [19] H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, Cambridge Univ. Press, Cambridge, UK, 1997. [20] M.B. Kennel, R. Brown, H.D.I. Abarbanel, Determining embedding dimension for phase-space reconstruction using a geometrical construction, Physical Review A 45 (1992) 3403–3411. [21] M.B. Kennel, H.D. Abarbanel, False neighbors and false strands: a reliable minimum embedding dimension algorithm, Physical Review E 66 (2002) 026209. [22] Y.C. Lai, M.A. Harrison, M.G. Frei, I. Osorio, Inability of Lyapunov exponents to predict epileptic seizures, Physical Review Letters 91 (2003) 068102. [23] Y.C. Lai, M.A. Harrison, M.G. Frei, I. Osorio, Controlled test for predictive power of Lyapunov exponents: their inability to predict epileptic seizures, Chaos 14 (2004) 630–642. [24] K. Lehnertz, C.E. Elger, Can epileptic seizures be predicted? Evidence from nonlinear time series analysis of brain electrical activity, Physical Review Letters 80 (1998) 5019–5022. [25] M. Le Van Quyen, J. Martinerie, V. Navarro, M. Baulac, F.J. Varela, Characterizing neurodynamic changes before seizures, Journal of Clinical Neurophysiology 18 (2001) 191–208. [26] M. Le Van Quyen, V. Navarro, J. Martinerie, M. Baulac, F.J. Varela, Toward a neurodynamical understanding of ictogenesis, Epilepsia 44 (2003) 30–43. [27] M. Le Van Quyen, J. Soss, V. Navarro, R. Robertson, M. Chavez, M. Baulac, J. Martinerie, Preictal state identification by synchronization changes in long-term intracranial EEG recordings, Clinical Neurophysiology 116 (2005) 559–568. [28] J. Martinerie, C. Adam, M. Le Van Quyen, M. Baulac, S. Clemenceau, B. Renault, F.J. Varela, Epileptic seizures can be anticipated by non-linear analysis, Nature Medicine 4 (1998) 1173–1176. [29] J.A. McCarthy, in: Goethe, Nietzsche, Grass (Eds.), Remapping reality: chaos and creativity in science and literature, Rodopi, Amsterdam & New York, 2006, p. 320. [30] P.E. McSharry, L.A. Smith, L. Tarassenko, Prediction of epileptic seizures: are nonlinear methods relevant? Nature Medicine 9 (2003) 241–242. [31] N. Moghim, D. Cornen, Evaluating bio-inspired approaches for advance prediction of epileptic seizures, in: Nature and Biologically Inspired Computing (NaBIC), 2011, pp. 540–545. [32] F. Mormann, K. Lehnertz, P. David, C.E. Elger, Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients, Physica D: Nonlinear Phenomena 144 (2000) 358–369.

337

[33] F. Mormann, T. Kreuz, R.G. Andrzejak, P. David, K. Lehnertz, C.E. Elger, Epileptic seizures are preceded by a decrease in synchronization, Epilepsy Research 53 (2003) 173–185. [34] F. Mormann, R.G. Andrzejak, T. Kreuz, C. Rieke, P. David, C. E Elger, Automated detection of a pre-seizure state based on a decrease in synchronization in intracranial EEG recordings from epilepsy patients, Physical Review E 67 (2003) 021912. [35] F. Mormann, T. Kreuz, C. Rieke, R.G. Andrzejak, A. Kraskov, P. David, C.E. Elger, K. Lehnertz, On the predictability of epileptic seizures, Clinical Neurophysiology 116 (2005) 569–587. [36] F. Mormann, R.G. Andrzejak, C.E. Elger, K. Lehnertz, Seizure prediction: the long and winding road, Brain 130 (2006) 1–20. [37] D.B. Murray, Forecasting a chaotic time series using an improved metric for embedding space, Physica D: Nonlinear Phenomena 68 (1993) 318–325. [38] A. Pikovsky, M. Rosenblum, J. Kurths, Synchronization a Universal Concept in Nonlinear Sciences, Cambridge Univ Press, 2001. [39] B. Schelter, M. Winterhalder, T. Maiwald, A. Brandt, A. Schad, A. Schulze-Bonhage, J. Timmer, Testing statistical significance of multivariate time series analysis techniques for epileptic seizure prediction, Chaos 16 (2006) 013108. [40] F. Shayegh, S. Sadri, R. AmirFattahi, A theoretical model for spontaneous seizure generation based on Markov chain process, in: 4th International IEEE/EMBS Conference on Neural Engineering, 2009, pp. 637–640. [41] F. Shayegh, R. Amirfattahi, S. Sadri, K. Ansari-Asl, M.H. Saraaee, Defining a new measure for synchronization of multichannel epileptic depth-EEG signals based on identification of parameters of a computational model, in: Proceedings of the IASTED International Conference Cambridge United Kingdom Computational Bioscience (CompBio 2011), 2011, pp. 344–350. [42] F. Shayegh, S. Sadri, R. Amirfattahi, K. Ansari-Asl, Proposing a two-level stochastic model for epileptic seizure genesis, Journal of Computational Neuroscience 1–15 (2013) (under publication). [43] M. Steriade, Neuronal Substrates of Sleep and Epilepsy, Cambridge University Press, 2003. [44] J. Theiler, Spurious dimensions from correlation algorithms applied to limited time-series data, Physical Review 343 (1986) 2427–2432. [45] F. Takens, Detecting strange attractors in turbulence, in: D.A. Rand, L.-S. Young (Eds.), Dynamical Systems and Turbulence Lecture Notes in Mathematics, 898, Springer-Verlag, 1981, pp. 366–381. [46] F. Wendling, F. Bartolomei, J.J. Bellanger, P. Chauvel, Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition, European Journal of Neuroscience 15 (2002) 1499–1508. [47] J.B. Wolf, H.L. Swift, J.A. Swinney, Vastano, Determining Lyapunov exponents from a time series, Physica D 16 (1985) 285–317.