Chapter 1 Introduction
1.3 Organization of the Thesis
This thesis is organized as follows. Chapter 2 introduces the theory of chaos analysis and related works. The proposed algorithm for seizure prediction is described in Chapter 3. Chapter 4 describes the implementation techniques of seizure circuit design. Finally, the experimental results and discussions are presented in Chapter 5.
Chapter 2
Fundamentals of Seizure Analysis Algorithm
This chapter will introduce the seizure analysis algorithm based on the Chaos Theory. First, we explain what the Electroencephalography (EEG) is, then use the Chaos Theory to model EEG for analysis. After this, two well-known algorithms as short-term Lyapunov exponential and correlation dimension are introduced.
2.1 Electroencephalography
Electroencephalography (EEG) is the measurement of electrical activity produced by the brain as recorded from electrodes. So-called scalp EEG is collected from tens to hundreds of electrodes positioned on different locations at the surface of the head. EEG signals shown in Fig. 2-2 are amplified and digitalized for post processing.
In some situations, such as epileptic studies, when deeper brain activity needs to be recorded with more accuracy than provided by scalp EEG, clinicians use an invasive form of EEG known as intracranial EEG (icEEG) where electrodes are placed directly inside the skull (see Fig. 2-1). In some cases, a grid of electrodes is laid on the external surface of the brain, on dura mater yielding epidural EEG but in other cases, a depth electrode known as subdural EEG (sdEEG) and electrocorticography (ECoG) is placed into brain structures, such as the amygdala or hippocampus. Because of the filtering characteristics of the skull and scalp, icEEG activity has a much higher spatial resolution than surface EEG.
Fig. 2-1 Electrodes positions: contacts in red are chosen from the seizure onset zone and contacts in blue are selected as not involved or involved latest during seizure spread.
Fig. 2-2 EEG recordings
Then, we have to classify the EEG by symptoms for the observation of the relations between the normal states and abnormal states of the brain, and evaluating our algorithm later. The time intervals of the EEG in different states are defined as follows:
(a) Pre-ictal: the period prior to the start of the seizure.
(b) Ictal: the seizure onset.
(c) Inter-ictal: the period between seizures.
(d) Post-ictal: the period after a seizure.
Fig. 2-3 Time intervals of EEG recordings
Seizure onset
2.2 Chaotic Modeling for EEG
Many studies have shown that non-linear analysis could characterize the dynamics of neural network underlying EEG which cannot be obtained with conventional linear approach. In the section, we will explain the properties of non-linear dynamics system and how to build an EEG model for chaotic analysis will be explained in detail.
2.2.1 Introduction to Chaos Theorem
Recall that a Newtonian deterministic system is a system whose present state is fully determined by its initial conditions [4] (at least, in principle), in contrast to a stochastic (or random) system, for which the initial conditions determine the present state only partially, due to noise or other external circumstances beyond our control.
For a stochastic system, the present state reflects the past initial conditions plus the particular realization of the noise encountered along the way. So, in view of classical science, we have either deterministic or stochastic systems.
For a long time, scientists avoided the irregular side of nature, such as disorder in a turbulent sea, in the atmosphere, and in the fluctuation of wild–life populations.
Later, the study of this unusual results revealed that irregularity, nonlinearity, or chaos was the organizing principle of nature.
A modern scientific term deterministic chaos depicts an irregular and unpredictable time evolution of many (simple) deterministic dynamical systems, characterized by nonlinear coupling of its variables. Given an initial condition, the dynamic equation determines the dynamic process, i.e., every step in the evolution.
However, the initial condition, when magnified, reveals a cluster of values within a certain error bound. For a regular dynamic system, processes issuing from the cluster
are bundled together, and the bundle constitutes a predictable process with an error bound similar to that of the initial condition. In a chaotic dynamic system, processes issuing from the cluster diverge from each other exponentially, and after a while the error becomes so large that the dynamic equation losses its predictive power.
For example, in 1960s, Ed Lorenz from MIT created a simple weather model in which small changes in starting conditions led to a marked changes in outcome, called sensitive dependence on initial conditions, or popularly, the butterfly effect (i.e., “the notion that a butterfly stirring the air today in Peking can transform storm systems next month in New York, or, even worse, can cause a hurricane in Texas”). Thus long–range prediction of imprecisely measured systems becomes impossibility.
The character of chaotic dynamics can be illustrated with the logistic map as follows [5] :
1 (1 )
n n n
x + =rx −x , (2.1)
a discrete-time analog of the logistic equation for population growth. Here, xn ≥ is 0 a dimensionless measure of the population in the nth generation, and r≥0 is the intrinsic growth rate. We restrict the control parameter r to the range 0≤ ≤r 4 so that (2.1) maps the interval 0≤ ≤x 1 into itself.
A. Period-Doubling
Suppose we fix r, choose some initial population x0, and then use (2.1) to generate the subsequent xn. For the growth rate, we can consider cases as follows.
(1) When r< , the population always goes extinct: 1 xn → as 0 n→ ∞. (2) When 1< <r 3 the population grows and eventually reaches a non-zero
steady state, called a period-1 cycle.
(3) Table 2-1 shows the results of logistic map of initial condition x0 =0.4 and x0 =0.8, and we can find that even though there is a huge difference between the initial conditions, the two series as shown in Fig. 2-4 converge
to the same value in a moment.
Table 2-1 Results of logistic map (r=2.8) Logistic growth equation (r=2.8)
X(0) X(1) X(2) X(3) … X(22) X(23) X(24) X(25) 0.4 0.672 0.617 0.661 … 0.642 0.643 0.642 0.642 0.8 0.448 0.692 0.596 … 0.643 0.642 0.643 0.642
Fig. 2-4 Logistic map of period-1 cycle
(4) When r=3.14, the population builds up again but now oscillates about the former steady state, alternating between a large population in one generation and a smaller population in the next. This type of oscillation, in which xn repeats every two iterations, is called aperiod-2 cycle. Table 2-2 shows the results of logistic map of initial condition x0 =0.1 and x0 =0.3, and the series as shown in Fig. 2-5 reach the same two states after a while.
Table 2-2 Results of logistic map (r=3.14) Logistic growth equation (r=3.14)
X(0) X(1) X(2) X(3) … X(22) X(23) X(24) X(25) 0.1 0.2826 0.6365 0.7264 … 0.5385 0.7803 0.5382 0.7804 0.3 0.6594 0.7050 0.6527 … 0.7792 0.5402 0.7799 0.5389
Fig. 2-5 Logistic map of period-2 cycle
Further period-doublings to cycles of period 8, 16, 32, . . . , occur as r increases.
Specifically, let rn denote the value of r where a 2n-cycle first appears. Then computer experiments reveal that r1 = , 3 r2 =3.449, r3 =3.54409, …, r∞ =3.5699. The convergence is essentially geometric: in the limit of large n, the distance between successive transitions shrinks by a constant factor (2.2).
1
In fact, the same convergence rate appears no matter what unimodal map is iterated. In this sense, the number δ is universal. It is a new mathematical constant, as basic to period-doubling as n is to circles.
When r r> , the answer turns out to be complicated: For many values of r , ∞ the sequence
{ }
xn never settles down to a fined point or a periodic orbit instead thelong-term behavior is aperiodic. Table 2-3 shows the results of logistic map. It is interesting to note that even though the difference of the initial conditions is only 0.0001, the two series as shown in Fig. 2-6 are totally divergent in a minute.
Table 2-3 Results of logistic map (r=3.9) Logistic growth equation (r=3.9)
X(0) X(1) X(2) X(3) … X(22) X(23) X(24) X(25) 0.4 0.936 0.23362 0.6982 … 0.9184 0.2919 0.8062 0.609 0.4001 0.93607 0.23336 0.6977 … 0.6759 0.8542 0.4856 0.9741
Fig. 2-6 Logistic map of r=3.9
A bifurcation diagram summarizes the above phenomenon (Fig. 2-7). The horizontal axis shows the values of the parameter r while the vertical axis shows the possible values of x.
At r=3.4, the attractor is a period-2 cycle, as indicated by the two branches.
As r increases, both branches split simultaneously, yielding a period-4 cycle. This splitting is the period-doubling bifurcation mentioned earlier. A cascade of further period-doublings occurs as r increases, yielding period-8, period-16, and so on, until at r r= ∞ ≈3.57, the map becomes chaotic and the attractor changes from a
finite to an infinite set of points.
Fig. 2-7 Bifurcation diagram for the Logistic map
B. Basic Terms of Nonlinear Dynamics
Recall that nonlinear dynamics is a language to talk about dynamical systems.
Here, brief definitions are given for the basic terms of this language [4].
z Dynamical system: A part of the world which can be seen as a self–contained entity with some temporal behavior. Mathematically, a dynamical system is defined by its state and by its dynamics.
z Phase space: In mathematics and physics, a phase space, introduced by Willard Gibbs in 1901, is a space in which all possible states of a system are represented, with each possible state of the system corresponding to one unique point in the phase space.
z Attractor: An attractor is a ‘magnetic set’ in the system’s phase space to which all neighboring trajectories converge. More precisely, we define an attractor to be a subset of the phase space with the following properties:
(1) It is an invariant set;
(2) It attracts all trajectories that start sufficiently close to it;
(3) It is minimal (it cannot contain one or more smaller attractors).
A strange attractor shown in Fig. 2-8 is defined to be an attractor that exhibits sensitive dependence on initial conditions. Geometrically, an attractor can be a point, a curve, a manifold, or even a complicated set with a fractal structure known as a strange attractor.
Fig. 2-8 A plot of Lorenz's strange attractor
z Fractal: Roughly speaking, fractals are complex geometric shapes with fine structure at arbitrarily small scales. Usually they have some degree of self-similarity. In other words, if we magnify a tiny part of a fractal, we will see features reminiscent of the whole. Sometimes the similarity is exact;
more often it is only approximate or statistical.
Fractals are of great interest because of their exquisite combination of beauty, complexity, and endless structure. They are reminiscent of natural objects like mountains, clouds, coastlines, blood vessel networks, and even broccoli, in a
way that classical shapes like cones and squares can't match.
z Embedding dimension: The number of variables needed to characterize the state of the system. Equivalently, this number is the dimension of the phase space.
z Fractal dimension: The strange attractors typically have fractal microstructure. The attractor dimension counts the effective number of degrees of freedom in the dynamical system, described by a non integer dimension.
C. The link between EEG and chaos
Within the context of brain dynamics [4], there are suggestions that “the controlled chaos of the brain is more than an accidental by–product of the brain complexity” and that “it may be the chief property that makes the brain different from an artificial intelligence machine”. Namely, Chaos drives the human brain away from the stable equilibrium, thereby preventing the periodic behavior of neuronal population bursting.
The EEG, being the output of a multidimensional system [6], has statistical properties that depend on both time and space. Components of the brain (neurons) are densely interconnected and the EEG recorded from one site is inherently related to the activity at other sites. This makes the EEG a multivariable time series. The analysis of such nonlinear dynamical systems from time series involves state space reconstruction, and we will introduce in the next section.
If prediction becomes impossible, it is evident that a chaotic system can resemble a stochastic system, say a Brownian motion. However, the source of the irregularity is quite different. For chaos, the irregularity is part of the intrinsic dynamics of the system, not random external influences. Usually, though, chaotic
systems are predictable in the short–term. This short–term predictability is useful.
Chaos theory has developed special mathematical procedures to understand irregularity and unpredictability of low–dimensional nonlinear systems. Lyapunov exponents and attractor dimension are some examples. Lyapunov exponents evaluate the sensitive dependence to initial conditions estimating the exponential divergence of nearby orbits, and we will discuss the method later. Correlation dimensions estimate the fractal dimension and will be described in Chapter 3.
2.2.2 Reconstruction of Attractors from Time Series
Roux et al. (in 1983) exploited a surprising data-analysis technique, now known as attractor reconstruction (Packard et al. 1980, Takens 1981). The claim is that for systems governed by an attractor, the dynamics in the full phase space can be reconstructed from measurements of just a single time series.
Construction of the embedding phase space from a data segment ( )x t of duration T is made with the method of delays [6]. The vectors X in the phase space i are constructed as
( ( ), ( ) ( ( 1) ))T
i i i i
X = x t x t +τ …x t + p− τ , (2.3)
where τ is the selected time delay between the components of each vector in the phase space, p is the selected dimension of the embedding phase space, and
ti∈[1,T -(p - 1)τ]. Obviously, the accuracy of computation depends on the sampling step Δt which decides the number of vectors Na within a duration T data segment:
0
( )
1ti = + −t i * tΔ , where i∈
[
1,Na]
, (2.4) where t0 is the initial time point of the fiducial trajectory and coincides with the timepoint of the first data in the data segment of analysis.
The embedding dimension p can be determined from (2.5) if the attractor dimension d is known.
2 1
p≥ d+ (2.5)
The choice of delay τ may also significantly affect the metric characteristics of an attractor. If τ is too small, the ith and the (i+1)th coordinates of a phase point are practically equal to each other. In this case, the reconstructed attractor is situated near the main diagonal of the embedding space, the latter complicating its diagnostics.
When a value for τ is chosen that is too large, the coordinates become uncorrelated, and the structure of reconstructed attractor is lost.
2.3 Related Works
A chaotic attractor is an attractor where, on the average, orbits originating from similar initial conditions (nearby points in the phase space) diverge exponentially fast (expansion process); they stay close together only for a short time. If these orbits belong to an attractor of finite size, they will fold back into it as time evolves (folding process). The Lyapunov exponents measure the average rate of expansion and folding that occurs along the local eigen-directions within an attractor in phase space. For an attractor to be chaotic, the largest Lyapunov exponent (LLE) must be positive.
As we mentioned before, a relevant time scale should always be used in order to quantify the physiological changes occurring in the brain. Furthermore, the brain being a nonstationary system, algorithms used to estimate measures of the brain dynamics should be capable of automatically identifying and appropriately weighing existing transients in the data.
Iasemidis et al. developed a method [7] for estimation of short-term Lyapunov exponents (STL), an estimate of LLE for nonstationary data. It is well-known and
widely used in many researches. Here we will take an epileptic seizure prediction system, proposed by L. D. Iasemidis et al. in the recent years, for example to explain the STL in detail.
The short-term Lyapunov exponent (STL) is defined as:
,
based on the reconstruction of attractors from time series, discussed in the last section, where: the evolution of this perturbation after time Δt.
z Δt is the evolution time for δXij, that is, the time one allows δXij to evolve in the phase space.
Fig. 2-9 Displacement vectors in the fiducial trajectory
L1
The crucial parameter is the adaptive estimation in time and phase space of the magnitude bounds of the candidate displacement vector to avoid catastrophic replacements. The improvement in the estimates of L can be achieved by using the proposed modifications.
z For L to be a reliable estimate of STL, the candidate vector ( )X tj should be chosen such that the previously evolved displacement vector
( 1),i j( )
X t
δ − Δ is almost parallel to the candidate displacement vector
, (0) angular separation between two vectors.
z For L to be a reliable estimate of STL, δXi j, (0) should also be small in magnitude in order to avoid computer overflow in the future evolution within very chaotic regions and to reduce the probability of starting up with points on separatrices. This means,
, (0) ( ) ( ) max
i j i j
X X t X t
δ = − < Δ (2.8)
with Δmax assuming small values.
A typical long-term plot of STL versus time, obtained by analysis of continuous EEG, is shown in Fig. 2-10. This figure shows the evolution of STL at a focal electrode site, as the brain progresses from interictal to ictal to postictal states. There is a gradual drop in STL over approximately 2 hours preceding this seizure. The seizure, 2 minutes in duration, is characterized by a sudden drop in STL values with a
consequent steep rise. Postictal STL values exceed preictal values and slowly approach interictal values. This behavior of STL indicates a gradual preictal reduction in chaoticity, reaching a minimum shortly after seizure onset, and a postictal rise in chaoticity that corresponds to the reversal of the preictal pathological state. There will be more discussions about this character in Chapter 3.
Fig. 2-10 Unsmoothed STL over time (140 min), including a 2-min seizure. [7]
Having estimated the STL temporal profiles at each electrode site, and as the brain proceeds toward the ictal state, the temporal evolution of the stability of each cortical site is quantified. However, since the brain is a system of spatial extent, information about the interactions of its spatial components should also be taken in consideration by the relations of the STL between different cortical sites.
The T- index at time t between electrode sites i and j is then defined as:
( ) ( )
{ }
i , j( )
i , j i j
T E STL t STL t t
N
= − ÷σ , (2.9)
where E
{ }
denotes the average of all absolute differences STL ti( )
−STL tj( )
within a moving window wt
( )
λ defined as:( )
1wt λ = for λ∈[t-N-1,t] and w
( )
λ =0 for λ∉[t-N-1,t]where N is the length of the moving window. Then, σi , j
( )
t is the sample standard deviation of the STL differences between electrode sites i and j within the moving window.A dynamical transition toward a seizure is announced at time t when the * T-indexes of sites over time transits from a value above threshold T1 at times t t'< , to a value below threshold T2 at time t , as shown in Fig. 2-11. *
Fig. 2-11 The T-index curves denoting entrainment 55 min before seizure SZ2 [7].
The method presented achieved amazing results with high prediction sensitivity, and pretty low false prediction rate. More details and comparison results with other algorithms and ours will be listed in Chapter 5.
Chapter 3
Wavelet-Correlation Dimension Seizure Prediction
In this chapter a real-time seizure prediction method based on correlation dimension analysis is presented, including the system architecture, data flow, and algorithms.
3.1 Architecture of Seizure Prediction
Before starting the prediction processing, first we observe the EEG signal whether exists any clue around the seizures. For example, as show in Fig. 3-1, for different clinical states including pre-ictal, ictal, and post-ictal states, the corresponding properties of intracranial EEG recordings are different.
Fig. 3-1 Typical EEG waveforms corresponding to epilepsy:
(a) pre-ictal, (b) ictal, and (c) post-ictal
In the pre-ictal state the EEG signal is of the chaotic nature. As we approach the epileptic seizure the signals are less and less chaotic and take the regular shape. These findings imply that seizures may represent spatiotemporal transitions of the epileptic brain from chaos-to-order-to-chaos. Therefore the chaoticity measure of the signal is a good prognostic of the incoming seizure. In fact, this phenomenon is confirmed by STL in the last chapter, but we try to use another method to prove it.
In this chapter, we would propose a Wavelet-Correlation Dimension based Seizure Prediction system, called WCDSP, as shown in Fig. 3-2. The system is consisted of three primary parts:
z Discrete Wavelet Transform analysis: Wavelet is used to decompose the EEG into several sub-bands.
z Chaos analysis: Correlation dimension is used to measure the EEG complexity.
z Feature extraction: The correlation coefficient is used to be the main feature for prediction rules and to decide the seizure states.
Phase Space Reconstruction Correlation Coefficient Prediction Rules
Fig. 3-2 System architecture of seizure prediction
3.2 Discrete Wavelet Transform
We would introduce the discrete wavelet transform (DWT) in this section, and briefly discuss the properties of the discrete Fourier transform (DFT), the short-time Fourier transform (STFT), and the wavelet transform.
We would introduce the discrete wavelet transform (DWT) in this section, and briefly discuss the properties of the discrete Fourier transform (DFT), the short-time Fourier transform (STFT), and the wavelet transform.