Exploration of the CP decomposition combined with SVD analysis for neonatal oscillatory seizure localization 1

1 of 4

Abstract— A common cause for damage to the neonatal brain is a shortage in the oxygen supply to the brain or asphyxia. Neonatal seizures are the most frequent manifestation of neonatal neurologic disorders. Multichannel EEG recordings allow topographic localization of seizure foci. In this paper we introduce a method based on higher order canonical decomposition combined with subsequent singular value decomposition (SVD) to objectively determine the seizure location in the brain. We illustrate the ability of the method in localizing the seizure activity from disturbing eye artefacts and burst activity in the EEG.

Keywords— CP decomposition, Neonatal EEG, Seizure localization,

I. INTRODUCTION

A common cause for damage to the neonatal brain is a shortage in the oxygen supply to the brain or asphyxia. This damage has important consequences for the development of the child. Neonatal seizures are the most frequent major manifestation of neonatal neurologic disorders [1]. Most seizures are subclinical, being detected only by EEG monitoring [2, 3, 4]. Interpretation of the neonatal EEG requires specific expertise, which is normally not present in the neonatal intensive care unit (NICU). Therefore, there is a need for automated monitoring techniques eliminating constant on-site supervision.

Research supported by a PhD grant of the Institute for the Promotion of Innovation through Science and Technology in Flanders (IWT).

- Research Council KUL:GOA-AMBioRICS, CoE EF/05/006 Optimization in Engineering (OPTEC), IDO 05/010 EEG-fMRI, IOF-KP06/11 FunCopt, several PhD/postdoc & fellow grants;

- Flemish Government: FWO: PhD/postdoc grants, projects, G.0407.02 (support vector machines), G.0360.05 (EEG, Epileptic), G.0519.06 (Noninvasive brain oxygenation), FWO-G.0321.06 (Tensors/Spectral Analysis), G.0302.07 (SVM), G.0341.07 (Data fusion), research communities (ICCoS, ANMMM); IWT: TBM070713-Accelero, TBM-IOTA3, PhD Grants;

- Belgian Federal Science Policy Office IUAP P6/04 (DYSCO, `Dynamical systems, control and optimization’, 2007-2011);

- EU: BIOPATTERN IST 508803), ETUMOUR 2002-LIFESCIHEALTH 503094), Healthagents (IST–2004–27214), FAST (FP6-MC-RTN-035801), Neuromath (COST-BM0601)

- ESA: Cardiovascular Control (Prodex-8 C90242)

Recently, we developed an automated neonatal seizure detection method [5]. This detection method divides the neonatal seizures into two types and for each type a different detection method was developed. The first type is the oscillatory type as it typically consists of a continuous oscillating pattern (Fig. 1B). The second type is the spike train type as it consists of a train of isolated repetitive spikes (Fig. 1A).

After detection, it is important to localize the seizures and determine the brain region involved. To do that, the seizure needs to be complete separated from all other concurrent EEG activity as this activity distorts the localization. In this

Figure 1: (A) Example of a spike train type seizure, (B) example of an oscillatory type seizure.

paper, we focus on extracting the oscillatory seizure activity from the EEG. How to deal with the spike train type seizures will be discussed in a follow-up paper.

We propose a method based on Canonical Decomposition (also known as Parallel Factor Analysis (PARAFAC), often referred to as the CP decomposition) combined with Singular Value Decomposition (SVD) to reliably localize the seizures. The CP decomposition for EEG analysis was introduced fairly recently [6]. Subsequent work [7, 8] has shown that the CP decomposition can provide a reliable

Exploration of the CP decomposition combined with SVD analysis

for neonatal oscillatory seizure localization

1

W. Deburchgraeve,

2

P.J. Cherian,

1

M. De Vos,

3

R.M. Swarte,

2

J.H. Blok,

2

G.H. Visser,

3

P. Govaert and

1

S. Van Huffel

1

Department of Electrical Engineering (ESAT), Katholieke Universiteit Leuven, Kasteelpark Arenberg 10, 3001 Leuven-Heverlee, Belgium.

2

Department of Clinical Neurophysiology, Erasmus MC, University Medical Center Rotterdam, ‘s-Gravendijkwal 230, 3015 CE, Rotterdam, The Netherlands.

3

Department of Neonatology, Sophia Children’s Hospital, Erasmus MC, University Medical Center Rotterdam, Dr. Molewaterplein 60, 3015 GJ, Rotterdam, The Netherlands.

2 of 4

estimate of the epileptic seizure onset zone. We start by introducing the CP decomposition. Next, we describe the combination between CP decomposition and SVD analysis and illustrate the application of the method.

II. DATASET

The data was recorded at the Sophia Children's Hospital (part of the University Medical Center Rotterdam, the Netherlands). The data is from a patient with perinatal asphyxia. The full 10-20 set of electrodes (17 electrodes: Fp1,2, F3,4, C3,4, Cz, P3,4, F7,8, T3,4, T5,6, O1,2) was used. Sampling frequency was 256 Hz. Before analysis, the data was filtered between 0.3 and 30Hz.

III. METHODS

The trilinear CP model for a three-way array

X

(I J K× × ) is given by:

1 R ijk ir jr kr ijk r

x

a b c

e

=

=

∑

+

(1) where R is the number of components used in the CP model and

e

_ijk are the residuals containing the unexplained variation. A pictorial representation of the CP decomposition with R atoms is given in Fig. 2. The CP model is a trilinear model: fixing the parameters in two modes,

x

_ijk is expressed as a linear function of the remaining parameters. The CP decomposition is usually computed by means of an Alternating Least Squares (ALS) algorithm [9]. This means that the least-squares cost function 2 1

( , , )

R r r r r

f A B C

X

A B

C

=

=

−

∑





(2) is minimized by means of alternating updates of one of its matrix arguments, keeping the other two matrices fixed. Because the CP decomposition is a multi-linear decomposition, each update just amounts to solving a classical linear least-squares problem.

CP decomposition is applied on a 3 dimensional tensor, but an EEG is only a two dimensional signal with dimensions channel and time. The third dimension of the tensor is constructed using the continuous wavelet transform of the EEG recording [6]. The tensor then consists of dimensions channel, time and frequency.

The wavelet transform is chosen for its optimal time-frequency resolution. Wavelets resolve high time-frequency components within small time windows and low frequencies in larger time windows. The continuous wavelet

C

at scale

a

and time

t

of a signal

x t

( )

is defined by

Figure 2: Pictorial representation of the CP decomposition with R atoms.

C a time

( ,

)

x t

( ) ( ,

φ

a time t dt

, )

∞

−∞

=

∫

(3)

with

φ

the chosen wavelet.

In this study, we used a biorthogonal wavelet with decomposition order 3. From the scale a of the wavelet, the frequency

f

of the signal can be estimated as:

f

f

c

a t

=

∆

(4)

with

f

_c the center frequency of the wavelet and

∆

t

the sampling period.

We took those frequency components with a center frequency

f

c ranging from 1Hz up to 30Hz with a step of 1Hz. The corresponding scale with each frequency component can be calculated using equation 4.

The length of a neonatal seizure is minimum 10s, but can easily last for several minutes. That is why we need to divide the EEG into several windows. We chose to take windows of 2 s with an overlap of 1 s between consecutive windows. A tensor is constructed of each window with a continuous wavelet as described above. Subsequently, the tensor is decomposed using a 2-component CP decomposition. By checking the Core Consistency Diagnostic, we found the optimal number of components to be two. This decomposition leads to two extracted time vectors (components B1 and B2, see Fig. 2), two localization vectors (A1 and A2) containing the spatial distribution of the extracted time series and two frequency distribution vectors (C1 and C2). For the purpose of localizing the seizures, we are particularly interested in the two localization vectors per window. Both localization vectors of each window are collected in one matrix which we call the A-matrix (Fig. 3). To extract the dominant sources, we perform an SVD analysis on the A-matrix and sort the extracted eigenvectors according to the power of the singular value. The right eigenvector corresponding to the singular value with the highest relative power (S1) is the

spatial distribution of the dominant source in the complete EEG. The relative power of the singular value can be calculated using: 2 2 1 i i N j j

s

S

s

=

=

∑

(5)

3 of 4

Figure 3: Schematic overview of the CP-SVD method. The A-matrix consists of 2 localization vectors for each EEG-window. An SVD analysis is performed on this

A-matrix. The spatial distribution vector in matrix V corresponding to the largest singular value s1 is the

spatial distribution of the dominant source in the complete EEG.

III. RESULTS

We will illustrate the application of the algorithm on the EEG displayed in Fig. 4. There are three important types of activity present in this EEG. There is an oscillatory seizure present on the occipital channel O1. There are horizontal eye movements on the frontal channels Fp1 and Fp2, starting from second 22. And finally, there is a burst (between seconds 10 to 15), with maxima on channels T4 and T6.

Figure 5 shows the result of the CP-SVD analysis described above. Only the three strongest sources are shown (A, B and C). The topographic plots are the spatial distributions corresponding to the first three right eigenvectors in Fig. 3. The value under each topographic plot is the relative power of the corresponding singular value and corresponds to the strength of the source (equation 5). The curves at the bottom are the left eigenvectors (matrix U in Fig. 3) of the SVD analysis of the A-matrix. They correspond to the distribution over the extracted spatial distributions in the A-matrix. The topographic plot in Fig. 5A corresponds to the seizure activity present on channel O1 in the EEG in Fig. 4. The seizure activity is nicely separated from all other activity leading to a good localization vector. The seizure activity is the most dominant source present in the EEG with a relative power of 43%. The curve at the bottom of Fig. 5A shows the distribution of this activity over the extracted spatial distributions in the A-matrix and gives an indication of the time occurrence of the detected activity. From this curve it can be seen that the seizure activity is lowest at the beginning and end of the window, with a dip in the middle. This dip corresponds to the burst activity seen in the EEG in Fig. 4.

Figure 4: Example of an oscillatory seizure predominantly present on channel O1. On the frontal channels Fp1 and Fp2 another oscillation is present due to horizontal eye movements. Between 10 and 15s a burst

of activity is present with maxima on the right temporal channels (T4 and T6).

The topographic plot in Fig. 5B corresponds to the eye activity on the frontal channels present in the EEG in Fig. 4. The distribution over the extracted spatial distributions in the A-matrix shows that the eye activity is strongest at the end of the EEG, high in the beginning and during the burst and weakest when the seizure activity is strongest. This is almost exactly the opposite of the distribution in Fig. 5A.

Figure 5: Results of the O-CP analysis of the EEG shown in Fig. 4. The three most dominant sources are shown here (A, B and C). The topographic plots are the spatial

distributions corresponding to the first three right eigenvectors in Fig. 3. The value under each topographic

plot is the relative power of the singular value and corresponds to the strength of the source. The curves at

the bottom, display the distribution over the extracted spatial distributions in the A-matrix.

4 of 4

The topographic plot in Fig. 5C corresponds to the burst activity on the temporal channels T4-T6. The distribution over the extracted spatial distributions in the A-matrix shows that this activity is indeed strongest during the burst.

IV. CONCLUSION

In this paper, we introduce a new approach to localize oscillatory neonatal seizures. The method uses a 3-dimensional tensor in which the first two dimensions are space and time and the third dimension is constructed by decomposing the EEG into different frequency components. The oscillatory type of seizure is defined to be continuous, meaning that the seizure activity is present during a longer period of time without interruption. The CP-SVD method makes intensive use of this characteristic. The complete seizure segment is divided into a number of overlapping windows with duration of 2 s. The CP-decomposition of the tensor of each window is sensitive for activity that is constant in time, frequency and location. An oscillatory type of seizure meets those requirements and thus the CP-decomposition will be sensitive for the seizure activity and isolate it into a separate trilinear component. The results of all windows are grouped into the A-matrix. The SVD analysis sorts the entries according to how dominant they are in the A-matrix. Because the oscillatory seizure activity is continuous in the complete EEG, the most dominant source of the SVD analysis will be the source of the seizure activity. Thus the O-CP method makes use of the continuity of the seizure at the local window level (CP-decomposition) and the global EEG segment level (dominant source of the SVD analysis).

The proposed CP-SVD method is not applicable to spike train type seizures as neonatal spikes are to short and thus in contrast to oscillatory seizures are not present during the complete 2s interval of analysis. How to deal with these spike trains will be discussed in a follow-up paper.

Future work is addressed to further validation of the method on a large set of patients by comparing the obtained topographic plots with the visual analysis of the EEG by an experienced neurophysiologist. We think these topographic plots can aid the electroencephalographer in reviewing the seizures.

REFERENCES

[1] Volpe, JJ., 2008. Neonatal seizures, Neurology of the newborn, 5th ed. Philadelphia: WB Saunders, 203-244.

[2] Hellstrom-Westas, L., Rosen, I., Swenningsen, N.W., 1985. Silent seizures in sick infants in early life. Diagnosis by continuous cerebral function monitoring. Acta. Paediatr. Scand. 74, 741-748.

[3] Scher, M.S., Painter, M.J., Bergman, I., Barmada, M.A., Brunberg, J., 1989. EEG diagnoses of neonatal seizures: clinical correlations and outcome. Pediatr Neurol 5, 17-24.

[4] Murray, D.M., Boylan, G.B., Ali, I., Ryan, C.A., Murphy, B.P., Connolly, S., 2008. Defining the gap between electrographic seizure burden, clinical expression and staff recognition of neonatal seizures. Arch. Dis. Child. Fetal Neonatal 93, 187-191.

[5] Deburchgraeve, W., Cherian, P.J., De Vos, M., Swarte, R.M., Blok, J.H., Visser, G.H., Govaert, P., Van Huffel, S., 2008. Automated neonatal seizure detection mimicking a human observer reading EEG. Clin. Neurophysiol. 119, 2447-2454.

[6] Miwakeichi, F., Martinez-Montes, E., Valdés-Sosa, P.A., Nishiyama, N., Mizuhara, H., Yamaguchi, Y., 2004. Decomposing EEG data into space-time-frequency components using parallel factor analysis. Neuroimage 22, 1035-1045.

[7] De Vos, M., Vergult, A., De Lathauwer, L., De Clercq, W., Van Huffel, S., Dupont, P., Palmini, A., Van Paesschen, W., 2007a. Canonical decomposition of ictal scalp EEG reliably detects the seizure onset zone. NeuroImage 37, 844-854.

[8] Acar, E., Aykut-Bingol, C., Bingol, H., Bro, R., Yener, B., 2007. Multiway analysis of epilepsy tensors. Bioinformatics 23, 10-18. [9] Smilde, A., Bro, R., Geladi, P., 2004. Multi-way Analysis with

applications in the Chemical Sciences. John Wiley & Sons.

Address of the corresponding author: Deburchgraeve Wouter, MS, Department of Electrical Engineering K.U.Leuven

Kasteelpark Arenberg 10, 3001, Heverlee Email: wouter.deburchgraeve@esat.kuleuven.be