Dynamic Network Connectivity Analysis to Identify Epileptogenic Zones Based on Stereo-Electroencephalography

1. Introduction

Focal epilepsy, in which the origin of epileptic seizures is limited to one hemisphere (Berg et al., 2010), is common and comprising more than 50% of patients with epilepsy (Hauser et al., 1996). Despite great developments in pharmacological treatment, about 30–50% of patients with focal epilepsy cannot be sufficiently controlled with antiepileptic drugs (Beleza, 2009). For these patients, surgical resection of the epileptogenic zones (EZs), the brain areas that are essential for the generation of epileptic seizures (Rosenow and Luders, 2001), may be the only way to suppress or reduce seizures.

EZs can sometimes be adequately identified by a combination of non-invasive techniques, such as analysis of ictal symptomatology, neurological examination, electroencephalography (EEG), magnetoencephalography (MEG), and magnetic resonance imaging (MRI) (Rosenow and Luders, 2001). However, when the EZs cannot be precisely identified, invasive intracranial EEG recordings are required. Stereo-EEG (SEEG) can record neural activities by stereotactic placement of intracranial electrodes within different brain regions, especially in deep areas (Bancaud and Talairach, 1973). In association with long-term video recordings, SEEG allows planning of tailored resections based on individual anatomical and electro-clinical characteristics (Varotto et al., 2012). SEEG is a relatively safe tool and has been considered the gold standard for EZ identification (Cossu et al., 2005).

In clinical practice, EZs are usually identified by visual inspection of ictal SEEG signals. Experienced neurologists can analyze the amplitude, frequency, and relation between signals from each recording channel to identify EZs. However, this is time-consuming and inevitably affected by subjective factors. Moreover, the rather high failure rate of surgical treatment of extratemporal epilepsies (Téllez-Zenteno et al., 2005) underlines that the precise identification of EZs is still an unsolved problem and requires more sophisticated methods to mine further information from SEEG signals.

In recent decades, many computational methods have been proposed to investigate SEEG signals, among which effective connectivity analysis is one of the most widely used. The effective connectivity is defined as the influence that one neuronal system exerts over another, either directly or indirectly (Friston et al., 1993). It is based on Granger causality (Granger, 1969) derived from multivariate autoregressive (MVAR) modeling of multichannel EEG signals and has been successfully used to study the flow of seizures in patients with focal epilepsy (Franaszczuk et al., 1994; Franaszczuk and Bergey, 1998; Varotto et al., 2012). One effective tool to evaluate the effective connectivity is the partial directed coherence (PDC), a frequency-domain measure based on the MVAR coefficients (Baccalá and Sameshima, 2001). In contrast to the directed transfer function (DTF), which is another frequency-domain measure and has also been widely used to study the flow of seizures (Kaminski and Blinowska, 1991; Franaszczuk et al., 1994; Franaszczuk and Bergey, 1998; Ding et al., 2007; Lu et al., 2012; Kim et al., 2014; Basu et al., 2015), PDC can distinguish indirect and direct causalities, the latter being important for localizing EZs. However, both PDC and DTF are based on the MVAR model and are limited to stationary signals. To overcome this limitation, the time-variant MVAR (TVMVAR) model has been introduced (Astolfi et al., 2008; Wilke et al., 2008), and time-variant seizures have been successfully studied using PDC and DTF based on this model, i.e., time-variant PDC (TVPDC) and time-variant DTF (TVDTF) (Wilke et al., 2008, 2009; Leistritz et al., 2013; van Mierlo et al., 2013; Coito et al., 2015). The TVMVAR model is usually estimated using a Kalman filter (Arnold et al., 1998; Schlogl et al., 2000; Milde et al., 2010). As the Kalman filter estimates the system states related to the TVMVAR model's coefficients iteratively, it takes a period of data to follow the real states during the initial estimation procedure (i.e., the adaptation process). Thus, the states estimated during the adaptation process are incorrect and should be discarded (van Mierlo et al., 2013). Actually, the adaptation determines how quickly the Kalman filter can precisely follow the system's dynamics, and an incorrect estimation of the adaptation time will lead to an incorrect estimation of states. In this paper, adaptation is quantitatively evaluated to ensure correct estimation of the time-variant connectivity during propagation of seizures.

Both DTF and PDC measure the information flow from one channel to another for all possible pairs in a network. However, the evaluation of interconnections between all these pairs in a particular frequency band will produce huge matrices of correlation data, which makes a statistical interpretation difficult (Panzica et al., 2013). Graph theory is a promising mathematical framework that is widely used to quantify topological properties of complex networks (Boccaletti et al., 2006), especially the epileptic neural network (Chiang and Haneef, 2014; Haneef and Chiang, 2014). Wilke et al. (2011) found that the betweenness centrality was correlated with the location of resected cortical regions in patients who were seizure-free following surgical intervention. Varotto et al. (2012) analyzed the connectivity pattern in the brain networks of patients with type II focal cortical dysplasia and found that the lesional nodes played a leading role in generating and propagating ictal activities in EEG recordings by acting as the hubs of the epileptic network. van Mierlo et al. (2013) found that the electrode contact with the highest out-degree always lay within the resected brain regions and that the patient-specific connectivity patterns were consistent over the majority of seizures. Sethi et al. (2016) analyzed a network constructed from functional MRI (fMRI) data from patients with polymicrogyria and refractory epilepsy, and found that the polymicrogyric nodes showed significantly increased clustering coefficients and characteristic path lengths compared with the normal contralateral homologous cortical regions. Thus, graph theory is an effective method to identify the highly interconnected hubs in the epileptic network, which is helpful in localizing the EZs in the network. As different graph indices define different types of hubs, the effectiveness of each graph index should be evaluated with a clinical dataset.

In this paper, TVPDC and graph theory are applied to analyze the SEEG signals of patients with refractory focal epilepsy to identify EZs. Three graph indices, namely, in-degree, out-degree, and betweenness centrality, are used to characterize the network dynamics during ictal onset and early propagation. The brain areas where the electrode contacts with relatively high index values are located are considered as the indicative EZs. These indicative EZs are compared with the conclusions from visual analyses of SEEG signals performed by epileptologists, and the effectiveness of each graph index is evaluated. It is found that the in-degree and the betweenness centrality are effective in localizing EZs, but the out-degree is not effective. As this framework requires little prior knowledge other than the seizure onset time, it may provide epileptologists with some objective guidance regarding the abnormal regions.

2. Materials and Methods

2.1. Time-Variant Effective Connectivity

The TVMVAR model is a parametric model that represents a multivariate dataset as a time-variant linear combination of previous samples plus white noise. It is usually used to describe non-stationary process and is characterized by

$\begin{array}{cc}X(n)={\displaystyle \sum _{r=1}^p{A}_r(n)X(n-r)}+E(n)& \left(1\right)\end{array}$

where $X(n)={\left[x_1(n),x_2(n),\dots ,x_d(n)\right]}^T$ is the data vector at time n, d is the number of recording channels, p is the model order, A_r(n) is a time-dependent d×d matrix of model coefficients for time delay r, and $E(n)={\left[e_1(n),e_2(n),\dots ,e_d(n)\right]}^T$ is a vector of white noise at time n.

The Kalman filter algorithm is an efficient method to estimate the coefficients of the TVMVAR model (Arnold et al., 1998; Schlogl et al., 2000; Milde et al., 2010). It takes the coefficients of the TVMVAR model as the state of a state model

$\begin{array}{cc}{S}_n=S_{n-1}+W_n& \left(2\right)\end{array}$

and the TVMVAR model is rewritten as a measurement equation

$\begin{array}{cc}{O}_n=H_n{S}_n+V_n& \left(3\right)\end{array}$

where O_n, S_n, and H_n connect variables of the TVMVAR model and are given by

$\begin{array}{cc}{S}_n=\left[\begin{array}{c}{A}_1{(n)}^T\\ ⋮\\ A_p{(n)}^T\end{array}\right]& \left(4\right)\end{array}$

$\begin{array}{cc}{O}_n=\left[\begin{array}{ccc}{x}_1(n)& \dots & x_d(n)\end{array}\right]=X^T(n)& \left(5\right)\end{array}$

$\begin{array}{cc}{H}_n=\left[\begin{array}{ccc}{O}_{n-1}& \dots & O_{n-p}\end{array}\right]& \left(6\right)\end{array}$

As can be seen, Equation (2) models the state process S_n as a random process with an additive noise W_n, and Equation (3) connects the state process with the observation via a transition matrix H_n and another additive noise V_n. Then the classical computational steps of a Kalman filter are used to estimate S_n, which consists of the TVMVAR coefficients as shown in Equation (4). The initial state for the first iteration is defined as S_p = 0, and the prediction error covariance as P_p = I_dp (the identity matrix of size d×p), and the following steps are then repeated for each recoding sample:

$\begin{array}{cc}{\widehat{S}}_n^{-}={\widehat{S}}_{n-1}& \left(7\right)\end{array}$

$\begin{array}{cc}{P}_n^{-}=P_{n-1}+\overline{W_n},\text{where}\overline{W_n}=E{W}_n{W}_n^T& \left(8\right)\end{array}$

$\begin{array}{cc}{K}_n=P_n^{-}{H}_n^T{(H_n{P}_n^{-}{H}_n^T+\overline{V_n})}^{-1},\text{where}\overline{V_n}=E{V}_n{V}_n^T& \left(9\right)\end{array}$

$\begin{array}{cc}{\widehat{S}}_n={\widehat{S}}_n^{-}+K_n(O_n-H_n{\widehat{S}}_n^{-})& \left(10\right)\end{array}$

$\begin{array}{cc}{P}_n=(I_{dp}-K_n{H}_n)P_n^{-}& \left(11\right)\end{array}$

where E denotes the expectation, ^{^} denotes the a posteriori estimate, and ⁻ denotes the a priori estimate. The covariance matrices V_n are computed recursively as follows (Milde et al., 2010):

$\begin{array}{cc}\overline{V_0}=I_d,\overline{V_n}=\overline{V_{n-1}}(1-\lambda )+\lambda {(O_n-H_n{S}_n)}^T(O_n-H_n{S}_n)& \left(12\right)\end{array}$

where I_d is the identity matrix of size d, and the factor λ works as an update coefficient and has to be chosen between 0 and 1. The matrix $\overline{W_n}$ can be chosen constantly as a weighted identity matrix $\overline{W_n}=\Lambda I_{dp}$ with a weighting constant Λ between 0 and 1. It is appropriate to choose identical λ and Λ (Milde et al., 2010).

In order to select a proper update coefficient (λ, equal to Λ) and determine the adaptation time of the Kalman filter, two Kalman filters with same settings but different start times of the data are used. The filter that starts with the earlier data frame should become adapted earlier, so if the filter that starts later has been adapted, the states estimated by the two filters at the same time should show little difference, and the adaptation time is then defined as the time period that the filter that starts later has taken to become adapted. Here, the difference of the model coefficients is measured using the relative squared error (RSE):

$\begin{array}{cc}\begin{array}{l}\text{RSE}(n)\\ =\frac{{\displaystyle \sum {(\text{vec}([A_1(n),\dots ,A_p(n)])-\text{vec}([B_1(n),\dots ,B_p(n)]))}^2}}{\sum {(|\text{vec}([A_1(n),\dots ,A_p(n)])|+|\text{vec}([B_1(n),\dots ,B_p(n)])|)}^2}\end{array}& \left(13\right)\end{array}$

where A and B are the respective matrices of the model coefficients estimated by the two Kalman filters, and vec refers to converting the matrix to a vector. As the adaptation time is defined by the Kalman filter that starts later, the lag in start time does not matter.

As described above, a set of time-dependent coefficients in the TVMVAR model are estimated. In order to examine the causal relation between signals from different channels in the spectral domain, the Fourier transform of Equation (1) is calculated as follows:

$A(f,n)X(f,n)=E(f,n)\left(14\right)$

where

$\begin{array}{cc}A(f,n)=-{\displaystyle \sum _{r=0}^p{A}_r(n)e^{-i2\pi fr}}& \left(15\right)\end{array}$

with A₀(n) = −I_d (the identity matrix of size d); A(f, n), X(f, n), and E(f, n) are the Fourier transforms of the coefficient matrices A_r(n), the time series X(n), and the residuals E(n), respectively.

Then a spectrum-weighted PDC (swPDC) from channel j to i in a frequency band (f₁, f₂) is calculated as follows (Astolfi et al., 2006; Plomp et al., 2014):

$\begin{array}{cc}{\text{swPDC}}_{ij}(n)=\frac{{\displaystyle {\sum }_{f=f_1}^{f_2}{\text{sPDC}}_{ij}(f,n)S_j(f,n)}}{{\displaystyle {\sum }_{f=f_1}^{f_2}{S}_j(f,n)}}& \left(16\right)\end{array}$

where the lower frequency bound f₁ is used to remove background electrical activity, and the upper frequency bound f₂ is usually limited by the resampling rate. In this study, f₁ and f₂ are set to be 3 and 40 Hz, respectively, which covers the majority of the power in SEEG signals around seizure onset (van Mierlo et al., 2013). S_j(f, n) is the time-variant power spectrum at the source region of j calculated by the short-time Fourier transform (STFT). The sPDC is calculated as follows (Astolfi et al., 2006):

$\begin{array}{cc}{\text{sPDC}}_{ij}(f,n)=\frac{|A_{ij}(f,n){|}^2}{{\displaystyle {\sum }_{k=1}^d|A_{kj}(f,n){|}^2}}& \left(17\right)\end{array}$

where A_ij(f, n) is the element in the ith row and jth column of the matrix A(f, n).

2.2. Graph Analysis

A graph is an abstract of a network consisting of a set of nodes (vertices) and edges. In our study, nodes represent brain regions at which the electrode contacts were located, edges denote interactions between these regions, and the strength of the connectivity between nodes is represented by the swPDC value, which ranges from 0 to 1.

Graphs can be characterized by various measures. One of the simplest quantitative indices of a node is its degree, which indicates the sum of weighted connections from (in-degree deg_in) or toward (out-degree deg_out) all of the other nodes, and represents the most common measure of centrality.

As the number of network nodes differs among patients, the in-degree and out-degree of each node are normalized by the number of network nodes d as follows (Varotto et al., 2012):

$\begin{array}{cc}{\text{deg}}_{\text{in}}(i)=\frac{{\displaystyle {\sum }_{j\ne i}{\text{swPDC}}_{ij}}}{d}& \left(18\right)\end{array}$

$\begin{array}{cc}{\text{deg}}_{\text{out}}(j)=\frac{{\displaystyle {\sum }_{i\ne j}{\text{swPDC}}_{ij}}}{d}& \left(19\right)\end{array}$

Another measure of centrality is the betweenness centrality (bc), defined as the ratio of the number of the shortest paths that pass through a specified node v to the total number of the shortest paths in the network (Wang et al., 2008):

$\begin{array}{cc}\text{bc}(v)={\displaystyle \sum _{i\ne j\ne v\in V}\frac{{\sigma }_{ij}(v)}{{\sigma }_{ij}}}& \left(20\right)\end{array}$

where σ_ij is the number of shortest paths between nodes i and j, and σ_ij(v) is the number of these shortest paths that pass through node v. The betweenness centrality is a measure of the “importance” of each node to information transmission in the network. The nodes that have a relatively high betweenness centrality act as “hubs” in a network, and removal of these nodes will change the network performance significantly (Wilke et al., 2011).

2.3. Patient Data

The SEEG dataset was obtained from seven patients who underwent presurgical evaluation at the Department of Functional Neurosurgery in Renji Hospital (Shanghai, China). The patients included in the study were selected based on the following criteria: focal ictal onset and subsequent resective surgery giving good outcome (with the seizure frequency reduced by at least 50%, and the seizure severity also reduced significantly) during a minimum follow-up of 8 months. The study was approved by the Ethics Committee of Renji Hospital, School of Medicine, Shanghai Jiao Tong University, and all patients gave written informed consent that their clinical data might be used for research purposes. The characteristics of the patients are described in Table 1.

TABLE 1

Table 1. Clinical patient characteristics.

Multi-lead electrodes (HKHS Healthcare, Beijing, China; 5–18 contacts each; diameter 0.8 mm and length 2 mm, and 1.5 mm apart for each contact) were implanted stereotactically. The SEEG signals were recorded using a common reference electrode (Nikon-Kohden system; 128 channels; sampling rate 512 or 1024 Hz) under video and clinical control.

2.4. Identification of the Epileptogenic Zones

For each patient, the ictal SEEG signals around the seizure onset time, which was marked by experienced epileptologists, were selected. As the multivariate causality measures were very sensitive to the data preprocessing (Florin et al., 2010), the dataset was decimated to 128 Hz without any filtering, and its mean was set to be zero. For each resampled epoch, the order of the TVMVAR model was empirically set to be 7.

To select a proper update coefficient and determine the adaptation time of the Kalman filter, a set of update coefficients (10⁻¹, 10⁻², …, 10⁻⁸) were used. The update coefficient with stable and fast adaptation was selected. The time lag of the two Kalman filters was set as 20 s. Then the adaptation time was determined as the time when the RSE reached 0.5%.

Once the update coefficient was selected and the adaptation time was determined, the TVMVAR coefficient matrices from 40 s before until 40 s after the marked seizure onset time were calculated. The coefficient matrices were smoothed with a sliding window (length 0.5 s, overlap 0.125 s). Then the swPDC was calculated in the frequency band (3–40 Hz) based on Equation (16). The three graph indices (in-degree, out-degree, and betweenness centrality) were then calculated based on the swPDC. Then the results within a window from 10 s before until 5 s after the marked seizure onset time were averaged to identify the EZs.

In order to identify the EZs quantitatively, 50% of the maximum was set as the threshold (Wilke et al., 2009, 2011), and the brain areas where the electrode contacts with values exceeding this threshold were located were considered as the indicative EZs. For each graph index, we compared the identified EZs with the epileptologists' results. An overview of the methods used in this study is shown in Figure 1.

FIGURE 1

Figure 1. Overview of the methods utilized in this study.

To evaluate the ability of each index used to identify EZs, we performed a receiver operating characteristic (ROC) curve analysis. Using the results of clinical diagnosis, we set the true positives to be the regions that were identified as EZs clinically and that had high graph values, and the false positives to be the regions that had high graph values but were not identified as EZs clinically. The area under the ROC curve (AUC) measures the effectiveness of the classifier, i.e., AUC ≅ 1 means that the corresponding graph index is an effective classifier to identify EZs.

3. Results

3.1. Adaptation of Kalman Filter

The adaptation of the Kalman filter was evaluated for each patient to optimize the update coefficient and determine the adaptation time. As an example, Figure 2 shows the results for Patient 1, with time 0 corresponding to the seizure onset time. As the figure shows, the highest update coefficient (10⁻¹, black line) leads the Kalman filter to be unstable (the line breaks at −22 s because of the appearance of the numerical value NaN). Taking account of the adaptation speed and the final RSE value, the optimal update coefficient of the Kalman filter for this patient is 10⁻³ (blue line). For the other update coefficients, the Kalman filter either adapts slowly or terminates with high RSE values.

FIGURE 2

Figure 2. Adaptation of the Kalman filter for Patient 1. The optimal update coefficient is 10⁻³ (blue line).

The adaptation time was defined as the time spent until the RSE value reached 0.5%. So, for this case, the adaptation time is 73 s. The optimal update coefficient and adaptation time for the seven patients are listed in Table 2.

TABLE 2

Table 2. The optimal update coefficients and adaptation time.

3.2. Case Study

For each patient, the dynamic networks were constructed by the time-variant method described in Section 2. The three graph indices (in-degree, out-degree, and betweenness centrality), were then calculated to identify the important nodes in the ictal epileptic networks.

For Patient 1, 10 depth electrodes were stereotactically implanted into the hippocampus (Ha, Hp), insula (ISa, ISm, ISp, IIa, IIp), cingulate gyrus (Ca, Cp), and amygdala (Am) (Figure 3A). Each electrode had 12 contacts, numbered from the inside out as 01, 02, …, and 12, and the outermost channel numbered 12 was removed because it was close to the skull. The SEEG signals used for the analysis are shown in Figure 3B.

FIGURE 3

Figure 3. Electrode placement and SEEG signals for Patient 1. (A) Electrode placement. The left panel shows the electrode trajectories on the 3D brain schemes. The right panel shows the trajectory of the electrode Ha superimposed on the sagittal, coronal, and axial MRI views, respectively. (B) The SEEG signals during one seizure. The seizure onset time is set as 0. As can be seen, the seizure started with typical low-amplitude fast activities, preceded by interictal (pre-ictal) spikes.

The corresponding time-varying graph indices (in-degree, out-degree, and betweenness centrality) are shown in Figure 4. In order to compare the in-degree and out-degree, the scales of the color bars are set to be equal. Around the seizure onset time, Hp 01 and Hp 02 showed high in-degree values. The in-degree of Hp 01 retained a high value for quite a long time in the pre-ictal period and increased remarkably at about 20 s before seizure onset. The in-degree of Ha 01 increased at seizure onset and sustained this value for 5 s, then gradually decreased and increased again together with ISa 08 at about 25 s after seizure onset. For the betweenness centrality, the values of Hp 01, Hp 02, and Ha 01 also increased before seizure onset, and the value of ISa 08 increased remarkably at 25 s after seizure onset. These results indicate that the network connection had already changed before seizure onset and that Hp 01, Hp 02, and Ha 01 played a key role in the seizure generation. Moreover, ISa 08 may have played an important role in seizure propagation. However, the out-degree did not reveal any useful information, because the values were almost the same for all the electrode contacts for the whole period studied.

FIGURE 4

Figure 4. Graph analysis results for Patient 1. For each subfigure, the left panel shows the in-degree (A), out-degree (B), and betweenness centrality (C). The middle panel indicates the names of the SEEG electrodes. The right panel shows the mean values of the in-degree (A), out-degree (B), and betweenness centrality (C) within a window from 10 s before to 5 s after seizure onset. The brain regions where the contacts with values exceeding the threshold (50% of maximum, red line) were located were considered to be the identified EZs. The clinically identified EZs are marked by violet shading.

To obtain a comprehensive result for identification of EZs, the mean values of each graph index within a 15 s window (10 s before and 5 s after seizure onset) were calculated, and the results are shown in the right panel of each subfigure in Figure 4.

To identify the EZs, the brain regions where the electrode contacts with values exceeding a threshold (50% of the maximum) were located (the red line in the right panel of each subfigure in Figure 4) were considered as the identified EZs. For Patient 1, the EZs identified by clinical epileptologists were Ha 01, Hp 01, and Hp 02, which are indicated by violet shading in Figure 4. The in-degree index identified Hp 01 and the betweenness centrality identified Ha 01, Hp 01, and Hp 02, which were located in the EZs identified clinically. However, the out-degree index did not indicate any EZs.

For all the other patients, the brain regions where the electrode contacts showing high in-degree and betweenness centrality values were located coincided closely with the EZs identified clinically (Figures 5, 6). However, the out-degree did not indicate any EZs for any of the patients.

FIGURE 5

Figure 5. Graph analysis results for Patients 2, 3, and 4. (A–C) Show the placement of SEEG electrodes for Patients 2, 3, and 4, respectively. (D–F) Show the graph analysis results for Patients 2, 3, and 4 respectively. For each subfigure (D–F), the left panel shows the SEEG signals during one seizure, and the right panel shows the mean values of the in-degree (ID), out-degree (OD), and betweenness centrality (BC) within a window from 10 s before to 5 s after seizure onset. The brain regions where the contacts with values exceeding the threshold (50% of maximum, red line) were located were considered to be the identified EZs. The clinically identified EZs are marked by violet shading.

FIGURE 6

Figure 6. Graph analysis results for Patients 5, 6, and 7. (A–C) Show the placement of SEEG electrodes for Patients 5, 6, and 7, respectively. (D–F) Show the graph analysis results for Patients 5, 6, and 7 respectively. For each subfigure in (D–F), the left panel shows the SEEG signals during one seizure, and the right panel shows the mean values of the in-degree (ID), out-degree (OD), and betweenness centrality (BC) within a window from 10 s before to 5 s after seizure onset. The brain regions where the contacts with values exceeding the threshold (50% of maximum, red line) were located were considered to be the identified EZs. The clinically identified EZs are marked by violet shading.

An ROC analysis was then performed, which allowed us to evaluate the ability of each graph index to identify EZs. As shown in Figure 7, the in-degree and betweenness centrality worked well in identifying EZs (with AUC values near 1). However, the values of the out-degree are all below the diagonal, indicating that the out-degree was not appropriate for identifying EZs in our cases.

FIGURE 7

Figure 7. ROC analysis results for Patients 1–7 (A–G). Both the in-degree (red lines) and betweenness centrality (blue lines) showed good results, whereas the out-degree (black lines) showed poor results (below the diagonal). TPR, true-positive rate; FPR, false-positive rate.

4. Discussion

In this study, we examined the applicability of a method based on time-variant effective connectivity and graph analysis to identify EZs based on SEEG signals of patients with focal refractory epilepsy. The SEEG dataset was obtained from seven patients who had undergone presurgical evaluation and for whom subsequent resective surgery gave a good outcome. By modeling the SEEG signals around seizure onset time with the TVMVAR model, the TVPDC was calculated to construct dynamic networks. Three graph indices, namely, in-degree, out-degree, and betweenness centrality, were calculated to identify EZs. The effectiveness of each graph index was evaluated by comparing the indicative EZs with clinical diagnoses obtained by epileptologists.

In all patients who had successful surgical outcomes, the abilities of the three graph indices to identify the EZs were different. Specifically, we found that the in-degree was most effective, but the out-degree could not identify EZs at all. This is inconsistent with previous studies suggesting that the out-degree is more suitable than the in-degree for identifying EZs (Wilke et al., 2009; Varotto et al., 2012; van Mierlo et al., 2013). The major reason for these conflicting conclusions might lie in the clinical cases that we used. To confirm this, further in-depth research is required. An ideal solution, which will be a future direction for our research, is to collect a large number of clinical cases and perform an analysis in which cases are categorized according to seizure type. Another issue is our use of the PDC rather than the DTF used in the previous studies: the PDC is more suitable for in-degree cases, whereas the DTF is more suitable for out-degree cases (Plomp et al., 2014). Although both DTF and PDC are derived from Granger causality (the MVAR or TVMVAR model coefficients), the normalization procedures are quite different. For column-wise normalized PDC (Baccalá and Sameshima, 2001), the sum of outflows is bounded by 1 for each node, which makes it insensitive to the out-degree since the latter becomes close to 1 if there are tens of nodes in a network (as shown in the right panel of Figure 4B). The same consideration applies for the in-degree calculated from the row-wise normalized DTF. In fact, we had tried both PDC and DTF, and all results indicated that the in-degree calculated from the column-wise normalized PDC matched the clinically diagnosed results best.

The betweenness centrality was also effective in identifying EZs, which is in line with the results of previous studies (Wilke et al., 2011; Varotto et al., 2012; Burns et al., 2014). Unlike the in-degree (or out-degree), which quantify the information that flows into (or out from) a node, the betweenness centrality quantifies information transition through a node. Although the degree measures and the betweenness centrality capture different aspects of connectivity, our results showed that the in-degree and betweenness centrality identified almost the same EZs. These results indicated that most of the connections in the network ended within the identified EZs (flow into) and also passed through the other identified EZs. So the EZs are not isolated, and interactions between the nodes may lead to the onset of seizures. Thus, for complicated cases in which there may be several EZs and in which these are difficult to identify, combining the betweenness centrality with the in-degree will provide more information about the network and could be useful in identifying the EZs.

The successful estimation of Granger causality (TVPDC in our case) depends on the success of the TVMVAR model, since all the necessary information is derived from the model coefficients. In practice, the model order is an important issue. If the model order is too low, the model will not capture the essential dynamics of the dataset. If the model order is too high, it will also capture unwanted components, leading to overfitting and instability. For the traditional MVAR model, the model order can be optimized by criteria like the Akaike information criterion (AIC) (Akaike, 1974) or the Bayesian information criterion (BIC) (Schwarz, 1978), both of which take the measurement noise covariance matrix into account. However, for non-stationary time series, it may be not accurate to estimate the model order by applying these criteria (Havlicek et al., 2010). In this study, we tried different model orders, including 3, 5, 7, 9, and 11. We found that if the model order was too small (i.e., 3 or 5), the results had very little similarity to the clinical diagnosis. For model orders larger than 7, we obtained the same results as for 7. So we set the TVMVAR model order as 7, similarly to a study with similar conditions and settings (van Mierlo et al., 2011a, 2013). This setting gave us acceptable computational requirements and sufficient precision in the results.

In this study, the adaptation process of the Kalman filter was first evaluated in a quantitative manner. As far as we know, this is the first report of such an evaluation in any of the relevant studies. With this evaluation, the adaptation time of the Kalman filter could also be determined and the update coefficients (λ and Λ) could be properly selected. Generally, the update coefficients determine the adaptation speed. The larger the update coefficients, the more rapidly the Kalman filter follows changes in the dataset. However, a larger update coefficient may lead the filter to become unstable. van Mierlo et al. (2011a) chose the update coefficient as 10⁻³ based on knowledge gained from a previous simulation study and discarded the coefficient matrices of the TVMVAR model from the first 5 s during the adaptation. In our study, the optimal update coefficients were the same as theirs, but the adaptation time was quite different (62 ± 6.68 s, mean ± SD, n = 7). Actually, the states estimated during the adaptation process were influenced by the selected start time of the Kalman filter, and the results are incorrect and should be discarded. To precisely identify EZs, it is necessary to ensure correct construction of the network around the seizure onset time. Thus, the adaptation time should be precisely estimated, and for this the evaluation of the adaptation of the Kalman filter is indispensable.

During identification of the EZs, we calculated the mean values of the graph indices within a 15 s window around the seizure onset time (10 s before and 5 s after). Other studies used different windows, either crossing the seizure onset time with a few seconds before and after (Guye et al., 2006; van Mierlo et al., 2011b) or starting from the seizure onset time with a length ranging from a few seconds to 20 s (Wilke et al., 2010; David et al., 2011; Lu et al., 2012). The underlying network characteristics are believed to change before seizure onset. Although the exact time point at which this change occurs is not clear, the network should contain valuable information on seizure generation within a few seconds before onset. On the other hand, the seizure may rapidly propagate from the seizure onset zones to other brain regions within a few seconds or even more rapidly. As the EZs are the brain regions involved in the generation and early propagation of seizures, the window should be located within the seizure generation and early propagation phase. Thus, the time period during which the seizure has already become widespread has not been included in the window.

We chose the threshold used in identifying the EZs as 50% of the maximum (Wilke et al., 2009, 2011). This choice of threshold is acceptable since it filters out most of the outstanding regions with relatively high graph index values. In clinical settings, epileptologists identify EZs from information collected in multiple ways. So, for them, this rough selection of a threshold is sufficient, because all the information provided by this framework marks out the outstanding regions around seizure onset. Epileptologists may then comprehensively analyze the results together with other clinical indices to identify EZs.

Because the spatial sampling of the SEEG recordings is discrete and limited, if the SEEG electrodes do not cover the real EZs, the framework may produce erroneous results. Therefore, care should be taken when interpreting results calculated from SEEG recordings. Burns et al. (2014) found that the inadequate coverage of the EZs could lead to inconsistent brain states. We analyzed several seizures for each patient (data not shown) and obtained similar results. Therefore, it can be confirmed that the electrode coverage was adequate for the clinical cases we studied.

In summary, the applicability of an effective connectivity and graph theory analysis to localize EZs from SEEG recordings around seizure onset time was investigated in seven patients. We found that the EZs defined by our method corresponded well to the results of clinical diagnosis performed by epileptologists. Furthermore, among the three graph measures investigated, the in-degree is the best at identifying EZs. The betweenness centrality can reveal more information about the epileptic network and can also be used to identity EZs, but the out-degree seems to be unable to identify any of the EZs in our cases.

Author Contributions

JM, XY, YL, PZ, and JX conceived and designed the experiments. JM, XY, and JX performed the experiments. JM, YL, PZ, and PL analyzed the data. JM, YL, and PZ contributed materials and analysis tools. JM, PZ, YL, XY, PL, and JX wrote the paper.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by the Key Basic Research Project of Science and Technology Commission of Shanghai (13DJ1400303).

References

Akaike, H. (1974). A new look at the statistical model identification. IEEE Trans. Automat. Control 19, 716–723. doi: 10.1109/TAC.1974.1100705

CrossRef Full Text | Google Scholar

Arnold, M., Miltner, W. H. R., Witte, H., Bauer, R., and Braun, C. (1998). Adaptive ar modeling of nonstationary time series by means of kalman filtering. IEEE Trans. Biomed. Eng. 45, 553–562. doi: 10.1109/10.668741

PubMed Abstract | CrossRef Full Text | Google Scholar

Astolfi, L., Cincotti, F., Mattia, D., De Vico Fallani, F., Tocci, A., Colosimo, A., et al. (2008). Tracking the time-varying cortical connectivity patterns by adaptive multivariate estimators. IEEE Trans. Biomed. Eng. 55, 902–913. doi: 10.1109/TBME.2007.905419

PubMed Abstract | CrossRef Full Text | Google Scholar

Astolfi, L., Cincotti, F., Mattia, D., Marciani, M. G., Baccalà, L. A., de Vico Fallani, F., et al. (2006). Assessing cortical functional connectivity by partial directed coherence: simulations and application to real data. IEEE Trans. Biomed. Eng. 53, 1802–1812. doi: 10.1109/TBME.2006.873692

PubMed Abstract | CrossRef Full Text | Google Scholar

Baccalá, L. A., and Sameshima, K. (2001). Partial directed coherence: a new concept in neural structure determination. Biol. Cybern. 84, 463–474. doi: 10.1007/PL00007990

PubMed Abstract | CrossRef Full Text | Google Scholar

Bancaud, J., and Talairach, J. (1973). Methodology of stereo EEG exploration and surgical intervention in epilepsy. Revue d'oto Neuroophtalmologie 45, 315–328.

PubMed Abstract | Google Scholar

Basu, I., Kudela, P., Korzeniewska, A., Franaszczuk, P. J., and Anderson, W. S. (2015). A study of the dynamics of seizure propagation across micro domains in the vicinity of the seizure onset zone. J. Neural Eng. 12:046016. doi: 10.1088/1741-2560/12/4/046016

PubMed Abstract | CrossRef Full Text | Google Scholar

Beleza, P. (2009). Refractory epilepsy: a clinically oriented review. Eur. Neurol. 62, 65–71. doi: 10.1159/000222775

PubMed Abstract | CrossRef Full Text | Google Scholar

Berg, A. T., Berkovic, S. F., Brodie, M. J., Buchhalter, J., Cross, J. H., van Emde Boas, W., et al. (2010). Revised terminology and concepts for organization of seizures and epilepsies: report of the ilae commission on classification and terminology, 2005-2009. Epilepsia 51, 676–685. doi: 10.1111/j.1528-1167.2010.02522.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Boccaletti, S., Latora, V., Moreno, Y., Chavez, M., and Hwang, D. U. (2006). Complex networks: structure and dynamics. Phys. Rep. Rev. Section Phys. Lett. 424, 175–308. doi: 10.1016/j.physrep.2005.10.009

CrossRef Full Text | Google Scholar

Burns, S. P., Santaniello, S., Yaffe, R. B., Jouny, C. C., Crone, N. E., Bergey, G. K., et al. (2014). Network dynamics of the brain and influence of the epileptic seizure onset zone. Proc. Natl. Acad. Sci. U.S.A. 111, E5321–E5330. doi: 10.1073/pnas.1401752111

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiang, S., and Haneef, Z. (2014). Graph theory findings in the pathophysiology of temporal lobe epilepsy. Clin. Neurophysiol. 125, 1295–1305. doi: 10.1016/j.clinph.2014.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Coito, A., Plomp, G., Genetti, M., Abela, E., Wiest, R., Seeck, M., et al. (2015). Dynamic directed interictal connectivity in left and right temporal lobe epilepsy. Epilepsia 56, 207–217. doi: 10.1111/epi.12904

PubMed Abstract | CrossRef Full Text | Google Scholar

Cossu, M., Cardinale, F., Castana, L., Citterio, A., Francione, S., Tassi, L., et al. (2005). Stereoelectroencephalography in the presurgical evaluation of focal epilepsy: a retrospective analysis of 215 procedures. Neurosurgery 57, 706–718. discussion: 706–718. doi: 10.1227/01.NEU.0000176656.33523.1e

PubMed Abstract | CrossRef Full Text | Google Scholar

David, O., Blauwblomme, T., Job, A. S., Chabardès, S., Hoffmann, D., Minotti, L., et al. (2011). Imaging the seizure onset zone with stereo-electroencephalography. Brain 134, 2898–2911. doi: 10.1093/brain/awr238

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, L., Worrell, G. A., Lagerlund, T. D., and He, B. (2007). Ictal source analysis: localization and imaging of causal interactions in humans. Neuroimage 34, 575–586. doi: 10.1016/j.neuroimage.2006.09.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Florin, E., Gross, J., Pfeifer, J., Fink, G. R., and Timmermann, L. (2010). The effect of filtering on granger causality based multivariate causality measures. Neuroimage 50, 577–588. doi: 10.1016/j.neuroimage.2009.12.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Franaszczuk, P. J., and Bergey, G. K. (1998). Application of the directed transfer function method to mesial and lateral onset temporal lobe seizures. Brain Topogr. 11, 13–21. doi: 10.1023/A:1022262318579

PubMed Abstract | CrossRef Full Text | Google Scholar

Franaszczuk, P. J., Bergey, G. K., and Kaminski, M. J. (1994). Analysis of mesial temporal seizure onset and propagation using the directed transfer function method. Electroencephalogr. Clin. Neurophysiol. 91, 413–427. doi: 10.1016/0013-4694(94)90163-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J., Frith, C. D., and Frackowiak, R. S. J. (1993). Time-dependent changes in effective connectivity measured with pet. Hum. Brain Mapp. 1, 69–79. doi: 10.1002/hbm.460010108

CrossRef Full Text | Google Scholar

Granger, C. W. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica 37, 424–438. doi: 10.2307/1912791

CrossRef Full Text | Google Scholar

Guye, M., Régis, J., Tamura, M., Wendling, F., McGonigal, A., Chauvel, P., et al. (2006). The role of corticothalamic coupling in human temporal lobe epilepsy. Brain 129, 1917–1928. doi: 10.1093/brain/awl151

PubMed Abstract | CrossRef Full Text | Google Scholar

Haneef, Z., and Chiang, S. (2014). Clinical correlates of graph theory findings in temporal lobe epilepsy. Seizure Eur. J. Epilepsy 23, 809–818. doi: 10.1016/j.seizure.2014.07.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Hauser, W. A., Annegers, J. F., and Rocca, W. A. (1996). Descriptive epidemiology of epilepsy: contributions of population-based studies from rochester, minnesota. Mayo Clin. Proc. 71, 576–586. doi: 10.4065/71.6.576

PubMed Abstract | CrossRef Full Text | Google Scholar

Havlicek, M., Jan, J., Brazdil, M., and Calhoun, V. D. (2010). Dynamic granger causality based on kalman filter for evaluation of functional network connectivity in fMRI data. Neuroimage 53, 65–77. doi: 10.1016/j.neuroimage.2010.05.063

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaminski, M. J., and Blinowska, K. J. (1991). A new method of the description of the information flow in the brain structures. Biol. Cybern. 65, 203–210. doi: 10.1007/BF00198091

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J. Y., Kang, H. C., Cho, J. H., Lee, J. H., Kim, H. D., and Im, C. H. (2014). Combined use of multiple computational intracranial eeg analysis techniques for the localization of epileptogenic zones in lennox-gastaut syndrome. Clin. EEG Neurosci. 45, 169–178. doi: 10.1177/1550059413495393

PubMed Abstract | CrossRef Full Text | Google Scholar

Leistritz, L., Pester, B., Doering, A., Schiecke, K., Babiloni, F., Astolfi, L., et al. (2013). Time-variant partial directed coherence for analysing connectivity: a methodological study. Philos. Trans. R. Soc. Math. Phys. Eng. Sci. 371:20110616. doi: 10.1098/rsta.2011.0616

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Y. F., Yang, L., Worrell, G. A., and He, B. (2012). Seizure source imaging by means of fine spatio-temporal dipole localization and directed transfer function in partial epilepsy patients. Clin. Neurophysiol. 123, 1275–1283. doi: 10.1016/j.clinph.2011.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Milde, T., Leistritz, L., Astolfi, L., Miltner, W. H. R., Weiss, T., Babiloni, F., et al. (2010). A new kalman filter approach for the estimation of high-dimensional time-variant multivariate ar models and its application in analysis of laser-evoked brain potentials. Neuroimage 50, 960–969. doi: 10.1016/j.neuroimage.2009.12.110

PubMed Abstract | CrossRef Full Text | Google Scholar

Panzica, F., Varotto, G., Rotondi, F., Spreafico, R., and Franceschetti, S. (2013). Identification of the epileptogenic zone from stereo-eeg signals: a connectivity-graph theory approach. Front. Neurol. 4:175. doi: 10.3389/fneur.2013.00175

PubMed Abstract | CrossRef Full Text | Google Scholar

Plomp, G., Quairiaux, C., Michel, C. M., and Astolfi, L. (2014). The physiological plausibility of time-varying granger-causal modeling: normalization and weighting by spectral power. Neuroimage 97, 206–216. doi: 10.1016/j.neuroimage.2014.04.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenow, F., and Luders, H. (2001). Presurgical evaluation of epilepsy. Brain 124, 1683–1700. doi: 10.1093/brain/124.9.1683

PubMed Abstract | CrossRef Full Text | Google Scholar

Schlogl, A., Roberts, S. J., and Pfurtscheller, G. (2000). “A criterion for adaptive autoregressive models,” in Engineering in Medicine and Biology Society, 2000. Proceedings of the 22nd Annual International Conference of the IEEE, Vol. 2 (Chicago, IL), 1581–1582. doi: 10.1109/iembs.2000.898046

CrossRef Full Text

Schwarz, G. (1978). Estimating the dimension of a model. Ann. Stat. 6, 461–464. doi: 10.1214/aos/1176344136

CrossRef Full Text | Google Scholar

Sethi, M., Pedersen, M., and Jackson, G. D. (2016). Polymicrogyric cortex may predispose to seizures via abnormal network topology: an fMRI connectomics study. Epilepsia 57, E64–E68. doi: 10.1111/epi.13304

PubMed Abstract | CrossRef Full Text | Google Scholar

Téllez-Zenteno, J. F., Dhar, R., and Wiebe, S. (2005). Long-term seizure outcomes following epilepsy surgery: a systematic review and meta-analysis. Brain 128, 1188–1198. doi: 10.1093/brain/awh449

PubMed Abstract | CrossRef Full Text | Google Scholar

van Mierlo, P., Carrette, E., Hallez, H., Raedt, R., Meurs, A., Vandenberghe, S., et al. (2013). Ictal-onset localization through connectivity analysis of intracranial EEG signals in patients with refractory epilepsy. Epilepsia 54, 1409–1418. doi: 10.1111/epi.12206

PubMed Abstract | CrossRef Full Text | Google Scholar

van Mierlo, P., Carrette, E., Hallez, H., Vonck, K., Van Roost, D., Boon, P., et al. (2011a). Accurate epileptogenic focus localization through time-variant functional connectivity analysis of intracranial electroencephalographic signals. Neuroimage 56, 1122–1133. doi: 10.1016/j.neuroimage.2011.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

van Mierlo, P., Carrette, E., Hallez, H., Vonck, K., Van Roost, D., Boon, P., et al. (2011b). “Epileptogenic focus localization through connectivity analysis of the intracranial EEG: a retrospective study in 2 patients,” in 5th International IEEE/EMBS Conference on Neural Engineering (Ner) (Cancún), 655–658. doi: 10.1109/ner.2011.5910633

CrossRef Full Text

Varotto, G., Tassi, L., Franceschetti, S., Spreafico, R., and Panzica, F. (2012). Epileptogenic networks of type II focal cortical dysplasia: a stereo-eeg study. Neuroimage 61, 591–598. doi: 10.1016/j.neuroimage.2012.03.090

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Hernandez, J. M., and Van Mieghem, P. (2008). Betweenness centrality in a weighted network. Phys. Rev. E 77:046105. doi: 10.1103/PhysRevE.77.046105

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilke, C., Ding, L., and He, B. (2008). Estimation of time-varying connectivity patterns through the use of an adaptive directed transfer function. IEEE Trans. Biomed. Eng. 55, 2557–2564. doi: 10.1109/TBME.2008.919885

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilke, C., van Drongelen, W., Kohrman, M., and He, B. (2009). Identification of epileptogenic foci from causal analysis of ecog interictal spike activity. Clin. Neurophysiol. 120, 1449–1456. doi: 10.1016/j.clinph.2009.04.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilke, C., van Drongelen, W., Kohrman, M., and He, B. (2010). Neocortical seizure foci localization by means of a directed transfer function method. Epilepsia 51, 564–572. doi: 10.1111/j.1528-1167.2009.02329.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilke, C., Worrell, G., and He, B. (2011). Graph analysis of epileptogenic networks in human partial epilepsy. Epilepsia 52, 84–93. doi: 10.1111/j.1528-1167.2010.02785.x

PubMed Abstract | CrossRef Full Text | Google Scholar