Skip Navigation

Bioinformatics 2007 23(13):i10-i18; doi:10.1093/bioinformatics/btm210
This Article
Right arrow Abstract Freely available
Right arrow FREE Full Text (Print PDF) Freely available
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Add to My Personal Archive
Right arrow Download to citation manager
Google Scholar
Right arrow Articles by Acar, E.
Right arrow Articles by Yener, B.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Acar, E.
Right arrow Articles by Yener, B.
Social Bookmarking
 Add to CiteULike   Add to Connotea   Add to Del.icio.us  
What's this?

© 2007 The Author(s)
This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/2.0/uk/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

Multiway analysis of epilepsy tensors

Evrim Acar 1,*, Canan Aykut-Bingol 2, Haluk Bingol 3, Rasmus Bro 4 and Bülent Yener 1

1Computer Science Department, Rensselaer Polytechnic Institute, NY, USA, 2Department of Neurology, Yeditepe University, 3Computer Engineering Department, Bogaziçi University, Istanbul, Turkey and 4Faculty of Life Sciences, Copenhagen University, Copenhagen, Denmark

*To whom correspondence should be addressed.


    ABSTRACT
 TOP
 ABSTRACT
 1 INTRODUCTION
 2 METHODOLOGY
 3 THREE-WAY ANALYSIS OF...
 4 DISCUSSIONS
 REFERENCES
 

Motivation: The success or failure of an epilepsy surgery depends greatly on the localization of epileptic focus (origin of a seizure). We address the problem of identification of a seizure origin through an analysis of ictal electroencephalogram (EEG), which is proven to be an effective standard in epileptic focus localization.

Summary: With a goal of developing an automated and robust way of visual analysis of large amounts of EEG data, we propose a novel approach based on multiway models to study epilepsy seizure structure. Our contributions are 3-fold. First, we construct an Epilepsy Tensor with three modes, i.e. time samples, scales and electrodes, through wavelet analysis of multi-channel ictal EEG. Second, we demonstrate that multiway analysis techniques, in particular parallel factor analysis (PARAFAC), provide promising results in modeling the complex structure of an epilepsy seizure, localizing a seizure origin and extracting artifacts. Third, we introduce an approach for removing artifacts using multilinear subspace analysis and discuss its merits and drawbacks.

Results: Ictal EEG analysis of 10 seizures from 7 patients are included in this study. Our results for 8 seizures match with clinical observations in terms of seizure origin and extracted artifacts. On the other hand, for 2 of the seizures, seizure localization is not achieved using an initial trial of PARAFAC modeling. In these cases, first, we apply an artifact removal method and subsequently apply the PARAFAC model on the epilepsy tensor from which potential artifacts have been removed. This method successfully identifies the seizure origin in both cases.

Contact: acare{at}cs.rpi.edu


    1 INTRODUCTION
 TOP
 ABSTRACT
 1 INTRODUCTION
 2 METHODOLOGY
 3 THREE-WAY ANALYSIS OF...
 4 DISCUSSIONS
 REFERENCES
 
Epilepsy is defined as spontaneous clinical seizures caused by paroxysmal, abnormally synchronous neuronal activity. The electrical symptoms of this abnormal activity are believed to uniquely define and reveal the mechanisms of the underlying abnormal neural function and structure. Localization of the initial seizure discharge is an attempt to find the region that generates the abnormal neural activity. Therefore, the analysis of ictal EEG (scalp or intracranial) is an effective standard for identification of an epileptic focus location.

The majority of the research devoted to automated detection of epileptic events concentrates around spike detection techniques. Although most of these techniques are based on single channel data, in Glover et al. (1989), context knowledge from 16-channel EEG data has been incorporated in building a detection system for epileptic sharp waves. Sharp wave source localization on multi-channel EEG data has also been applied in Flanagan et al. (2002) to determine the areas of interest with epileptic activity.

The main challenge in focus localization is the contamination of EEG with artifacts. Ictal EEG, EEG data recorded during the seizure period of an epileptic patient, is often contaminated with signals originating from eye blinks, eye movements and muscle artifacts. These artifacts undermine the efforts to localize epileptic foci and understand the characteristics of a seizure. Commonly used approaches for artifact removal in ictal EEG are simple filtering techniques and statistical methods such as independent component analysis (ICA) and, lately, canonical correlation analysis (CCA). The filtering methods eliminating EEG activity within certain frequencies may result in loss of significant information about seizure structure in the cases where epileptic signals and artifacts overlap in the frequency domain. Therefore, statistical approaches, especially methods based on ICA, are quite common in artifact removal literature. Some of these statistical techniques rely on the widely accepted assumption of independence between artifacts and epileptic brain signals. Considering EEG recorded at each electrode as a linear mixture of signals originating from independent sources, independent components are extracted using ICA (Delorme et al., 2001; Greco et al., 2005; Urrestarazu et al., 2004; Zhou and Gotman, 2005). The components corresponding to artifacts are later identified by visual inspection (Urrestarazu et al., 2004; Zhou and Gotman 2005) or a semi-automated/automated artifact identification technique based on high-order statistics, i.e., kurtosis, entropy (Delorme et al., 2001; Greco et al., 2005). As an alternative to ICA, a CCA-based artifact removal approach has recently been proposed (Clercq et al., 2005). This technique is similar to ICA-based approaches except for the independence assumption. In Clercq et al. (2005), the underlying idea is the mutual non-correlation between artifacts and epileptic signals.

Artifact removal approaches mentioned so far focus on multi-channel EEG data arranged as a two-way dataset of recordings collected at several electrodes at different time samples. Two-way analysis methods on multi-channel EEG data, however, only allow us to capture temporal and spatial signatures, such as the ones identified by ICA and CCA-based techniques. In order to capture frequency domain information, these methods require one more step, e.g. feature extraction as in LeVan et al. (2006), where several features based on the spectral information content of a component are extracted along with some other features to classify components as artifacts or seizures.

On the other hand, multiway analysis enables us to inspect the information content of the signal in time, frequency and electrode domain simultaneously. In neuroscience, multiway models have been previously employed in studying the effect of a drug on brain activity (Estienne et al., 2001). Results in this study demonstrate that significant information is successfully extracted from a complex drug dataset by a Tucker3 model (Tucker, 1964, 1966) rather than two-way models such as principal component analysis (PCA). Multiway models have become more popular in neuroscience with the idea of decomposing EEG data into space-time-frequency components Miwakeichi et al. (2004). In Miwakeichi et al. (2004), continuous wavelet transformation (CWT) is applied on the signals recorded at each electrode and wavelet-transformed data is arranged as a three-way array with modes Formula Formula x Formula x Formula . The three-way array is then analyzed using a PARAFAC model (Harshman, 1970). Factors in the first, second and third component matrices are used to represent the temporal, spectral and spatial signatures in EEG data, respectively. PARAFAC models with non-negativity constraints are later used in another study on event-related potentials (ERP) to find the underlying structure of brain dynamics (Mørup et al., 2006). These studies have also motivated the application of multiway models in understanding the structure of epileptic seizures in our previous study (Acar et al., 2006), which, to our knowledge, has been the first to analyze epileptic EEG data using a multiway model. We have constructed an epilepsy tensor by rearranging the multi-channel ictal EEG data as a third-order tensor with modes, i.e. time samples, scales (frequency) and electrodes and studied the performance of non-linear and multilinear models in capturing specific epilepsy dynamics. Furthermore, in Acar et al., (2006), the performance comparison of two-way and three-way models in localizing a seizure origin has been presented and the results have suggested better focus localization using multiway models.

1.1 Our contributions
In this extended study, we again construct a third-order tensor in the way we have demonstrated in Acar et al., (2006). While we are particularly interested in localizing the focus of a seizure, in addition to that, our goal is to identify spatial, spectral and temporal signatures of an epileptic seizure as well as an artifact. With a goal of analyzing seizure in these three domains, i.e. time, frequency and electrode, our contributions in this article are as follows:

  1. Epileptic focus localization: We model epilepsy tensors using a PARAFAC model and use PARAFAC components in time, frequency and electrode domain to define a seizure. We localize seizure origins based on the spatial signature of a seizure extracted by a PARAFAC model and identified by a neurologist. We have previously applied a Tucker3 model with orthogonality constraints on the component matrices as a generalization of PCA to three-way arrays in Acar et al., (2006) However, the justification of orthogonality constraints is unclear in neuroscience. In addition to that, interpretation of a Tucker3 model is much harder than that of PARAFAC components due to the flexibility of Tucker3 model.
  2. Artifact extraction: We extract artifacts using a PARAFAC model and use PARAFAC components as spectral, spatial and temporal signatures of an artifact in order to define an artifact.
  3. Artifact removal: Through multilinear subspace analysis, we remove artifacts such as eye movements so that the remaining data does not contain any activity correlated with the artifact.
  4. Dataset: We extend the set of patients and focus on the analysis of 10 seizures from 7 patients with different etiological pathologies.
After a brief introduction on multiway arrays and multilinear models in Section 2, the construction and three-way analysis of an epilepsy tensor are presented in Section 3. We also demonstrate the proposed artifact extraction and removal methods in Section 3. Results, interpretations and future steps are discussed in Section 4.


    2 METHODOLOGY
 TOP
 ABSTRACT
 1 INTRODUCTION
 2 METHODOLOGY
 3 THREE-WAY ANALYSIS OF...
 4 DISCUSSIONS
 REFERENCES
 
Multiway data analysis is an exploratory analysis tool, which captures the multilinear structures in a dataset. Standard two-way methods, e.g. singular value decomposition (SVD) (Gloub and Loan, 1996), commonly applied on matrices often fail to find the underlying structures in multiway arrays. Therefore, data is rearranged as a multiway array and analyzed using multiway models in numerous disciplines including chemometrics, neuroscience, computer vision and social network analysis. In this section, we briefly introduce multiway arrays and common multiway models applied in this article.

2.1 Background
Multiway arrays, often referred to as tensors, are higher-order generalizations of vectors and matrices. Higher-order arrays are represented as Formula , where the order of Formula is N (Formula ) while a vector and a matrix are arrays of order 1 and 2, respectively.

A multiway array can be rearranged as a two-way array by unfolding the slices in a certain mode as shown in Figure 1. This operation is called matricization. Through matricization, matrix multiplication is generalized to multiway arrays. The n-mode product of a multiway array Formula by a matrix Formula is denoted by Formula , and it is a multiway array of size Formula . The n-mode product is defined in Lathauwer et al., (2000) as:


Formula 1

(1)
where, Formula and Formula represent the entries of an N-way and a two-way array, respectively.


Figure 1
View larger version (6K):
[in this window]
[in a new window]
[Download PowerPoint slide]
 
Fig. 1. Matricization of a three-way array in the first mode. A three-way array Figure 1 is unfolded in the first mode and a matrix of size Figure 1, denoted by Figure 1 is formed. Subscript in Figure 1 indicates the mode of matricization.

 
2.2 Multilinear models
The most common multiway models in literature are PARAFAC and Tucker3 models. These models both capture the multilinear structure in data by extracting components that are linear combinations of the original variables. These components are then used to interpret the underlying information content of the data. We briefly discuss the similarity and differences of these two models.

2.2.1 Parallel factor analysis (PARAFAC)
is the extension of bilinear factor models to multilinear data. Mathematically, a PARAFAC model can be represented as the decomposition of a tensor as the linear combination of rank-1 tensors. Let Formula be a three-way array. An R-component PARAFAC model on Formula is given by


Formula 2

(2)
where Formula , Formula and Formula indicate the ith column of component matrices Formula , Formula and Formula , respectively. Formula is a three-way array containing the residuals. The symbol {whitebullet} denotes the outer product of vectors. Vector outer product is defined as follows. Let x, y and z be column vectors of size Formula and Formula and Formula and Formula is a tensor of size Formula , then Formula if and only if Formula . The illustration of a 2-component PARAFAC model on a three-way dataset is also shown in Figure 2.


Figure 2
View larger version (8K):
[in this window]
[in a new window]
[Download PowerPoint slide]
 
Fig. 2. Top: two-component PARAFAC model, where a three-way array Figure 2 is expressed as the sum of two rank-1 tensors and error terms. Figure 2, Figure 2 and Figure 2 are the ith components in the first, second and third mode, respectively. Figure 2 is a three-way array containing the residuals. Bottom: (P, Q, R)-component Tucker3 model, where a three-way array Figure 2 is modeled with component matrices Figure 2, Figure 2 and Figure 2 in the first, second and third mode, respectively. Figure 2 is the core array and Figure 2 contains the error terms.

 
2.2.2 Tucker3
is a more flexible multiway model compared to PARAFAC. The flexibility of a Tucker3 model is due to the existence of a higher-order core array, which enables the interaction of each component in each mode with all components in other modes. On the other hand, in a PARAFAC model, a component in a certain mode can be related to a single component in another mode, i.e. the one with the same index. A (P, Q, R)-component Tucker3 model on a three-way array Formula is given by:


Formula 3

(3)
where Formula , Formula and Formula are the component matrices corresponding to the first, second and third modes, respectively. Formula is the core array and Formula contains the residuals. The illustration of a Tucker3 model on a three-way array is given in Figure 2.

In a PARAFAC model, we extract the same number of components in each mode. When one component is identified as an artifact, for instance in the analysis of an epilepsy tensor, that particular component shows the signature of an artifact in time domain in the first mode, in frequency domain in the second mode and in electrode domain in the third mode. Therefore, we actually identify an artifact using the rank-1 tensor corresponding to it. However, a Tucker3 model decomposes the data using a full core array, Formula , which makes the interpretation of a Tucker3 model more difficult than a PARAFAC model. Any component can interact with any component in another mode in a Tucker3 model, e.g. a component identified as an artifact using the component in electrode mode can have any of the frequency signatures captured by the components in the second mode. This relation is quantified by the core elements and can be interpreted using the core array, which is rather complicated than the interpretation of a PARAFAC model. In addition to that, the main problem might be the rotational ambiguity in a Tucker3 model. Unlike PARAFAC, a Tucker3 model cannot determine component matrices uniquely. When a component matrix is rotated by a rotation matrix, it is possible to apply the inverse of the rotation matrix to the core and still obtain the same model fit. Therefore, a Tucker3 model can determine component matrices only up to a rotation. Consequently, a PARAFAC model is a much more restricted and a simpler model with certain uniqueness properties compared to a Tucker3 model.

2.2.3 Determining the number of components
It is important to extract the right number of components in a multilinear model in order to capture the true underlying structure in data. There are several techniques for determining the number of components, i.e. residual analysis, visual appearance of loadings, the number of iterations of the algorithm, core consistency, etc. In this study, we mostly rely on the core consistency diagnostic (Bro and Kiers, 2003) for finding the number of components of a PARAFAC model.

The core consistency quantifies the resemblance between a Tucker3 core and a PARAFAC core, which is a super-diagonal core or in other words, a vector of coefficients. This diagnostic suggests whether a PARAFAC model with the specified number of components is a valid model for the data. Let Formula be a super-diagonal PARAFAC core such that Formula if Formula , otherwise Formula . Let Formula be a Tucker3 core, where Formula can be non-zero for all i, j, k. Then core consistency diagnostic is defined as follows:


Formula

In Bro and Kiers (2003), a general rule on core consistency diagnostic is introduced. Core consistency above 90% is often used as an indication of the trilinear structure in data and suggests that a PARAFAC model with the specified number of components would be an appropriate model for the data. A core consistency value close to or lower than 50%, on the other hand, demonstrates that PARAFAC-like model would not be appropriate. This diagnostic has been commonly applied in the neuroscience-multiway literature (Estienne et al., 2001; Miwakeichi et al., 2004) often together with other diagnostic tools in order to determine the component number.


    3 THREE-WAY ANALYSIS OF EPILEPSY TENSORS
 TOP
 ABSTRACT
 1 INTRODUCTION
 2 METHODOLOGY
 3 THREE-WAY ANALYSIS OF...
 4 DISCUSSIONS
 REFERENCES
 
Multi-channel EEG consists of a set of signals recorded at several electrodes located on the scalp. We construct an Epilepsy Tensor from multi-channel ictal EEG data for each seizure and model using multiway models for artifact extraction, artifact removal and seizure localization. During our analysis, we use PLS_Toolbox (by Eigenvector Research Inc.) for multiway models and EEGLab (Delorme and Makeig, 2004) for topographic maps across the scalp.

3.1 Dataset
We studied scalp EEG recordings of 10 seizures from 7 patients with different pathology substrates: tumor 4; mesial temporal sclerosis 2; cortical dysplasia 1. Ictal EEG recordings were done with long term video EEG monitoring with scalp electrodes in the epilepsy monitoring units of Yeditepe University Hospital and Marmara University. The recordings were obtained by Telefactor Beehave and Millenium video-EEG monitoring systems with 32 channels. The recording of EEG with referential electrode Formula was used for computational analyses.

The duration of ictal EEG corresponding to each seizure, sampling frequencies and the number of electrodes are summarized in Table 1.


View this table:
[in this window]
[in a new window]

 
Table 1. Dataset of multi-channel ictal EEG

 
3.2 Epilepsy tensor construction and preprocessing
Multi-channel EEG data originally forms a matrix of time samples by electrodes. We center across and scale within the electrode mode before we proceed with the analysis. Then we apply CWT on the signals recorded at each electrode in order to identify the frequency component available at each time sample. As a mother wavelet, we make use of Mexican-hat wavelet. Our selection of the mother wavelet is based on a previous work (Latka et al., 2003) showing that a Mexican-hat wavelet captures epileptic events well. We have not studied the performance of other wavelets such as complex Morlet used in modeling brain signals in Miwakeichi et al., (2004) and Mørup et al., (2006), and it is a future research direction to be explored. After computing wavelet coefficients, we downsample time samples by a certain factor in order to reduce the space complexity of the analysis. The downsampling factors used in the analysis are given in Table 1.

Wavelet transformation of a signal from a single electrode forms the frontal slice corresponding to a particular electrode. Similar to the way described in Miwakeichi et al., (2004), we construct our dataset as a three-way array, Formula with modes: time samples, scales and electrodes (Fig. 3). Each entry of Formula denoted by Formula represents the square of the absolute value of wavelet coefficient at Formula electrode for Formula scale at Formula time sample. Scales and frequencies are often used interchangeably in this article. However, they are different. Scale mode reveals the frequency information but scales are inversely proportional to frequencies.


Figure 3
View larger version (7K):
[in this window]
[in a new window]
[Download PowerPoint slide]
 
Fig. 3. Epilepsy Tensor. Figure 3 represents the multi-channel ictal EEG data which is transformed by CWT using a Mexican-hat wavelet and arranged as a three-way array. Each entry of Figure 3, Figure 3, corresponds to the square of the absolute value of a wavelet coefficient at ith time sample, jth scale and kth electrode.

 
Before multiway analysis of Formula , the data is scaled in scales mode in order to capture the activity in all frequencies rather than only at low frequencies with relatively much higher energy than higher frequencies. Scaling a three-way array within one mode is different than scaling in two-way datasets. Unlike matrices where columns or rows are scaled, in three-way case, whole matrices have to be scaled (Bro and Smilde, 2003).

3.3 Artifact extraction
Once the three-way array Formula with modes: time samples, scales and electrodes, is constructed and preprocessed, we model Formula using an R-component PARAFAC model as in Equation (2). PARAFAC is originally based on Cattell's principle of Parallel Proportional Profiles (Cattell, 1944). The idea behind Parallel Proportional Profiles is that if the same factors are present in two samples under different conditions, then each factor in the first sample is expected to have the same pattern in the second sample but profiles of the factors will be scaled depending on the conditions. When we take a closer look at the idea of parallel proportional profiles, we can observe that a signal from an electrode can be referred to as a sample. These samples are generated by certain underlying sources with spectral, spatial as well as temporal signatures specific to the sources. Each electrode, thus, has a coefficient representing the contribution of the source to the signal (or sample) recorded at that particular electrode. Our aim is to identify the sources, such as an eye artifact, a muscle artifact or an epileptic activity generating a seizure, based on these signatures and relative coefficients of electrodes.

An R-component PARAFAC model on Formula extracts the components Formula , Formula and Formula , for Formula , where these components indicate the signatures of sources in time, frequency and electrode domain, respectively as shown in Figure 4. Consequently, a PARAFAC model can serve as an artifact extraction method by identifying patterns indicative of artifacts. In Figure 4, signatures captured by the first component characterize an eye-artifact. Formula indicates the times the artifact takes place. Formula shows that eye-artifact observed at the specified times (those with high coefficients in Formula ) has a high-scale signature indicating a low-frequency content (1.25–2.5 Hz). Finally, Formula localizes the artifact around electrodes Formula and Formula . Based on the visual analysis of EEG recordings, two neurologists identify the time and the location of an artifact. We observe that the signatures extracted by the model match with the clinically identified time and location of an artifact. The model also gives us information about the spectral properties of an artifact. In our study, we observe that most artifacts have low-frequency content. Other examples of extracted artifact signatures are also given in Figure 6.


Figure 4
View larger version (30K):
[in this window]
[in a new window]
[Download PowerPoint slide]
 
Fig. 4. Two-component PARAFAC model on an epilepsy tensor Figure 4 for a particular seizure. We demonstrate the modeling of an epilepsy tensor by a 2-component PARAFAC model, where the first component corresponds to an eye-artifact while the second component represents a seizure. Top: Temporal (Figure 4), spectral (Figure 4) and spatial signatures (Figure 4) of an eye-artifact. Figure 4 represents the coefficients of time samples, Figure 4 represents the coefficients of scales. Since there is a peak in higher scales on the plot of Figure 4, it indicates that this artifact takes place at lower frequencies. Figure 4 contains the coefficients of electrodes. These coefficients are demonstrated on a colormap using EEGLab (Delorme and Makeig, 2004). Bottom: Temporal (Figure 4), spectral (Figure 4) and spatial signatures (Figure 4) of a seizure. Similar to the first component, Figure 4 represents the coefficients of time samples, Figure 4 represents the coefficients of scales. There is a peak in lower scales on the figure corresponding to Figure 4, which indicates that seizure takes place at higher frequencies. Finally, Figure 4 contains the coefficients of electrodes, which are used to localize the seizure around electrodes T4 and T6.

 
3.4 Seizure origin localization
By pursuing the same discussion on Cattell's idea, when one of the underlying sources in the signals recorded by the electrodes is an epileptic seizure, we can argue that one or more of PARAFAC components can model a seizure in the same way an artifact is modeled. Similar to an artifact, a seizure also has a signature in time, frequency and electrode domain. Once these signatures are extracted using a PARAFAC model, the signature of a seizure in the electrode domain can be used to localize the seizure origin.

Therefore, we can also employ PARAFAC as a model for localizing a seizure origin. We observe in Figure 4 that the second component in time samples mode, Formula , shows an ongoing activity in ictal period. When the second component in the second mode, Formula , is examined, we detect that this ongoing activity in ictal period takes place in low scales indicating a rather high-frequency content (12.5–25 Hz) compared to that of the first component. Eventually, Formula suggests that the activity with described characteristics takes place particularly around electrodes T4 and T6. In fact, this activity is a seizure and the component of a PARAFAC model in the electrode mode localizes the seizure origin. These conclusions are also drawn based on the clinically identified seizure location and time. Since the seizure origin is identified by neurologists as T4 and T6, we expect to observe high coefficients corresponding to these electrodes in the spatial signature of the seizure (we illustrate more examples of seizure origin localization in Fig. 6). Furthermore, the temporal signature corresponding to this activity should have an ongoing activity characterized by high coefficients all through seizure period. The components extracted by the model have these characteristics and therefore, they are considered as the signatures of a seizure. Seizures are also often observed to have relatively higher frequency content compared to artifacts.

3.5 Artifact removal
When artifacts account for most of the variation during ictal period, seizure signatures cannot be captured by PARAFAC components. In this case, we suggest that the variation due to artifacts are removed from the data and Formula contaminated with less artifacts is modeled using a PARAFAC model.

In order to understand the underlying structure of data, we model Formula using a Tucker3 model because a Tucker3 model, unlike the PARAFAC model, is known to reflect the main subspace variation in each mode assuming a multilinear structure. We fit a Tucker3 model as in Equation (3) with large number of components in each mode such that we extract enough components to capture most of the variation in data (over 75%). Using a Tucker3 model with orthonormality constraints in each mode, we model the data with component matrices A, B and C corresponding to time samples, scales and electrode mode, respectively and having orthonormal columns. Components in all modes are extracted in decreasing order of captured variance just like in SVD on matrices. Then based on visual inspection of the components in the electrode mode, first N components with characteristics of a potential artifact are identified. Our goal is to remove the activity associated with these potential artifacts. Similar to the underlying idea in interference subtraction based on subspace analysis in Parra et al., (2005), we make use of multilinear subspace analysis to remove the artifacts. We project the data onto the nullspace of the space spanned by the components characterizing an artifact. The steps of the artifact removal method are described more formally as follows:

  1. Fit a Tucker3 model to Formula with component numbers large enough to capture most of the variation in data. (Suppose that modes of Formula are as given in Fig. 3.)
  2. Pick N components in electrodes mode, which are identified as potential artifacts by visual inspection.
  3. Form matrix Formula with N columns using the N components picked in Step 2.
  4. Construct an orthogonal projector, Formula that projects onto the nullspace of Q, namely onto Formula , where M is an N-dimensional subspace of Formula and the columns of Q are the bases for M:


    Formula 4

    (4)
    where I denotes the identity matrix and Formula represents the pseudoinverse of Q.

  5. Compute Formula , which is the projection of Formula onto Formula as:


    Formula 5

    (5)
    where x3 denotes the product of tensor Formula with matrix Formula in electrode mode.

This artifact removal scheme takes out the effect of an artifact across all frequencies during ictal period from Formula . After removing the artifacts, we remodel Formula using a PARAFAC model and use PARAFAC components to identify the seizure origin and inspect spatial, spectral and temporal signatures of the remaining artifacts and seizure. While the artifact removal process enables the localization of seizure, after artifact removal, signatures of seizures in scale mode suggest that seizures have very low-frequency component. On the other hand, we consistently observe seizure activity at high frequencies (12.5–50 Hz) for the seizures where artifact removal is not needed. We summarize the whole process of multiway analysis of multi-channel ictal EEG in Figure 5.


Figure 5
View larger version (11K):
[in this window]
[in a new window]
[Download PowerPoint slide]
 
Fig. 5. Multiway analysis of multi-channel ictal EEG. After the collection of multi-channel EEG data from epileptic patients, we normalize the data and construct a three-way Figure 5 called an Epilepsy Tensor through wavelet transformation. Figure 5 is then downsampled and scaled in scales mode before multiway analysis. Preprocessed three-way array is modeled using a PARAFAC model for artifact extraction and localization of epileptic focus. Finally, PARAFAC components are compared with clinical findings of epilepsy patients. In the case of artifact removal, preprocessed three-way array is first modeled using a Tucker3 model to detect potential artifacts. Figure 5 is formed as a result of artifact removal and Figure 5 is then modeled using a PARAFAC model to extract the signatures of an artifact and a seizure.

 

Figure 6
View larger version (41K):
[in this window]
[in a new window]
[Download PowerPoint slide]
 
Fig. 6. Some illustrative examples of artifact extraction and seizure origin localization. We present our results corresponding to the electrode mode for four of the seizures when they are modeled using a PARAFAC model. Color scales in the figures are the same as the scale in Figure 4. c1, c2 and c3 stands for the first, second and third component in the electrode mode, respectively. We have not shown the signatures in other modes because of space limitations. The number of components in a PARAFAC model is determined based on core consistency diagnostic (Bro and Kiers, 2003). (A) Seizure 1. First component represents an eye-artifact while the second component localizes the seizure. (B) Seizure 8. First component shows the seizure origin and the second component corresponds to an artifact, which has a low-frequency signature. The third component cannot be visually identified. (C) Seizure 7. This is one of the examples where artifact removal is applied. The components are the PARAFAC components extracted after artifact removal. While the first and third components are the artifacts, the second component represents the seizure. (D) Seizure 10. The first component localizes the seizure around F7 and C3 while the second component corresponds to an artifact.

 
3.6 Parameter selection
In this section, we inspect whether core consistency diagnostic serves as a reliable tool for determining the optimal number of components (optimal in terms of interpretation of the data). We model the Epilepsy Tensor for each seizure with a PARAFAC model using R components, where Formula , until the core consistency drops considerably.

Table 2 demonstrates the core consistency values corresponding to different number of components. We often extract as many components as possible until the core consistency drops because additional components help us capture more variation in data (as long as the captured variation with additional components is significant). For instance, for seizure 10, we can fit a 2-component PARAFAC model with confidence since core consistency indicates that more than two components would not be suitable. On the other hand, for seizure 2, we can choose three or four components. In that case, as we mentioned above, the increase in the explained variation by extracting more components may be used in order to determine the component number. We use all four components and observe that the last two components correspond to artifacts. The first two components are already enough to capture the seizure dynamics and account for most of the variation in the data. Consequently, extracting three components instead of four would not change the seizure origin. Table 2 shows that core consistency gives fairly good indications of the right number of components. However, it is still important to take into account other diagnostics, e.g. residual analysis, visual appearance of the loadings, etc. in order to determine the component number.


View this table:
[in this window]
[in a new window]

 
Table 2. Core Consistency for different component numbers

 

    4 DISCUSSIONS
 TOP
 ABSTRACT
 1 INTRODUCTION
 2 METHODOLOGY
 3 THREE-WAY ANALYSIS OF...
 4 DISCUSSIONS
 REFERENCES
 
We construct an epilepsy tensor for each seizure shown in Table 1 and analyze using the process summarized in Figure 5. When the components extracted by a PARAFAC model are visually inspected, we observe that:
  1. In patients with tumors, seizure localization is restricted to a smaller area and the concordance with visual analysis is high. Artifact extraction on ictal EEG of these patients is also well correlated with clinical observations.
  2. In patients with mesial temporal sclerosis, lateralization is well defined but localization is more widespread.
  3. Artifact removal helps to localize the seizure in the patient with cortical dysplasia.

The visual identification of a seizure onset and its localization has some difficulties. Ictal EEG's are frequently contaminated with movement and muscle artifacts that complicate the analysis of seizure localization and onset. Our multilinear approach and multiway-model based method give us promising results in seizure analysis. Not only the detection of the artifacts but also the localization of all seizures are correlated with the visual analysis. Another observation in this study is that the results of multiway analysis of ictal EEG of tumor patients are more precise.

The development of an automated system capable of localizing epileptic focus would strongly affect the outcome of epilepsy surgeries. Removing or extracting artifacts and exploring the underlying brain dynamics in a seizure are also as crucial as seizure origin localization. They would improve not only the accuracy of the focus localization but also the understanding of the complex structure of epilepsy, which has not yet been fully discovered. Although scalp EEG recordings have limitations in detection and localization of seizure onset, we should advance scalp recording techniques with computer analysis. Multilinear models on scalp EEG recordings show promising results in terms of seizure focus localization and artifact extraction. Experimental results demonstrate that using multilinear models, we can get the definitions of an artifact and a seizure. These definitions are formed by the spectral, spatial and temporal signatures extracted by multiway analysis of multi-channel EEG data arranged as a three-way, named as an Epilepsy Tensor. Based on the accuracy of multilinear models, we conclude that the epilepsy tensor represents the multilinear structure of ictal EEG.

Nevertheless, there still exist many research directions for improving and generalizing the results of this study. First of all, we have so far studied whether multilinear models can be used to model an epileptic seizure using the clinical feedback from neurologists. On the other hand, for a fully automated system, automated detection of artifact and seizure components are essential. Secondly, we need to investigate what other information is removed during an artifact removal process via multilinear subspace analysis in addition to the artifact itself and how this process affects the signatures of the remaining artifacts and a seizure. Another investigation area is the correlation between different pathological substrates and the performance of the proposed model, e.g. more precise results on tumor patients. All these future studies are based on the extension of this study on a larger set of patients with different etiological pathologies. Finally, we would like to compare the performance of PARAFAC-based seizure localization and artifact removal approach with ICA-based methods.

Conflict of Interest: none declared.


    REFERENCES
 TOP
 ABSTRACT
 1 INTRODUCTION
 2 METHODOLOGY
 3 THREE-WAY ANALYSIS OF...
 4 DISCUSSIONS
 REFERENCES
 

    Acar E, et al. Computational analysis of epileptic focus localization. ( (2006) ) Proceedings of the Fourth IASTED International Conference on Biomedical Engineering. 317–322..

    Bro R, Kiers HAL. A new efficient method for determining the number of components in parafac models. J. Chemom, ( (2003) ) 17, : 274–286.[CrossRef].

    Bro R, Smilde AK. Centering and scaling in component analysis. J. Chemom, ( (2003) ) 17, : 16–33.[CrossRef].

    Cattell RB. Parallel proportional profiles and other principles for determining the choice of factors by rotation. Psychometrika, ( (1944) ) 9, : 267–283.[CrossRef].

    Clercq WD, et al. A new muscle artifact removal technique to improve the interpretation of the ictal scalp electroencephalogram. ( (2005) ) Proceedings of the 27th Annual Conference of the IEEE Engineering in Medicine and Biology. 944–947..

    Delorme A, Makeig S. Eeglab: an open source toolbox for analysis of single-trial eeg dynamics including independent component analysis. J. Neurosci. Methods, ( (2004) ) 134, : 9–21.[CrossRef][ISI][Medline].

    Delorme A, et al. Automatic artifact rejection for eeg data using high-order statistics and independent component analysis. ( (2001) ) Proceedings of the 3rd International Workshop on ICA. 457–462..

    Estienne F, et al. Multi-way modelling of high-dimensionality electroencephalographic data. Chemom. Intell. Lab. Syst, ( (2001) ) 58, : 59–72.[CrossRef].

    Flanagan D, et al. Computer-aided spatial classification of epileptic spikes. J. Clin. Neurophysiol, ( (2002) ) 19, : 125–135.[CrossRef][ISI][Medline].

    Glover J, et al. Context-based automated detection of epileptogenic sharp transients in the eeg: elimination of false positives. IEEE Trans. Biomed. Eng, ( (1989) ) 36, : 519–527.[CrossRef][ISI][Medline].

    Golub GH, Loan CFV. Matrix Computations, ( (1996) ) Baltimore, MD: The Johns Hopkins University Press..

    Greco A, et al. Semi-automatic artifact rejection procedure based on kurtosis, renyi's entropy and independent component scalp maps. Int. J. Biomed. Sci, ( (2005) ) 1, : 12–16..

    Harshman RA. Foundations of the parafac procedure:models and conditions for an ‘explanatory’ multi-modal factor analysis. UCLA work. pap. phonetics, ( (1970) ) 16, : 1–84..

    Lathauwer LD, et al. A multilinear singular value decomposition. SIAM J. Matrix Anal. Appl, ( (2000) ) 21, : 1253–1278.[CrossRef].

    Latka M, et al. Wavelet analysis of epileptic spikes. Phys. Rev. E Stat. Nonlin. Soft Matter Phys, ( (2003) ) 67, : 105–122..

    LeVan P, et al. A system for automatic artifact removal in ictal scalp eeg based on independent component analysis and bayesian classification. Clin. Neurophysiol, ( (2006) ) 117, : 912–927.[CrossRef][ISI][Medline].

    Miwakeichi F, et al. Decomposing eeg data into space-time-frequency components using parallel factor analysis. NeuroImage, ( (2004) ) 22, : 1035–1045.[CrossRef][ISI][Medline].

    Mørup M, et al. Parallel factor analysis as an exploratory tool for wavelet transformed event-related EEG. NeuroImage, ( (2006) ) 29, : 938–947.[CrossRef][ISI][Medline].

    Parra LC, et al. Recipes for linear analysis of eeg. NeuroImage, ( (2005) ) 28, : 326–341.[Medline].

    Frederiksen N. The extension of factor analysis to three-dimensional matrices. In: Contributions to Mathematical Psychology, ( (1964) ) New York: Holt, Rinehart and Winston. 279–311..

    Tucker LR. Some mathematical notes on three-mode factor analysis. Psychometrika, ( (1966) ) 31, : 279–311.[CrossRef][ISI][Medline].

    Urrestarazu E, et al. Independent component analysis removing artifacts in ictal recordings. Epilepsia, ( (2004) ) 45, : 1071–1078.[CrossRef][ISI][Medline].

    Zhou W, Gotman J. Removing eye-movement artifacts from the eeg during the intracarotid amobarbital procedure. Epilepsia, ( (2005) ) 46, : 409–414.[CrossRef][ISI][Medline].


Add to CiteULike CiteULike   Add to Connotea Connotea   Add to Del.icio.us Del.icio.us    What's this?



This Article
Right arrow Abstract Freely available
Right arrow FREE Full Text (Print PDF) Freely available
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Add to My Personal Archive
Right arrow Download to citation manager
Google Scholar
Right arrow Articles by Acar, E.
Right arrow Articles by Yener, B.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Acar, E.
Right arrow Articles by Yener, B.
Social Bookmarking
 Add to CiteULike   Add to Connotea   Add to Del.icio.us  
What's this?