A. Preprocessing

In the following, we explain the preprocessing steps. More details about these steps can be found in [7].

B. Main Process

1. Source Separation

In this step, we use a reference model of activation, and by maximizing its relative power ratio, we extract the most similar sources to the reference model using GEVD.

a. GEVD principles

GEVD of one pair of symmetric and positive definite matrices (R ¹,R ²) can be stated as follows:

R¹W = R²WΛ

(1)

where Λ is the diagonal matrix of generalized eigenvalues. The generalized eigenvectors build the columns w of matrix W. The ratio

$$\lambda (\mathbf{w})=\frac{{\mathbf{w}}^{\prime }{\mathbf{R}}^1\mathbf{w}}{{\mathbf{w}}^{\prime }{\mathbf{R}}^2\mathbf{w}}$$

(2)

where (·)′ denotes transpose, is known as the Rayleigh quotient, and its maximization is usually performed by GEVD. In fact, the first derivative of λ whose roots are the extremum points of λ writes

$$\nabla \lambda (\mathbf{w})=\frac{2}{{\mathbf{w}}^{\prime }{\mathbf{R}}^2\mathbf{w}}({\mathbf{R}}^1\mathbf{w}-\lambda {\mathbf{R}}^2\mathbf{w}).$$

(3)

Equation (3) shows that the greatest eigenvalue in (1) will be the maximum value of the Rayleigh quotient [11].

Defining various pairs of matrices (R ¹,R ²) and using GEVD for maximizing the Rayleigh quotient has been applied to different problems [12]. Various pairs of matrices (R ¹,R ²) are then computed according to the different criteria. For example, if we are interested in enhancing the signal-to-noise ratio, we can define R ¹ and R ² as the correlation matrix of data during time windows in which the signal has high power, and low power, respectively. Therefore, the temporal sources obtained by GEVD have the maximum power during the first-class time windows and are eventually less noisy.

In [8,13], and [14], it was shown that nonstationary blind source separation (BSS) also can be solved by using GEVD, where matrices R ¹ and R ² are computed on two distinct time windows. Following these works, we may note that common spatial pattern methods in a two-class problem [15] are similar to these methods as well (see [16, Sec. I and II]).

In all these methods, we find the spatial filters using GEVD, but for different criteria. Therefore, the correlation matrices can be defined depending on our objective, and using GEVD, we obtain the temporal sources which have the maximum power related to the first state, class, or time window opposed to the second one [8].

b. Data

Let us denote the observation data as X = [x ₁, …, x _N]′ ^{N ×T}, where x _j = [x_j[1], …, x_j[T]]′, j = 1, …, N is a zero-mean T × 1 matrix corresponding to the T samples recorded on electrode lead j, and N is the number of channels (electrode leads).

Assuming a linear model, X = AS, let us denote A = [a ₁, …, a _N], a _i = [a_i[1], … a_i [N]]′ and S = [s ₁ … s _N]′. a _i (columns of matrix A) and s _i (rows of matrix S) represent the ith spatial pattern and temporal sources, respectively. a_i[j] shows the contribution of the ith temporal source, s _i, to the jth electrode lead.

c. Reference-based source separation method

Here, our objective is to extract the sources related to a reference activation model or IED event. We consider two states, denoted C ¹ and C ², which correspond to the reference and nonreference activations like IED and non-IED, respectively. Denoting , = 1, 2, the set of time samples related to each state, we can define each column of the corresponding segment matrix, X ^{N ×M}, as

χ_i = [x₁[i], …, x_N[i]]′, i ∈ 𝒯^ℓ, M^ℓ = card(𝒯^ℓ)

(4)

where card(.) indicates the cardinality of a set. The correlation matrix of data for each state can be estimated as

$${\widehat{\mathbf{R}}}^{\ell }=\frac{1}{M^{\ell }}{\mathbf{X}}^{\ell }{\mathbf{X}}^{\ell \prime }.$$

(5)

The spatial filters, W, for which the temporal sources S = W′X have maximum similarity with the reference activation state, i.e., maximum variance in the reference state compared to the other state, is computed as

$$\underset{\mathbf{w}}{max}\frac{{\mathbf{w}}^{\prime }{\widehat{\mathbf{R}}}^1\mathbf{w}}{{\mathbf{w}}^{\prime }{\widehat{\mathbf{R}}}^2\mathbf{w}}\text{s}.\text{t}.||\mathbf{w}||=1.$$

(6)

Using (3) for solving (6) leads to GEVD of ( ¹, ²):

$${\widehat{\mathbf{R}}}^1\mathbf{W}={\widehat{\mathbf{R}}}^2\mathbf{W}\mathbf{\Lambda }.$$

(7)

Using W, the spatial patterns, A = (W′)⁻¹, and the temporal sources, S = W′X, are extracted. As explained before, the maximum eigenvalue in (7) is related to the maximum power ratio in (6). We rank the eigenvalues in decreasing order. This implies ranking of the estimated temporal sources according to their resemblance to the reference activation state. The scheme depicted in Fig. 1 demonstrates the block diagram of the R-SS method. The L ¹ reference (or IED) and L ² non-IED segment matrices are denoted as ${\mathbf{X}}_m^{\ell }\in {\mathbb{R}}^{N\times m^{\ell }}$, where m = 1, …, L is the index of each segment. Each column of ${\mathbf{X}}_m^{\ell }$ is , $i\in {\mathcal{T}}_m^{\ell },M_m^{\ell }=\mathit{card}({\mathcal{T}}_m^{\ell })$. ${\mathcal{T}}_m^{\ell }$ indicates the set of time samples of the mth segment related to each state. is equal to 1 and 2 for indicating IED and non-IED, respectively. The correlation matrix for each of these IED and non-IED segment matrices is calculated using (5) by substituting X with ${\mathbf{X}}_m^{\ell }$ and M with $M_m^{\ell }$. These correlation matrices are denoted as ${\mathbf{R}}_m^{\ell }\in {\mathbb{R}}^{N\times N}$. The averages of L ¹ IED correlation matrices and L ² non-IED correlation matrices denoted as R ¹ and R ² are the inputs of GEVD.

Open in a separate window

Fig. 1

R-SS method. ${\mathbf{R}}_m^1$ and ${\mathbf{R}}_m^2$, m = 1, …, L denote the correlation matrices of different IED and non-IED time intervals, respectively. The averages of these matrices denoted as R ¹ and R ² are the inputs of GEVD. Using GEVD, the discriminative temporal sources, s _i, i = 1, …, N, are estimated. These sources are ranked according to their similarity to the IED class.

The output of GEVD is the ranked spatial patterns and temporal sources according to their similarity to reference or IED state.

C. Feature Extraction

In the previous step, the source members of IED class are identified. Now, for the identification of IED electrode leads, we need to transfer this information from source space to the observation space.

According to the linear relationship between observations and sources and to the constraint ||w|| = 1 considered in (2), the probability of activation of each source, s _i, in each observation, x _j (associated with electrode lead j), corresponds to its relative power contribution to x _j

Using (11) and (12), we define the membership probability of each observation for the IED class, p(x _j C ¹), as follows:

One can estimate the IED regions by maximizing (13) as a single-objective function, i.e.,

where e indicates the electrode lead number involved in the IED event. Estimating the IED electrode leads through (14) has the following two problems. First, (14) is a single-objective optimization problem and provides one single IED lead, while most often there exists a network of structurally and functionally connected brain regions [17]; thus, we expect a set of IED leads close to these regions. However, to obtain a set of IED leads, one can select the leads whose IED class membership probability (p(x _j C ¹)) is greater than a given threshold, which provides threshold dependent results. Second, in (13), the probability of activation of different sources of IED class (s _i, i = 1, …, i ^*) in each electrode lead is averaged, while we are interested in optimizing the activation of each source in each electrode. These problems of single-objective optimization methods are addressed in [18] and [19]. Historically, Pareto introduced multiobjective optimization methods in which all the objective functions are considered simultaneously.

A. Simulated Data

The efficiency of the proposed method is evaluated using computer simulations. As illustrated in Fig. 2, three multilead depth electrodes (A, B, and C) are considered, and positioned parallel to each other. We consider eight brain sources: two epileptic sources are placed 1) between electrodes A and B (e₁), and 2) close to electrode C (e₂), and six nonepileptic sources are randomly placed in the volume around the electrodes. The locations of the electrodes and sources are kept constant in simulations and, for simplicity, are restricted to belong to a 2-D plane. The electrical activity of each brain source is represented by a current dipole [10]. The orientations of the six nonepileptic source dipoles are randomly chosen in the 2-D plane and kept constant in simulations, while the orientations of two epileptic source dipoles are assumed as one of the three possible orientations: a tangential orientation (i.e., along the electrode axis), a radial orientation (i.e., orthogonal to the electrode axis), and a “mixed” orientation (i.e., with a 45° angle with the electrode axis). Three simulations (see Fig. 2) are performed where the orientations of two epileptic dipoles are assumed as 1) both orthogonal (D₀), 2) both tangential (D₁), and 3) both mixed (D₂).

Open in a separate window

Fig. 2

Schematic diagram of simulated data. Three electrodes (A, B, and C), six nonepileptic and two epileptic (e₁ and e₂) sources in three different orientations (D₀, D₁, and D₂) are demonstrated. The squares indicate the electrode leads. Each electrode consists of ten equally spaced recording contacts (3.5-mm intercontact spacing). The thick and thin frames demonstrate the first and second layer nodes of Pareto, respectively.

The time-varying magnitudes of the source dipole moments are assumed to represent the global neuronal activity in the source regions. They are obtained from a neural mass model (modified version of the model proposed in [21]). Model parameters (related to neuronal excitability) are tuned such that epileptic source dipoles are assigned an epileptic time-course (spiking activity), while nonepileptic source dipoles are assigned a “normal” time-course (background activity). For epileptic activity, it is assumed that spikes are originated in source e₁ and propagated to source e₂ with a delay of 30–50 ms. The EEG signals, produced by brain sources at all depth-electrode leads, are simulated by solving the “EEG forward problem” (see [10] for details). Briefly, we assume an infinite homogeneous medium with conductivity σ = 33 × 10⁻⁵ S/mm. At each electrode lead , the electric potential V produced at time t by a current dipole is V (t) = (m(t).u _r)/(4πσr ²), where m(t) = m(t)d is the dipole moment (with magnitude m(t) and orientation d), u _r is a unit vector oriented from the source to the electrode lead, r is the distance between them, and the potentials produced by individual sources add up linearly.

The duration of the simulations is 600 s (sampling frequency: f_s = 512 Hz). In practice, for each simulation, 100 IED time intervals (length: 300 samples each), centered on the peaks of IEDs, are extracted. Each IED time interval includes a single IED (here, a spike). For each of the three simulations, 100 non-IED time intervals (same length as the IED intervals) are also extracted.

The method is applied on the three simulated data. The results are shown in Fig. 2, using thick and thin frames for the first and second Pareto layers, respectively. It can be seen that the union of first and second Pareto layers include the closest electrode leads to at least one of the epileptic sources for the three simulated data. There is no major difference between the results related to the three data of different oriented epileptic sources. Electrode lead A₂ is not selected as there are other nonepileptic sources close to this electrode; thus, the contribution of epileptic source e₁ to A₂ is less than its contribution to A₀ and A₁. More discussion on the results in terms of Pareto layers can be found in Section IV.

The proposed method has been evaluated using more sophisticated simulations in 3-D and more number of electrodes (not reported in this paper due to lack of place), which provided congruent results.

To study the effect of signal-to-interference ratio² (SIR) on the proposed method, the simulations as explained in Section III are repeated for different SIR values. Different SIRs are generated by changing the contribution of the nonepileptic sources (generating background activity) to the simulated iEEG. For SIR [−2, +20] dB, we get congruent results in all directions, although the IEDs are not visible in the iEEG signals with the lowest SIR value. More specifically, the identified electrode leads are identical for all SIR values with some changes in the assignment to the first or the second Pareto layer. Indeed, for example, in the simulated data with source orientation D0, for the highest SIR value, all the identified electrode leads (see Fig. 2) are assigned to the first layer. Decreasing SIR causes A2 to be reassigned from the first layer to the second layer, and for the lowest SIR value, B0 is also assigned to the second layer.

Using simulated data presented in this paper, we qualitatively demonstrated that the proposed algorithm works properly. However, a complete study including the analysis of different parameters and their related quantitative evaluations will be done in our future work.

B. Comparison With Other Methods on Actual Data

The iEEG recordings were obtained from five patients suffering from focal epilepsy. The patients underwent presurgery evaluations with iEEG recordings. They are seizure free after resective surgery. Eleven to fifteen semirigid multilead intracerebral electrodes with 0.8-mm diameter were bilaterally implanted in suspected seizure origins based on clinical considerations. The multilead electrodes (Dixi, Besançon, France) include 5, 10, 15, or 18 leads. Each lead has 2-mm length and is evenly spaced with interspace of 1.5 mm. The iEEG were recorded with an audio–video–EEG monitoring system (Micromed, Treviso, Italy) with a maximum of 128 channels and sampled at 512 Hz. The electrode leads were recognized on the patient’s implantation scheme and localized in the Montreal Neurological Institute atlas. Bipolar derivations were considered between adjacent leads within each electrode. The 50 Hz is removed by a fifth-order notch Butterworth filter with 3-dB cutoff frequencies equal to 48 and 52 Hz.

We compare the IED regions estimated by the R-SS method with two methods: visually inspected SOZ (vSOZ) by the epileptologist and leading IED regions (IED) [7].

1. Comparison With vSOZ

In Table I, the IED regions estimated by the R-SS method and Pareto optimization method are reported. Here, we experimentally choose ε = 0.3 from our data, i.e., the second Pareto layer can be considered as the neighborhood of the first Pareto layer if d _max is smaller than or equal to 30% of ||p _z|| (15). The values of d _max/||p _z|| for our five patients are as follows: 0.16, 0.09, 0.01, 0.17, and 0.6. In our experiments, d _max ≤ ε ||p _z|| is satisfied for patients 1–4. Therefore, the Pareto optimal solutions reported in Table I for these patients are the union of the first and second nondominated layer members while, for patient 5, the regions are associated with the first nondominated layer members. Patient 5 is a particular patient with very focused vSOZ, which may explain why the Pareto optimal solution is limited to only the first nondominated layer and eventually the related IED regions are very focused.

TABLE I

Comparison Between the Visually Inspected SOZ (vSOZ) By the Epileptologist, the Leading IED Regions (IED Regions) Estimated Based on dDCG and the IED Regions Estimated By the R-SS Method

P1	antHC	postHC	pHcG	amyg	mTP
vSOZ	×	×	×	×	×
IED	×	×	×	×
R-SS	×	×	×	×	×
P2	antHC	postHC	pHcG	amyg
vSOZ	×	×	×	×
IED	×
R-SS	×
P3	antHC	postHC	pHcG
vSOZ	×	×	×
IED	×	×
R-SS	×	×	×
P4	antHC	postHC	amyg	mTP	entC
vSOZ	×	×	×	×	×
IED	×	×	×		×
R-SS	×	×			×
P5	midInsG
vSOZ	×
IED	×
R-SS	×

Open in a separate window

In the following, we present the results in terms of regions where the electrode leads are located since it is more significant. The abbreviations used in Table I are as follows: amygdala (amyg); anterior/posterior/internal/superior (ant/post/int/sup); entorhinal cortex (entC); hippocampus (Hc); parahippocampal gyrus (pHcG); temporal (T); temporal pole (TP); mesial (m); gyrus (G); middle short gyrus of insula (midInsG); and patient i(P_i ).

The comparison reported in Table I shows the congruency between vSOZ and the estimated IED regions by our method except for the patient P₂. For P₂, there are vSOZs which are not estimated by the R-SS method. It is important to mention that the suggested vSOZs, which we considered as ground truth, are based on EEG and extra clinical information such as semiology. Moreover, since all patients are seizure-free after resective surgery, we deduce that the removed regions included the necessary regions for creating the seizures, but the removed regions might be more than required. These issues make the interpretation of false negatives or sensitivity challenging. However, the estimation of IED regions including the vSOZ (zero false positive, or precision = 100%; see Table II) is valuable since vSOZs are always included in the removed brain regions. Precision and sensitivity are defined later in this section.

TABLE II

Quantitative Comparison Between IED Regions Based on dDCG and IED Regions Based on R-SS

	Precision		Sensitivity		dis (mm)		% ovp		% ovp2
	R-SS	IED	R-SS	IED	R-SS	IED	R-SS	IED	R-SS	IED
P1	100	100	100	80	1.6	8	100	67	100	70
P2	100	100	25	25	0	0	100	100	85	85
P3	100	100	100	67	3.2	4	100	100	100	100
P4	100	100	80	80	2	11.5	100	67	79	79
P5	100	100	100	100	5.3	5.3	100	100	100	100
mean	100	100	81	70	2.4	5.8	100	87	93	87

Open in a separate window

The optimum value in each row for each measure is in bold.

B. Advantages of Using GEVD

In most usual BSS methods, identification of sources related to the desired state is done in two steps. The sources are first estimated, and then, the sources of interest are selected usually by correlation with the reference. On the contrary, in GEVD, the two tasks are done in a unique step, since the estimated sources are ranked according to their similarity to the reference model. Therefore, we only need to identify the number of related sources. In addition, GEVD does not require the assumption of either independence or non-Gaussianity as it is done in independent component analysis (ICA) methods. However, it would obviously lead to decorrelated sources. Since in our application, the sources are not necessarily independent, GEVD is more flexible. Finally, solving GEVD is very fast and only based on second-order statistics contrary to ICA which requires higher order statistics.

The authors gratefully acknowledge P. Kahane, L. Vercueil, O. David, L. Minotti (from Grenoble Institute of Neuroscience, U836 INSERM, and/or Neurology Department, CHUG, Grenoble, France), and F. Wendling (from Inserm 1099, Signal and Image Processing Laboratory, University of Rennes 1) for providing the data, interpretations, and their helpful assistance in this study.