The SEEG dataset was obtained from 10 patients (5 males) at Sanbo Brain Hospital of Capital Medical University in Beijing. All were selected based on the following criteria: refractory focal epilepsy and neurosurgical resection rendering the patient seizure-free during at least 2 years. This study was carried out in accordance with the recommendations of the Ethics Committee of the Sanbo Brain Hospital of Capital Medical University with written informed consent from all subjects. The basic characteristics of these recruited patients are described in Table 1. The age is presented as a range to avoid identifiable patient information.

The locations of electrodes were decided on the findings of non-invasive investigations, and the implanting operation aided by an advanced robot ROSA was completed accurately. Each electrode had multiple contacts, each of which was regarded as a channel. In addition, the average monitoring period was 1 week, and the sampling frequency was set to 512 Hz.

For each patient, we randomly selected two episodes (patient No. 8 had only one) containing epileptic seizures. These episodes were taken from 60 s before the pre-ictal stage to 30–60 s after seizure onset. It is worth mentioning that the pre-ictal period (IP) was divided due to the presence of high-frequency oscillations, also referred to as rapid discharges (22). If patients did not have a definite pre-ictal stage, the time series would begin at 60 s before clinical seizure.

The Granger causality is a classical technique to evaluate the statistical interdependence of multiple simultaneous time series (23), the strength of which can be quantified using multivariate autoregressive model (MVAR). To track the properties of non-stationary SEEG, the model parameters should be adaptive rather than fixed. Consequently, time-variant autoregressive model (TVAR) can be represented by the following equation: (1)X(n)=∑i=1pAi(n)X(n−i)+E(n), where X(n) = [x₁(n), …, x_k(n)]T is a vector indicating the recordings of K channels at time n, A_i(n) is the K × K coefficient matrix for delay i relative to time n, E(n) is the prediction error assumed to be white noise, and p is the model order which can be estimated by Akaike information criteria (24).

In fact, the identification of the time-variant coefficient matrix is an ill-posed mathematical issue. Ordinary methods (e.g., least squares, maximum likelihood) are unable to solve it unless time windows are adopted, which may spoil the data. Kalman filtering algorithm based on the state space model is entirely different and can estimate the coefficient matrix at any time. The details of algorithm implementation were described in Arnold et al. (11). An example is presented in Figure 1, wherein the prediction error and the fluctuation of a coefficient are depicted.

At time n, the p matrices jointly determine the causality relationships between different channels, and ADTF is then developed to integrate these into the frequency domain (25) (2)A(f,n)=I−∑i=1pAi(n)e−u.2πfi, (3)ADTFij(f,n)=|Hij(f,n)|2, where u is the imaginary unit, A (f, n) is the Fourier transformation of the coefficient matrix, and H_ij (f, n) is the element of H (f, n), the inverse matrix of A (f, n).

Since both the frequency and time can determine the value of H_ij (f, n), it may be more convenient to investigate the evolution with time if the information is integrated in a specific frequency band. Meanwhile, the ADTF is generally normalized with respect to the incoming flow, which leads to the spectrum-weighted ADTF (swADTF) (26) (4)swADTFij(n)=∑f=f1f2|Hij(f,n)|2∑l=1K|Hjl(f,n)|2∑k=1K∑f′=f1f2(|Hik(f′,n)|2∑s=1K|Hks(f′,n)|2).

The swADTF is calculated in the frequency band (3, 45 Hz) to exclude the interference of background electrical activity and include the effect of low gamma rhythm simultaneously. Notably, the values of swADTF have a slight fluctuation even for stationary signals. Therefore, it is necessary and reasonable to take the average value every 0.5 s.

Effective connectivity is defined as the influence that one node exerts over another, In contrast to functional connectivity, it can describe the directionality of interactions and the path of information flow. The abnormal activity of the epileptogenic area can stimulate the connected neurons and cause sporadic or periodic spikes. The EZ often maintains robust causality relationships with neuronal populations. Therefore, the top K non-diagonal elements of the causality matrix are selected to determine the dynamic threshold th (n). The connection matrix L is defined as follows: (5){swADTFij(n)≥th(n) and i≠j→Lij(n)=1swADTFij(n)≤th(n) or i=j→Lij=0, where L_ij(n) = 1 indicates that there exists an edge from node j pointing to node i. The process of constructing dynamic effective connectivity is described in the left branch of Figure 2.

The out-degree is the number of outgoing edges from one node. A large out-degree means that this node dominated the network at that time. Its role is extremely similar to the impact of the EZ on other neural ensembles in the production of synchronous discharges. Therefore, the out-degree is selected as the principal indicator of the epileptic area.

In 2002, Wendling et al. proposed an exquisite single neuronal population model on the basis of previous research results (18, 27, 28). In this model, physiological insights have been considered as follows: (i) the non-uniform alteration of GABAergic inhibition in experimental epilepsy (29); (ii) the possible depression of GABA_fast circuit activity by GABA_slow inhibitory postsynaptic currents (30). This model includes three subsets of neurons, namely the main cells (e.g., pyramidal cells), the dendritic-projecting inhibitory interneurons (GABA_A,slow receptors) and the somatic-projecting inhibitory interneurons (GABA_A,fast receptors). Excitatory and inhibitory connections between subsets are depicted in Figure 3A, and their intrinsic mechanism is described in the form of block diagram (Figure 3B).

Neuronal populations are connected and cooperate to achieve specific functions. To account for this organization, a coupled neuronal population model is established (Figure 4), and the connection strength is represented by the following equation: (6){Wij(n)=γswADTFij(n);i≠jWij=0;i=j where γ is the amplification factor that neutralizes the impact of normalization. Its value depends on the size of the network, which generally takes 320–360.

According to the block diagram, the following set of 12 × K differential equations that describe the state of every population can be established: (7){y˙1i=y6iy˙6i=EXC⋅aS(y2i−y3i−y4i)−2ay6i−a2y1iy˙2i=y7iy˙7i=EXC⋅a[pi(t)+C2S(C1y1i)+∑j=1,j≠ikWijy11i]−2ay7i−a2y2iy˙3i=y8iy˙8i=SDI⋅bC4S(C3y1i)−2by8i−b2y3iy˙4i=y9iy˙9i=FSI⋅gC7S(C5y1i−C6y4i)−2gy8i−g2y4iy˙5i=y10iy˙10i=SDI⋅bs(C3y1i)−2by10i−b2y5iy˙11i=y12iy˙12i=EXC⋅adS(y2i−y3i−y4i)−2ady12i−ad2y11ii=1,2…K.

The parameters’ interpretation and standard values are listed in Table 2 (18, 28). In the process of simulation, it is assumed that all populations have the same parameters values so that only the relationship between seizures and networks is investigated.

It is also possible to perform virtual surgery and observe outcomes by removing the specified population from the model. For instance, if population u is selected for removal, the coupling matrix W(n) can be changed into W′(n) (8){Wij′(n)=Wij(n);i≠uandj≠uWij′(n)=0;i=uorj=u

Thus, it is easy to investigate the outcomes of different operation schemes, including the removal of the EZ based on theoretical calculation.

The comparison between the out-degree calculation results and clinical conclusions is shown in Table 3. The electrode contacts with the maximum out-degree in various stages are regarded as the theoretical epileptogenic focus. SEEG reports accomplished by several experienced electrophysiologists indicate which channels are involved in the onset of epileptiform discharges. All of the epileptogenic channels in the SEEG reports are located in the resection region. Fine surgical outcomes can ensure that the clinical conclusions are credible enough to evaluate the theoretical results. In addition, the pre-ictal stage is combined with the ictal process since it has long been recognized as the sign of seizure onset (22).

The accuracy is the proportion of true positives to all calculated results. The detection rate measures corresponding percentage in clinical conclusions. They describe the influence of false positives and false negatives, respectively. According to Table 3, SEEG reports list 34 epileptogenic channels, while the calculated results show 35, 29 of which are actual positives, meaning that the accuracy is 82.86%, and the detection rate is 85.29%. In addition, the results from eight patients demonstrate that the driving hubs (the channels with the maximum out-degree) of the interictal activity overlap with those in the ictal stage. This implies that the cause of seizures may lurk in the interictal activity.

Particularly, the out-degree of each channel with time is depicted for patient No. 2 in Figure 5A. The stages are divided according to the features of the SEEG data and clinical symptoms. Although the salient characteristics of the brain networks are relatively stable within a certain stage, they differ from one stage to another. To clearly reveal the transient process, a compromise is struck. Only the dots whose out-degree is above 30 are marked. Meanwhile, the networks in the same stage are summed and then averaged. The distribution of the mean out-degree is shown to interpret the difference of connectivity patterns (Figure 5B). In addition, the results are in good agreement during the two episodes.

It is shown that E04, the fourth contact of electrode E, governs the whole networks in both the interictal stage and the pre-IP. Similarly, J14 and M08 dominate activities in the ictal stage. The SEEG report verifies that intermittent spikes can be observed in E04 during the interictal period (IIP). Electrophysiologists deem that the abnormal discharges spread from E04 to electrode J and M. The calculated results are included within the resection area (Figures 5C,D). Fine surgical outcome also proves that this method is effective and accurate for localizing the EZ.

Spikes are one of the most common biomarkers for indicating abnormal brain activity. In the simulation process, sporadic spikes are seen in the interictal stage and rhythmic spikes are observed during the IP (Figure 6B), which are consistent with clinical findings (Figure 6A). Virtual surgical resections are then carried out, and the outcomes of different operations are observed (Figures 6C,D). These simulations in three conditions (no resection, random resection, and removing the epileptogenic channel) are repeated 30 times, respectively. Compared with random resection, the removal of the EZ is more effective to reduce the amount of spikes (at least 75%) (Figure 6E). It is worth mentioning that the clinical resection range is wider than single channel, so sporadic spikes can be further eliminated.

However, two cases of simulations are not in accordance with the real situation. In Figure 6F, periodic spikes occupy almost the entire IIP, while there are only partial spikes during the ictal process. Nevertheless, the channels with the maximum out-degree in the interictal stage were also the targets for surgical resection. These spikes can be eliminated after removing the calculated source. Such interictal networks are also pathological and may well be the cause of seizures.

Epilepsy is regarded as a network-level phenomenon (32). This type of network disorder may have already existed before clinical seizures. Intermittent spikes during this period can be regard as implications. Previous studies have demonstrated that interictal epileptic events were helpful for confirming the range of abnormal areas (33, 34). Wilke et al. concluded that the sources of interictal spike activity might be the generators underlying the clinical seizures (13). In Table 3, eight patients demonstrated that the sources of the interictal activity overlapped with the locations that the ictal activity originated from. In addition, some networks in the interictal stages exhibited epileptogenic behaviors and rhythmic spikes were observed in the simulations. This organization that is prone to producing spikes may push brain function to the edge of disorder. In other words, eliminating the interictal spikes or correcting these pathological networks should be part of epilepsy treatment.

It is increasingly understood that seizures arise from epileptogenic networks, through which epileptiform discharges are produced and propagated (35). However, it is still unclear that which type of network can be considered as epileptogenic connectivity. Epileptogenic networks are not equal to the networks in the ictal stage. Some networks in the interictal stage were also epileptogenic. Varotto et al. reached a similar conclusion on the basis of SEEG recordings from patients with type II focal cortical dysplasia (36). Brain functional connectivity identified by the fMRI also showed an abnormal pattern in the absence of epileptiform discharges (37).

We have known that the channel with the maximum out-degree varies from stages. However, the occurrence of seizures is not only due to the transfer of important nodes. Qualitative change of connectivity patterns must follow. Coupled neuronal population model was used to simulate the time course in different patterns. An increased number of spikes indicate that this connectivity pattern is more likely to generate ictal activity. The nature of epileptogenic network should be a specific organization which likely enhances local excitability and facilitates synchronization. Because the epileptogenic region affects other neural ensembles and enables them to synchronize, it should be responsible for the formation of these pathological networks. Virtual surgery also suggested that the removal of the EZ could significantly reduce the amount of spikes. Sinha et al. also concluded that the resection of these regions could reduce the overall likelihood of seizures after conducting simulations with a bi-stable model (21). Hence, it is helpful to explain why seizures are reduced or suppressed after accurate resection.

The distribution of the out-degree could aid identification of the EZ and deepen our understanding of the difference between connectivity patterns. However, it failed to describe the qualitative changes of network organization. Coupled neuronal population model was used to demonstrate the network behavior and judge whether the network was epileptogenic or not. It was also convenient to observe the outcomes of different operation schemes. There are two reasons to carry out virtual surgery. First, contradiction may appear when the EZ overlaps with the brain function area. Hypothetical surgery provides new options for controlling seizures with regard to preserving major function. Second, virtual surgery is beneficial in terms of balancing the resection range and therapeutic effect (38, 39). It is well known that a smaller resection region means lower physiological cost. By removing one or more channels, it is possible to assess the corresponding surgical effect, and the optimal scope can be determined. In a word, the treatment for refractory focal epilepsy should be a comprehensive and systematic project instead of just achieving seizure freedom.

This study consists of the EZ localization and coupled neuronal population model. The former aims at identifying the electrode contacts triggering the abnormal discharge. However, seizures are often associated with abnormal changes in brain synchronization mechanisms. This model is used to find epileptogenic connectivity patterns and evaluate the role of the EZ on pathological network. The results suggested that the presence of the EZ made the network more easily synchronized. Accordingly, if the EZ is removed, pathological network can be corrected and seizures are reduced or suppressed. It is helpful to understand the relationship between the EZ and epileptogenic networks.