Journal of Neural Engineering

Related content

PAPER

Dynamics of regional brain activity in epilepsy: a cross-disciplinary study on both intracranial and scalp-recorded epileptic seizures

- Discrimination between pre-seismic electromagnetic anomalies and solar activity effects G Koulouras, G Balasis, I Kiourktsidis et al.

To cite this article: George Minadakis et al 2014 J. Neural Eng. 11 026012

- A fractal model of earthquake occurrence: Theory, simulations and comparisons with the aftershock data Pathikrit Bhattacharya, Bikas K Chakrabarti and Kamal

- Self-organized criticality Donald L Turcotte

View the article online for updates and enhancements.

Recent citations - Evaluation of dynamic scaling of growing interfaces in EEG fluctuations of seizures in animal model of temporal lobe epilepsy Claudia Lizbeth Martínez-González et al

This content was downloaded from IP address 128.59.222.107 on 14/09/2017 at 18:54

Journal of Neural Engineering J. Neural Eng. 11 (2014) 026012 (16pp)

doi:10.1088/1741-2560/11/2/026012

Dynamics of regional brain activity in epilepsy: a cross-disciplinary study on both intracranial and scalp-recorded epileptic seizures George Minadakis 1 , Errikos Ventouras 2 , Stylianos D Gatzonis 3 , Anna Siatouni 3 , Hara Tsekou 3 , Ioannis Kalatzis 2 , Damianos E Sakas 3 and John Stonham 1 1

Department of Electronic and Computer Engineering, Brunel University Uxbridge, Middlesex, UB8 3PH, UK 2 Department of Biomedical Engineering, Technological Educational Institute of Athens, Ag. Spyridonos, Egaleo, Athens 12210, Greece 3 Department of Neurosurgery, Athens Medical School, Epilepsy Surgery Unit, Evangelismos Hospital, 45–47 Ipsilantou Str, Athens 10676, Greece E-mail: [email protected] Received 15 November 2013, revised 2 January 2014 Accepted for publication 23 January 2014 Published 10 March 2014 Abstract

Objective. Recent cross-disciplinary literature suggests a dynamical analogy between earthquakes and epileptic seizures. This study extends the focus of inquiry for the applicability of models for earthquake dynamics to examine both scalp-recorded and intracranial electroencephalogram recordings related to epileptic seizures. Approach. First, we provide an updated definition of the electric event in terms of magnitude and we focus on the applicability of (i) a model for earthquake dynamics, rooted in a nonextensive Tsallis framework, (ii) the traditional Gutenberg and Richter law and (iii) an alternative method for the magnitude–frequency relation for earthquakes. Second, we apply spatiotemporal analysis in terms of nonextensive statistical physics and we further examine the behavior of the parameters included in the nonextensive formula for both types of electroencephalogram recordings under study. Main results. We confirm the previously observed power-law distribution, showing that the nonextensive formula can adequately describe the sequences of electric events included in both types of electroencephalogram recordings. We also show the intermittent behavior of the epileptic seizure cycle which is analogous to the earthquake cycles and we provide evidence of self-affinity of the regional electroencephalogram epileptic seizure activity. Significance. This study may provide a framework for the analysis and interpretation of epileptic brain activity and other biological phenomena with similar underlying dynamical mechanisms. Keywords: epileptic seizures, nonextensive statistics, self-affinity, intermittent criticality (Some figures may appear in colour only in the online journal)

EEG EEG-EQ EM-EQ G-R

Abbreviations ES EQ

epileptic seizure earthquake

1741-2560/14/026012+16$33.00

1

electroencephalogram electric earthquake fracto-electromagnetic earthquake Gutenberg and Richter © 2014 IOP Publishing Ltd

Printed in the UK

G Minadakis et al

J. Neural Eng. 11 (2014) 026012

SCP SOC IC

standard Boltzmann–Gibbs (B-G) entropy, and the physical quantity a, which is a constant of proportionality between the energy released during the fracture of fragment and its size r. Herein, we firstly give an updated definition of the electric event, namely the EEG-EQ and we show that the latter nonextensive formula can adequately describe the population of electric events included in the seizure part of both scalprecorded and intracranial EEG recordings. The results are further verified in terms of three additional equations that lead to the power-law distribution of magnitudes expressing the fractal nature of the system under study (Sator and Hietala 2010, Turcotte 1986, Sammis et al 1986, Kaminski and Jaupart 1998, Carpinteri and Pugno 2005): (i) the traditional G-R law (Gutenberg and Richter 1954), (ii) the Sarlis et al (2010) equation which links the nonextensive formula with the traditional G-R law, and (iii) an alternative method (Utsu 1964, Aki 1965) which is based on the maximum likelihood estimation of the b-value included in the traditional G-R formula. The second approach this study draws on refers to the spatiotemporal behavior of the parameters included in the nonextensive formula from the perspective of self-affinity and criticality. First, we apply temporal analysis on the EEG-EQs contained in long-time recordings that include ES, showing the intermittent critical behavior of the regional ES cycle. Applying different thresholds of EEG-EQ magnitudes contained in both types of EEG recordings under study, we show that the same behavior exists on both types of EEG recordings. Such analysis further verifies the common scalefree nature of such phenomena and provides a preliminary indication of the self-affine nature of the regional ES activity in the brain, in terms of nonextensive statistical physics.

Sotolongo-Costa and Posadas (2004) self-organized-criticality intermittent criticality

1. Objective Strong analogies between EQ dynamics and neurodynamics have been reported in the last two decades, suggesting that seismic cycles and neural reverberations can be analyzed within similar mathematical frameworks (Hopfield 1994, Herz and Hopfield 1995, Rundle et al 2002, Kapiris et al 2005, Osorio et al 2010, Eftaxias et al 2012). Specifically, Osorio et al (2010), have shown that a dynamical analogy exists between ESs and EQs using five ‘scale-free’ statistics: the G-R distribution of event sizes, the distribution of intervals, the Omori and inverse Omori laws and the conditional waiting time until the next event. In this direction, Eftaxias et al (2012) using concepts of nonextensive statistical physics showed the underlying analogy not only for the populations of EQs which occurred in different faults and ESs which occurred in different patients (Osorio et al 2010), but also for the populations of (i) fracto-electromagnetic pulses rooted in the fracture of the backbone of strong entities distributed along a single fault sustaining the system, and (ii) the electric pulses included in a single scalp-recorded ES. The similarity of their results has been extended up to the laboratory seismicity, solar flares and magnetic storms. They also showed the existence of common power-law distribution of burst lifetime (duration) in the populations of (i) electric pulses included in a single scalp-recorded ES, (ii) fracto-electromagnetic pulses rooted in the activation of a single fault, and (iii) acoustic pulses in a laboratory. Against these arguments, recent studies have shown that neuronal activity is not consistent with critical states (Bedard et al 2006, Schevon et al 2013), suggesting that the complex extracellular space with its complex resistance and capacitance structure, can account for the apparently 1/ f relationship observed in the EEG. The dynamics of brain networks are typically close to criticality, but depart from the critical state during ESs (Meisel et al 2012). In this line of thought, this work endeavors to further verify the dynamical analogy between EQs and ESs and to elucidate the potential existence of critical ‘scale-free’ behavior of regional brain activity during ESs, building on two different approaches. First, we examine whether the same concepts of nonextensive statistical physics and power-law behaviors exist not only for the population of electric events included in the seizure part of scalp-recorded EEGs, as shown by Eftaxias et al (2012), but also for those included in the seizure part of intracranial ones. We focus on a recently introduced model for EQ dynamics (SotolongoCosta and Posadas 2004) that derives from a generalized statistical mechanics formalism proposed by Tsallis (1988) that describes physical systems characterized by either longrange interactions, long-term memories or a multi-fractal nature. This model leads to a nonextensive G-R type formula for EQ dynamics that includes two parameters: the entropic index q that refers to the deviation of Tsallis entropy from the

2. Approach 2.1. Data collection

Analysis here focuses on six scalp-recorded and two intracranial EEG recordings of ESs, deriving from human subjects diagnosed with epilepsy. The data have been recorded at the Department of Neurosurgery, University of Athens, ‘Evangelismos’ Hospital. The study was approved by the Ethical Committee of ‘Evangelimos’ Hospital by number 20a/18–1–12 and undertaken with the understanding and written consent of each subject. The study conforms with The Code of Ethics of the World Medical Association (Declaration of Helsinki), printed in the British Medical Journal (18 July 1964). The data collection took place in the Epilepsy Telemetry Unit during long-term video/EEG for presurgical evaluation, using a sampling rate of 400 Hz, with 12 bits A/D resolution. The EEG recording apparatus was a GrassTelefactor Beehive Millennium System. The amplifiers used (AS40) had a band-pass of 0.5–100 Hz. No additional filtering was applied. For both scalp and intracranial recordings, the reference was at a scalp electrode located approximately 1 cm behind Cz (‘midline’-like referential montage). The seizures were selected and retrieved by one person from the seizures data bank of Epilepsy Surgery Unit as shown in table 1. 2

G Minadakis et al

J. Neural Eng. 11 (2014) 026012

Table 1. Presentation of patients characteristics and EEG data. m—male; f—female. Seizure types: complex partial (CP).

A/A

Specimen

Sex

ES type

Duration (s)

ES duration (s)

Contact/electrodes used

I1 I2 S1 S2 S3 S4 S5 S6

Intracranial Intracranial Scalp-recorded Scalp-recorded Scalp-recorded Scalp-recorded Scalp-recorded Scalp-recorded

m m f m m m m f

CP CP CP CP CP CP CP CP

325.00 325.00 125.00 250.00 175.00 210.50 094.03 125.00

85.59 80.14 64.27 75.61 74.89 111.54 37.81 48.10

96 of 96 96 of 96 1 of 46 1 of 46 1 of 46 1 of 46 1 of 46 1 of 46

However, all these seizures have been studied, characterized and measured during presurgical discussion by an epilepsy surgery team. Onset and offset times of each seizure are results of a consensus measuring of the team of expert raters. The seizures included in this study were selected retrospectively from a large number of seizures because of their clear start and end.

values of the nonextensive entropic index q characterize the strong correlations developed within the system under study (Tsallis 1988, 2009). The quantification of dynamic changes of the complexity of the system is mainly obtained by time variations of the Tsallis entropy (Sq ), for a given q-parameter. Lower Sq values characterize systems with lower complexity (or higher organization). In seismology, a widely used magnitude–frequency relationship for EQs is the G-R scaling relation (Gutenberg and Richter 1954), given by

2.2. A fragment–asperity model for earthquakes coming from a nonextensive Tsallis formulation

log N(>M) = a − bM,

What we call today as ‘B-G Statistical Mechanics’, is a microscopic expression for the thermodynamic entropy, introduced by Boltzmann and later complemented by Maxwell and Gibbs, aiming to establish a direct link between the mechanical laws and classical thermodynamics. A certain property of the B-G entropy (S) in the context of classical thermodynamics is ‘extensivity’, expressing the proportionality with the number of elements of the system, if the subsystems are statistically (quasi-) independent or alternatively if the correlations within the system are essentially local. In cases where the situation is not of this type and correlations may be far from negligible at all scales, the system is called ‘nonextensive’. Inspired by multi-fractal concepts and structures, Tsallis (1988) proposed a generalized expression of the entropy S in information theory of the B-G statistical mechanics, by introducing an entropic expression characterized by an index q which leads to a nonextensive statistics, as follows:   W  1 q pi , (1) 1− Sq = k q−1 i=1

(2)

where N (>M) is the cumulative number of EQs with a magnitude greater than M; b and a are constants. Recently, Sotolongo-Costa and Posadas (2004) (SCP) introduced a model for EQ dynamics consisting of two rough profiles interacting via fragments filling the gap. According to this model, the motion of the fault planes can be hindered not only due to the irregularities of the two overlapping profiles but also due to the eventual relative position of the fragments filling the gap. These fragments mainly come from the residues of the breakage of the tectonic plates, from where the faults have originated. The mechanism of triggering EQs is established through the combination of the irregularities of the fault planes and the fragments filling the gap between them. The authors proposed a new G-R type formula, assuming that ε ∝ r, where ε is the released relative energy and r the size of the fragment. This was later revised by Silva et al (2006). In the revised version, the fragment–asperity interaction model of SCP was re-analyzed in view of the q-definition of the mean value in the context of Tsallis nonextensive statistics, introduced in the study of Abe and Bagci (2005). A new scale law, ε ∝ r3 , between the released relative energy ε and the size r of fragments was proposed, which is in full agreement with the standard theory of rupture (the well-known seismic moment scaling with rupture length (Vilar et al 2007)). Finally, their approach leads to the following G-R type law for the magnitude distribution of EQs:

where pi are the probability estimations associated with the microscopic configurations, W is their total number, k is Boltzmann’s constant and q is a real number. q → 1 corresponds to the standard extensive B-G statistics. Indeed, = e(q−1) ln(pi ) ∼ 1+(q−1) ln(pi ) in using p(q−1) i the limit q → 1, we recover the usual B-G entropy: S1 = −k W i=1 pi ln(pi ). The entropic index q characterizes the degree of non-additivity reflected in the following pseudo-additivity rule: Sq (A + B) = Sq (A)+Sq (B)+(1−q)Sq (A)Sq (B). The cases q > 1 and q < 1 correspond to sub-additivity, or super-additivity, respectively. It should be stressed that the q-parameter is not a measure of the complexity of the time series but measures the degree of nonextensivity of the corresponding system. Higher observed

log [N (> M)] = log N       2M  2−q 1−q 10 + , (3) log 1 − 1−q 2−q a2/3 where α is the constant of proportionality between the EQ energy, ε, and the size of fragment, r. q is the entropic index that describes the deviation of Tsallis entropy from the traditional Shannon one. 3

G Minadakis et al

J. Neural Eng. 11 (2014) 026012

where M is the lowest binned magnitude in the sample. It has been shown that the latter equation gives unique, unbiased and accurate estimation of b-values (Utsu 1972). Characteristically, Monte Carlo simulations applied on equation (7) have shown that this method is superior to the least-squares estimation of the b-value in the case of N = 50–100 samples at least (Utsu 1964, 1965, 1967a, 1967b). Another strength is that it does not require any fitting process.

The latter equation is not a trivial result, but incorporates the characteristics of nonextensivity into the distribution of EQs by magnitude. Specifically, it has been successfully applied to various fields of study, providing excellent fit on the experimental data related to (i) seismicity generated in various large geographic areas (Sotolongo-Costa and Posadas 2004, Silva et al 2006, Vilar et al 2007, Telesca 2010b, 2010c, 2011, Matcharashvili et al 2011), (ii) the precursory sequence of electromagnetic EQs associated with the activation of a single fault (Papadimitriou et al 2008, Kalimeri et al 2008, Minadakis et al 2012, Eftaxias et al 2012), (iii) solar flares and magnetic storms (Balasis and Eftaxias 2009, Balasis et al 2011), (iv) scalp-recorded ESs (Eftaxias et al 2012). The similarity of q-values obtained from these studies (distributed within the range 1.6–1.8) indicates the universal nonextensive character of such phenomena. It should be stressed that this similarity is not based on an empirical G-R law but on an analytically derived G-R type formula deduced in the framework of the nonextensive Tsallis statistical theory (Tsallis 1988, 2009), which can better describe non-equilibrium systems characterized by great variability. Moreover, Sarlis et al (2010), stated that above some magnitude threshold, the qparameter included in the nonextensive formula (equation (3)) is associated with the b parameter of the G-R formula (equation (2)), by the relation   2−q . (4) b=2× q−1

2.4. The notion of electric EEG earthquake (EEG-EQ)

The notion of EM-EQ, that has been validated in previous studies (Karamanos et al 2006, Eftaxias et al 2000, Papadimitriou et al 2008, Minadakis et al 2012), states that the recorded radiation must be emerging clearly from the EM background. This means that the detected EM radiation has not been significantly absorbed by conducting layers of the crust or the even more conductive sea, implying that useful EM precursor must be associated with an on-land seismic event which is both strong, i.e., with magnitude 6 or greater, and shallow. As opposed to the notion of EM-EQ, the brain is a system that always produces electric events even at the idle state of patient’s activity. Thus, in a manner analogous to the notion of EM-EQ, in this study we seek to identify specific ‘electric events’ in the EEG recording, corresponding to specific electric ‘pulses’ included in seizure, preseizure or true interictal dynamics. According to this line of thought, in this study we will use the term EEG-EQ and the corresponding ‘electric emission’, referring to the electric potential pulses that can be characterized as ‘events’ of interest. A characteristic feature in EEG recordings is that the electric field balances around a specific (‘inactive’) reference point (Ptr ) (Speckmann and Elger 2011), usually forming a Gaussian (normal) distribution (Elul 1969). Indeed, for the EEG recordings under study, this feature can also be observed in figure 1. Specifically, figures 1(a), (b) and (c) suggestively show the histograms of the fluctuations of the electric potential, recorded during the idle time behavior of three different patients, while figures 1(d), (e) and (f) depict the histograms of the fluctuations of the electric potential recorded during the time of ES. The samples I1 and I2 refer to intracranial EEG recordings while the S1 sample refers to scalp-recorded one. It is observed that although the distribution of the electric potentials included in the seizure EEG recordings differs from that of the idle ones, the most frequent occurrence still remains within a zone around the peak of the distribution. Since ac amplifiers were used, it is expected that the peak of the distribution should be at or very near the zero point. We define this zone as a ‘balance zone’. On these grounds, two criteria have been considered in order to establish the start and the end of each candidate electric EEG emission: (i) the electric emission starts when the pulse leaves the balance zone and ends the time where the pulse enters the zone, and (ii) the time duration where the pulse stays in the balance zone is considered as a criterion of state-changing (start-end). Note that for an ideal case (without noise), this balance zone could be limited very close to the ‘inactive’ reference point (Ptr ). Drawing from this feature, we regard as amplitude A of a

2.3. The magnitude–frequency relation for earthquakes: an alternative approach

Utsu (1965) proposed an alternative method for determining the b-value of the G-R formula (2), showing that it is inversely proportional to the mean magnitude M, given the minimum magnitude, Mmin , of the EQ magnitudes. More precisely, if the magnitude–frequency distribution of EQs obey the G-R law, we get the following approximate relation for the moment of order n: N  N (n + 1) , (5) (Mi − Mc )n ≈ (b ln 10)n i=1 where N is the total number of EQs and Mc is the smallest magnitude in the sample. For the case of n = 1, equation (5), is given by log10 e N log10 e = (6) b = N M − Mc M − NM i c i=1 where M is the average magnitude and Mc is the minimum magnitude in a sample. At the same year, Aki (1965) showed that equation (6) is equivalent to the maximum likelihood estimation, and provided the confidence limits for this estimation from a given sample. For the case where there is no uncertainty in magnitude, the approximate 95% confidence √ limit of b is ±1.96(b/ N ). Later on, Utsu (1966) suggested slight modification of equation (6) addressing the use of the binned magnitudes. The final equation reads as follows: log10 e , (7) butsu = M − (Mc − M ) 2 4

G Minadakis et al

200 0 −200

200

El. field (a.u)

400 El. field (a.u)

El. field (a.u)

J. Neural Eng. 11 (2014) 026012

200 0 −200

0 −200

−400 1

2

3

4

5

6

1

7

2

3

4

5

6

7

8

9

1

10

2

3

4

5

6

7

8

9

x 10

x 10

x 10

10 4

4

5

700

7000 500

Balance Zone

Balance Zone

400

N

N

4000

500 400

300

N

5000

Balance Zone

600

6000

300

3000 200

200

2000 100

1000 0

−200

−100

0 electric field

100

0 −400

200

100

−300

0

1.575

1.58

1.585

1.59

200

300

−200

400

1.595

1.6

500 0 −500 −1000 1.59

1.595

1.6

1.605

1.61

1.615

1.62

6

100

200

0 −100

1.625

1.25

1.3

1.35

1.4

1.45

1.5

Balance Zone

1.55 5

x 10 500

350 Balance Zone

300

250

300

100

x 10

300

Balance Zone 400

250

200

300

N

N

200 N

0 electric field

6

x 10

150

150 100

200

100 100

50 0 −1000

−100

(c) S1 scalp-recorded

El. field (a.u)

1000

−1000 1.57

−100 0 100 electric field

(b) I2 intracranial

El. field (a.u)

El. field (a.u)

(a) I1 intracranial

−200

50 −500

0 electric field

500

(d) I1 intracranial seizure

1000

0 −600

−400

−200

0 electric field

200

400

600

(e) I2 intracranial seizure

0 −200

−150

−100

−50

0 50 electric field

100

150

200

(f) S1 scalp-recorded seizure

Figure 1. Histogram distribution of the electric potentials recorded during idle (interictal) ((a), (b), (c)) and seizure ((d), (e), (f)) time windows of three different patients. The EEG recordings are shown above the histograms (x-axis are in samples). The range between the vertical lines indicates the ‘balance zone’ around the ‘inactive’ reference point (Ptr ).

capacitance changes for scalp electrodes, as well as similar phenomena for intracranial electrodes and not by physiological sources, namely the constant conflict between simultaneous excitation and inhibition of neuronal populations. The EEGvideo recorder used, as is the case for the overwhelming majority of commercial EEG recorders used by Epilepsy Units around the world, has as lowest cut-off frequency 0.5 Hz. Our effort was to address at final stage a tool for every day clinical practice in seizure detection. Such a tool ought to be based on capacities and restrictions of commercially available and widely used EEG machines in clinical settings.

candidate electric EEG time series the sequence of k successive units that exceed a threshold (±PtBG ) around the balance point (Ptr ): Aeeg > | ± PtBG (ti )|, i = 1, . . . , k. Since the squared amplitude of the ‘electric emission’ is proportional to their power, the magnitude M of the candidate EEG-EQ is given by the relation

 (8) M = log ε ∼ log [Aeeg (ti )]2 . At this point, it should be taken into account that the ability to record the epileptic patients’ EEG signals with direct current amplifiers, i.e., with a high-pass cut-off frequency at the range of, or lower than, 0.01 Hz, could have provided a window for processing in detail baseline shifts and their possible relation to the mechanisms investigated in the current work (Kim et al 2009, Miller et al 2007, Rampp and Stefan 2012, Vanhatalo et al 2003). Nevertheless, additionally to the fact that in current practice EEG recordings comprise frequencies higher than 0.5 Hz (Kim et al 2009), the experience of the Epilepsy Surgery Unit has shown that EEG recording in the standard presurgical settings with direct current amplifiers, despite precautionary measures, might lead to frequent amplifier saturation, because of uncontrolled direct current shifts caused by electrode

3. Main results 3.1. Evidence of nonextensivity on both scalp-recorded and intracranial epileptic EEG recordings

Contributing to and extending previous observations (Eftaxias et al 2012) in this line of study, first we examine whether equation (3) can describe the distribution of energy of the electric fluctuations included in a single intracranial ES as well as scalp-recorded EEG-ES. Indeed, figures 2(a) and (b) 5

G Minadakis et al

J. Neural Eng. 11 (2014) 026012 Balance zone = ± 0.10 μV

El. field (a.u)

El. field (a.u)

Balance zone = ± 0.10 μV 1000 0

−1000

500 0 −500

−1000

0

50

100

150

200

250

0

50

100

150

time (sec)

200

250

300

time (sec)

nonextensive formula

Silva et. al. formula

0

0 −0.2 −0.4 −0.6

log(G(>M))

log(G(>M))

−0.5

q = 1.8838 ± 0.0004 11 α = (2.2 ± 0.2) x 10

−1

b

est

= 0.26

−0.8

−1.2

q = 1.8281 ± 0.0009 13 α = (2.9 ± 0.3) x 10

−1.4

best = 0.42

−1

−1.6 −1.8 −1.5 2

4

6

8

−2

10

2

4

M

(a)

I1 intracranial (seizure)

(b)

200 0 −200 50

8

10

100

150

I2 intracranial (seizure)

Balance zone = ± 0.10 μV, Sig = 300.00 sec

El. field (a.u)

El. field (a.u)

Balance zone = ± 0.10 μV, Sig = 250.00 sec 400

0

6

M

200 0 −200

200

0

50

100

time (sec)

150

200

250

time (sec)

Silva et. al. formula

Silva et. al. formula

0

0

−0.2 −0.4

−0.8

−0.5

q = 1.2430 ± 0.0192 α = (5.5 ± 0.6) x 1019

log(G(>M))

log(G(>M))

−0.6

−1 −1.2

q = 1.3163 ± 0.0122 18 α = α = (8.3 ± 0.8) x 10 −1

−1.5

−1.4 −1.6 −1.8

−2

−2 1

2

3

4

5

6

7

8

9

1

2

M

(c)

3

4

5

6

7

8

M

I1 intracranial (interictal)

(d)

I2 intracranial (interictal)

Figure 2. (a), (b) Equation (3) was used to fit the distribution electric events included in two intracranial EEG recordings. The part of the

signal included between the vertical solid lines refers to the selected and processed seizure part. (c), (d) Equation (3) was used to fit the distribution of electric events included in 250 and 300 s of interictal parts.

Marquardt 1963) was used for optimizing the fitting process. Herein, N is the total number of the detected EEG-EQs, N(M >) the number of EEG-EQs with magnitude larger than M, G(> M) = N(M >)/N the relative cumulative number of EEG-EQs with magnitude larger than M and α the constant

suggestively show that the distribution of magnitudes of electric pulses included in the ictal part of two intracranial EEG recordings is also described by the formula (3), yielding q = 1.8838 ± 0.0004 and 1.8281 ± 0.0009, respectively. The Levenberg–Marquardt (LM) fitting method (Levenberg 1944, 6

G Minadakis et al

2

2

1.9

1.9

1.8

1.8

1.7

1.7

1.6

1.6

1.5

1.5

q

q

J. Neural Eng. 11 (2014) 026012

1.4

1.4

1.3

1.3 Interictal (120.0s): q ∈ [1.256 ∼ 1.648], qμ =1.485, serr ∈ [0.006 ∼ 0.030]

1.2

Preictal (120.0s): q ∈ [1.389 ∼ 1.670], q =1.505, s μ

1.1

err

Seizure (82.7s): q ∈ [1.510 ∼ 1.905], q =1.758, s μ

1

10

20

30

40

50

err

60

80

Preictal (120.0s): q ∈ [1.363 ∼ 1.764], q =1.571, s μ

1.1

∈ [0.002 ∼ 0.019]

70

Interictal (120.0s): q ∈ [1.153 ∼ 1.653], qμ =1.468, serr ∈ [0.006 ∼ 0.031]

1.2

∈ [0.006 ∼ 0.027]

μ

1

90

10

20

30

#EL

(a)

I1 intracranial

(b)

1.8

1.8

1.7

1.7

1.6

1.6

1.5

1.5

(c)

60

70

80

90

I2 intracranial μ

1.9

preseizure

50

∈ [0.004 ∼ 0.026]

∈ [0.002 ∼ 0.023]

q



interictal

40

err

#EL

1.9

1.4

err

Seizure (85.0s): q ∈ [1.492 ∼ 1.896], q =1.744, s

1.4

seizure

I1 intracranial

interictal

preseizure

(d)

seizure

I2 intracranial

Figure 3. Panels (a), (b) depict the calculated q-parameter included in equation (3), for all the analyzed electrodes. Panels (c), (d) depict the estimated mean qμ parameter (± standard deviation).

the electrodes used in the aforementioned intracranial EEG recordings. The red-squares refer to the seizure part, the blackbullets refer to 120 s of an interictal brain activity recorded approximately half an hour before the seizure and the bluerhombus refer to 120 s right before the start of the seizure. It is observed that the mean nonextensive parameter q is much higher for the seizure part (qμ ≈ 1.75) in contrast to the interictal and preseizure ones. In addition, there is a relative small difference between interictal and preseizure parts. Nevertheless, the results of the investigation concerning the differentiation between the interictal and preseizure (i.e., pre-ictal) states should be considered with caution, and only tentatively, since there is an uncertainty in clearly differentiating the two states. Specifically, it is difficult to decide when a time segment that will be deemed as preseizure has to start, so that its initial part will not belong to the preceding interictal state. Consequently, very short time series are used for the investigation of the preseizure state, therefore precluding the delineation of this state with a sufficient degree of certainty. In figure 4, for repeatability reasons, we provide further evidence showing that the distribution of magnitudes of electric pulses included in the ictal part of all the six

of proportionality between the energy released and the size of fragment (Silva et al 2006). Second, of crucial interest would be to examine whether the nonextensive statistics differentiate seizure, preseizure or true interictal dynamics. However, in order to give a natural meaning for such analysis, we have to consider that analogous to the earth observations, brain is a system that always produce electric events even in its idle state, independently to whether these events can be differently interpreted. Following this line of thought we further used 250 and 300 s of an interictal brain activity, approximately half an hour before the occurrence of the ESs depicted in figures 2(a) and (b). These time-frames have been considered as the relative idle state of the patient’s activity. Figures 2(c) and (d) depict the distribution of magnitudes of electric pulses included in the underlying interictal parts. The calculated q parameters are q = 1.2430 ± 0.0192 and 1.3163 ± 0.0122, respectively. A clear difference is observed here in contrast to those that refer to the ictal state. Extending this type of analysis, we further apply a comparative study between the three underlying states of brain activity: seizure, preseizure and interictal. Figures 3(a) and (b) show the calculated nonextensive parameter q for all 7

G Minadakis et al

J. Neural Eng. 11 (2014) 026012

100 0 −100

Balance zone = ± 0.10 μV

El. field (a.u)

Balance zone = ± 0.10 μV

El. field (a.u)

El. field (a.u)

Balance zone = ± 0.10 μV 200 0

40

60

80

100

120

0

50

100

time (sec) Silva et. al. formula 0

−0.2

−0.2

−0.4

−0.4

est

log(G(>M))

log(G(>M))

b

−1

250

0

50

100

= 0.40

−1.2

150

200

time (sec) Silva et. al. formula 0 −0.2 −0.4

q = 1.720 ± 0.002 11 α = (4.8 ± 0.6) x 10

−0.6

q = 1.832 ± 0.002 α = (2.0 ± 0.3) x 108

−0.8

200

Silva et. al. formula

0

−0.6

150

time (sec)

−0.8

b

−1

est

log(G(>M))

20

0

−1000

−200 0

1000

= 0.78

−1.2

q = 1.801 ± 0.001 14 α = (2.3 ± 0.3) x 10

−0.6 −0.8

b

est

= 0.50

−1

−1.4 −1.4 −1.6 −1.8

−1.6

−1.2

−1.8

−1.4

−2

−1.6

−2 2

3

4

5

6

7

8

9

1

2

3

4

M

(b)

4

100

150

(c)

−100 0

20

40

Silva et. al. formula

60

100 0 −100

80

0

20

40

60

q = 1.715 ± 0.002 α = (1.2 ± 0.1) x 1010

log(G(>M))

−1.5

0

−0.2

−0.2

−0.6

100

120

Silva et. al. formula

0

−0.4

−0.4

best = 0.80

80

time (sec)

Silva et. al. formula

−0.5

10

Scalp-recorded 3

time (sec)

0

8

Balance zone = ± 0.10 μV

0

200

6

M

100

time (sec)

log(G(>M))

2

El. field (a.u)

El. field (a.u)

El. field (a.u)

−100

−1

8

Balance zone = ± 0.10 μV

0

50

7

Scalp-recorded 2

Balance zone = ± 0.10 μV 100

0

6

M

Scalp-recorded 1

(a)

5

q = 1.842 ± 0.002 12 α = (2.5 ± 0.4) x10

log(G(>M))

1

best = 0.37 −0.8

q = 1.747 ± 0.003 14 α = (3.4 ± 0.3) x 10

−0.6 best = 0.68 −0.8 −1

−2

−1 −1.2

−2.5

−1.2 1

2

3

4

5

6

7

8

1

2

3

4

M

(d)

Scalp-recorded 4

5

6

7

8

1

2

M

(e)

Scalp-recorded 5

3

4

5

6

7

8

M

(f)

Scalp-recorded 6

Figure 4. Equation (3) was used to fit the distribution of electric events included in six scalp-recorded EEGs. The part of the signal included

between the vertical solid lines refers to the selected and processed seizure part.

statistics than the one described by either equation (3) or the traditional G-R law. Herein, for the seizure parts depicted in figures 3(a) and (b), the percentage of calculated electrodes with q  1.6 was found to be 86.3% and 90.5%, respectively. This is not an unexpected result since the great variability in the morphology and mode occurrence of seizures, the great variability of non-ictal EEG patterns and the richness and complexity of pathological dynamics (Khan and Gotman 2003, Polychronaki et al 2010, Osorio et al 2010, Eftaxias et al 2012) have already been emphasized. Besides, the variation range of the q-parameter (q  1.6) for most of the analyzed electrodes indicates the universal nonextensive character of such phenomena. Thus, in the following we further examine the consistency of results using three equations related to the traditional G-R formula for EQs.

scalp-recorded EEGs is also described by the formula (3). It is observed that the nonextensive q-parameter varies between 1.715 and 1.842 for all the six recordings depicted in figure 4, with a standard error varying between 0.001 and 0.003 indicating the excellent approximation of the fitting process. It should be noted that the fitting of some analyzed electrodes depicted in figures 3 and 4 was relatively poor: approximately 25% of the calculated electrodes were found with serr > 0.01. Analogous observation has been also verified by Eftaxias et al (2012) who analyzed 100 time intervals of 4096 samples each included in the ictal part of 100 different patients, in terms of nonextensive statistical physics. They found that 84% of the analyzed time intervals can be fitted by equation (3), while the rest seemed to follow quite different 8

G Minadakis et al

J. Neural Eng. 11 (2014) 026012 2

10

1

N

(>M)

N(>M)

10

1

10

Data Points Fitting: M =3.69

Data Points Fitting: Mc=4.17

c −0.26 ± 0.06

N(>M) ~ M

0

10

N

(>M)

0

0

2

4

6

8

10

10

M

(a)

0

~ M−0.41 ± 0.08

2

4

6

8

10

M

I1 intracranial

(b)

I2 intracranial

Figure 5. (a), (b) Equation (2) was used to fit the distribution electric events included in two intracranial EEG recordings depicted in

figure 3. The straight line represents the fitting to the plotted points on a log(N>M ) versus M diagram by the method of least squares. Table 2. Comparative table between equations (2), (4) and (7).

3.2. Evidence by means of Gutenberg and Richter law

Herein, we examine the applicability of equations (2), (7) and (4) to the study of epileptic EEG recordings. Such analysis aims to subsequently enhance the validity of the nonextensive formula, on the one hand, and further elucidate the suggested dynamical analogy between the two phenomena (EQs and ESs), on the other. In figure 5, we present the analysis of the intracranial recordings by means of traditional G-R law (equation (2)) above some magnitude threshold (Mc ). The experimental data were directly fitted to the log(N>M ) versus M distribution of events with magnitude larger than Mc = 3.69 and Mc = 4.17, yielding bGR = 0.26 ± 0.06 and bGR = 0.41 ± 0.08, respectively. The latter observation is consistent with recent studies on premonitory behavior of b-value before EQs which have shown that foreshock temporal EQ series and main shocks are characterized by a much smaller b-exponent (b < 1.0) compared to aftershock sequences (Hainzl et al 2003, Knopoff 2000, Papadimitriou et al 2008). Herein, we note that the value of b in the traditional G-R equation (2) for a given group of EQs is determined usually from the slope of a straight line fitted to the plotted points on a log(N>M ) versus M diagram by the method of least squares. However, such determination of b-value may give inaccurate results that are subjected to the following criticisms (Utsu 1964, 1972). Firstly, the ordinary least-squares method gives too heavy weight to the points for large magnitude. Secondly, different b-values are obtained according to the choice of the length of interval M of magnitude in classifying EQs. Further criticisms can be raised concerning the need for the determination of the magnitude cut-off completeness threshold, which is also a crucial choice for estimating the slope that refers to the linear part of the log(N>M ) versus M distribution of the experimental data. On these grounds the nonextensive equation (3) seems to be superior to the traditional G-R one, since the characteristic value that governs the overall system is the one that corresponds to the smaller magnitude threshold (Minadakis et al 2012, Eftaxias et al

Specimen

bGR , equation (2)

best , equation (4)

butsu , equation (7)

Intracranial 1 Intracranial 2 Scalp-recorded 1 Scalp-recorded 2 Scalp-recorded 3 Scalp-recorded 4 Scalp-recorded 5 Scalp-recorded 6

0.26 ± 0.06 0.41 ± 0.08 0.43 ± 0.11 0.82 ± 0.14 0.46 ± 0.07 0.85 ± 0.09 0.40 ± 0.10 0.59 ± 0.09

0.26 0.42 0.40 0.78 0.50 0.80 0.37 0.68

0.27 0.33 0.38 0.76 0.40 0.64 0.38 0.47

2012). Table 2 shows a comparative list of the calculated b values derived from equations (2), (4) and (7) applied on both the experimental data of the intracranial EEGs depicted in figure 3 and the scalp-recorded ones depicted in figure 4. It is observed that the later equation (7) provides similar b-values indicating its consistency with the results obtained from equations (2) and (4). Expanding on the relevance of equations (2), (7) and (4) to the study of ESs, the analysis of table 2 is applied on all the electrodes used for each one of the two intracranial EEG recordings under study, at two different stages of brain activity: interictal and seizure. More precisely, the relation between the nonextensive q-parameter and the b-value associated with the G-R formula (equation (2)) has been given by Sarlis et al (2010) (see equation (4)). Figure 6 shows six scatter plots of the obtained b-values and q-parameters where the solid lines depict this theoretical relationship. Herein, the bullet-marks refer to the calculated b-values associated to the EEG-EQs contained in the interictal part, while the red-squares refer to those included in the seizure parts. Note that for the estimation of each b-value that derives from equations (2) and (7), we used a cut-off completeness threshold Mc which is mainly determined by the experimental data that refer to the linear part of the log(N>M ) versus M distribution of the G-R law representation. From figure 6, it is observed that the experimental results are 9

G Minadakis et al

J. Neural Eng. 11 (2014) 026012 Scatter plot between q−parameter and b

Scatter plot between q−parameter and b

G−R

G−R

4

4 Theoretical Interictal: 1%(b

Dynamics of regional brain activity in epilepsy: a cross-disciplinary study on both intracranial and scalp-recorded epileptic seizures.

Recent cross-disciplinary literature suggests a dynamical analogy between earthquakes and epileptic seizures. This study extends the focus of inquiry ...
3MB Sizes 0 Downloads 3 Views