# Localizing targets for neuromodulation in drug-resistant epilepsy using intracranial EEG and computational model

- Department of Biomedical Engineering, School of Electrical Engineering, Shenyang University of Technology, Shenyang, China

Neuromodulation has emerged as a promising technique for the treatment of epilepsy. The target for neuromodulation is critical for the effectiveness of seizure control. About 30% of patients with drug-resistant epilepsy (DRE) fail to achieve seizure freedom after surgical intervention. It is difficult to find effective brain targets for neuromodulation in these patients because brain regions are damaged during surgery. In this study, we propose a novel approach for localizing neuromodulatory targets, which uses intracranial EEG and multi-unit computational models to simulate the dynamic behavior of epileptic networks through external stimulation. First, we validate our method on a multivariate autoregressive model and compare nine different methods of constructing brain networks. Our results show that the directed transfer function with surrogate analysis achieves the best performance. Intracranial EEGs of 11 DRE patients are further analyzed. These patients all underwent surgery. In three seizure-free patients, the localized targets are concordant with the resected regions. For the eight patients without seizure-free outcome, the localized targets in three of them are outside the resected regions. Finally, we provide candidate targets for neuromodulation in these patients without seizure-free outcome based on virtual resected epileptic network. We demonstrate the ability of our approach to locate optimal targets for neuromodulation. We hope that our approach can provide a new tool for localizing patient-specific targets for neuromodulation therapy in DRE.

## 1 Introduction

Epilepsy is a neurological disease caused by disorder of the brain network (Terry et al., 2012; Lam et al., 2016). It has the characteristics of recurrent seizures, which often bring irreversible brain damage and affect the normal life of patients with epilepsy (Trinka et al., 2015). About 70% of patients can be cured by taking antiepileptic drugs, and 30% of them will develop drug-resistant epilepsy (DRE) (Kwan and Brodie, 2000; Stephen et al., 2006; Brodie et al., 2012). Patients with DRE can be treated with surgery (Choi et al., 2008) or neuromodulation (Schulze-Bonhage, 2017; Davis and Gaitanis, 2020; Sisterson and Kokkinos, 2020), such as transcranial magnetic stimulation (TMS) (Davis and Gaitanis, 2020), transcranial focused ultrasound (tFUS) (Lin et al., 2020; Zou et al., 2020). In neuromodulation therapy, different brain regions or nerves can be chosen as target, such as the vagus nerve (Stern et al., 2021), thalamus (Ryvlin and Jehi, 2021), hippocampus (Abouelleil et al., 2022), or localized epileptogenic zone (Tsuboyama et al., 2020; Rincon et al., 2021). Neurologists use intracranial EEG (iEEG), MRI and other methods combined with clinical experience to define the brain regions responsible for seizure generation and resect these regions to prevent seizure. About 30% of patients with DRE failed to achieve seizure freedom after surgical intervention (Janszky et al., 2005; de Tisi et al., 2011). Most of them are not suitable for further surgery because the suspected brain regions have been damaged. Neuromodulation is a promising technique for these non-seizure-free patients. However, to our best knowledge, the reports of localizing neuromodulatory targets for patients who failed to achieve seizure freedom are few.

EEG is widely used in the diagnosis of epilepsy. Comparing with scalp EEG, iEEG electrodes need to be embedded in the patient’s skull. The intracranial electrodes are closer to the epileptogenic zone (Kovac et al., 2017), which facilitates subsequent resection of the epileptogenic area. IEEG recording techniques include subdural grids, strips, and depth electrodes. For different epilepsy patients, different iEEG techniques need to be selected (Kovac et al., 2017).

IEEG recordings reflect the characteristics of epileptic networks and have the function of localizing epileptogenic tissues in epilepsy patients (van Diessen et al., 2013; Taylor et al., 2015; Sinha et al., 2017). Network methods can be used to extract the epileptic network, such as Pearson correlation, Granger causality (Coben and Mohammad-Rezazadeh, 2015; Sinha et al., 2017). The coefficients of Pearson correlation represent the correlation between variables. Sinha et al. (2017) calculated the coefficients between different EEG channels and used them as undirected connectivity of the epileptic model. Granger causality explores direct or indirect relationships between variables. Directed transfer function (DTF) computes interactions between input signals in frequency domain (Franaszczuk et al., 1985). A variety of network characteristics can be quantified on the extracted network matrix. One single network feature cannot fully explain all the properties of the network (van Diessen et al., 2013). Different features of network were used in epilepsy studies (Sinha et al., 2017, 2021; Paldino et al., 2019). Seizure is a dynamic process. It is difficult to explore dynamical behaviors of the original signals based on the extracted network matrix. The quantified features of epileptic networks cannot comprehensively describe the dynamics of seizure onset and termination.

Neural computational model can be used to better simulate the dynamical process of seizure (Proix et al., 2018; Saggio et al., 2020; Sip et al., 2022). Richardson proposed a method of combining dynamics and connectomics to explain the abnormal dynamics of epileptic networks (Richardson, 2012). Lopes da Silva et al. (2003) modeled transition states between normal and epileptic states for predicting epileptic pathways. Creaser et al. (2002) modeled the node dynamics and the coupling relationship between nodes, and obtained the transient dynamics during epileptic seizures. Numerous models have been used to explain the physiology of epilepsy or epileptic activity (Wendling et al., 2016). In computational models, the state of the system is commonly changed by adjusting either excitatory or inhibitory parameters, such as Z6 model (Benjamin et al., 2012) and Epileptor model (Proix et al., 2018). In the Z6 model, the values of different excitatory parameters determine whether the system is in normal or epileptic state.

In this paper, we propose a novel approach for localizing targets for neuromodulation in patients with DRE, especially for patients without achieving seizure freedom after surgery. The patient-specific epileptic network is reconstructed using multi-unit computational model. The most effective node for neuromodulation in preventing seizure is localized by introducing external stimulation. The effectiveness of our proposed approach is validated on a multi-variate autoregressive (MVAR) model. Then we validate the approach on a group of DRE patients with iEEG recordings. Finally, the candidate targets for neuromodulation are provided using the proposed approach and virtual resected network of those DRE patients.

## 2 Methods

### 2.1 Data and subject description

IEEG recordings from 11 patients were analyzed in this study. The datasets were obtained from the IEEG public website (http://www.ieeg.org), and all patients with DRE had received surgical treatment. Three patients are in seizure-free group with good outcome, who were scored as international league against epilepsy (ILAE) 1 (completely seizure-free) or 2 (no seizures, only auras) (Wieser et al., 2001). The other eight patients are in non-seizure-free group with poor outcome, who were scored as ILAE 3–6 (non-seizure-free). Interictal iEEG recordings of 10 min duration are chosen several hours away from any seizure. The iEEG data are divided into segments of 1 s duration. Each segment overlaps the previous one by 0.5 s. The sampling rate of recordings is 500 Hz. We evaluate the overlapping between the localized target nodes and the resected regions.

### 2.2 Directed transfer function with surrogate analysis

The directed transfer function (DTF) is a multi-channel directional measurement method based on Granger causality and autoregressive models (Franaszczuk et al., 1985). This method calculates the causal connection matrix between multi-channel EEG signals and measures the causal relationship between channels. The multichannel EEG process in the framework of autoregressive model (AR) can be described by the following equation (Franaszczuk et al., 1985; Kaminski and Blinowska, 1991):

where *x*(*t*) = [*x*_{1}(*t*), *x*_{2}(*t*), …, *x*_{N}(*t*)] is the vector of EEG *N*-channel process, *p* is the order of the model. *A*_{0} is identity matrix, *A*_{1}, *A*_{2}, … , *A*_{p} are the *N* × *N* matrices of model coefficients, *w*(*t*) = [*w*_{1}(*t*), *w*_{2}(*t*), …, *w*_{N}(*t*)] is the vector of multivariate zero mean uncorrelated white noise process. We use the order selection criteria of Akaike’s Final Prediction Error (FPE) criterion implemented in ARFIT toolbox (Akaike, 1971; Schneider and Neumaier, 2001).

The coefficients *A*_{j} can be obtained from (1) by multiplying its both sides by *x*^{T} is transposed vector of *x*. We get following equation (Kaminski and Blinowska, 1991):

where *s* for the vector *x*, *E* means expectation value. Applying the z-transform to the both sides of Eq. 1 (Franaszczuk et al., 1985), we have

where *H*(*z*) is the transfer function. Set *z*^{−1} = *e*^{−i2πfΔt}, where *f* is the frequency, Δ*t* is the sampling interval. Then we get *H*(*f*), where *H*_{ij}(*f*) is the directed causal relationship from the node *j* to the node *i*.

The directional characteristic of the information flow from node *j* to node *i* is defined as following:

Note that the value of

Surrogate data is a statistical method of analyzing nonlinear signals that facilitates the analysis of EEG signals (Dolan and Spano, 2001). We generate surrogate signals by assigning the phase of the EEG signal randomly in 200 times. The strongest 5% of the total possible causal connection are kept for further analysis. The network characteristics in high frequency gamma band are most closely correlated with improved postsurgical outcome (Wilke et al., 2011). Our preliminary study on seizure-free group also showed similar results. In this study, the network analysis focuses on gamma rhythm (31–80 Hz).

### 2.3 Other methods to build brain network

Besides DTF, there are other ways to build brain networks. The Pearson correlation coefficient (PCC) reflects the linear correlation between iEEG channels. We divide the iEEG into 1 s data segments. The PCC calculates the degree of linear correlation between two variables. It is the ratio of the covariance and standard deviation between two signals, as shown in the following:

where *P*_{a,b} represents the degree of linear correlation between *n* dimensional signal *a* and *b*.

Partial directed coherence (PDC) analyzes the connectivity between multi-channel signals, which is also based on Granger causality. The calculation method of *p* and *A*_{j} is the same as that of DTF. The transfer function

where *I* is the identity matrix.

Isolated effective coherence (iCoh) is similar to PDC. This method computes the interrelationships of directly related nodes, but zeros out all other indirect causal relationships (Pascual-Marqui et al., 2014). When we compute the causal relationship from node *j* to node *i*, other nodes except node *j* and node *i* are called irrelevant nodes. Node *j* is the relevant nodes of node *i*.

The weighted phase lag index (wPLI) measures the phase correlation between signals by weighting the cross-spectrum of the phases of the two signals. First, we need to calculate the phase lag index (PLI) between the signals (Li et al., 2021):

where sgn represents the sign function, |⋅| is absolute value, Δ*θ*(t) represents the instantaneous phase difference between the input signal *s*_{1} and the output signal *s*_{2}. Then, wPLI is calculated to quantify the phase agreement between the signals (Li et al., 2021),

where *A*_{1} and *A*_{2} are the corresponding amplitudes of the *s*_{1} and *s*_{2}, respectively.

Relative entropy is an asymmetric measure (Kullback and Leibler, 1951), also known as Kullback-Leibler divergence (KLDIV), quantifies the difference between two signals. *x*(*t*) is divided by 1 s and overlapped by 50%, and then subjected to short-time Fourier transform to obtain *X* (*n*, *f*). The normalized spectrogram is shown in Eq. 9.

Suppose *W*_{y}(*n*, *f*) is the normalized spectrum of the signal *y*(*t*), the KLDIV from *W*_{x} (*n*, *f*) to *W*_{y}(*n*, *f*) is as follows:

KLDIV is asymmetric. When its value is higher, the difference between the signals is larger.

### 2.4 MVAR model

The electrophysiological activity in short duration can be viewed as a MVAR process. In this study, the MVAR model is written in the following differential form (Baccala and Sameshima, 2001),

where *X*_{i}, *i* = 1, 2, … , 5, represents the *i*th node of the network, and *w*_{i}, *i* = 1, 2, … , 5, is the white noise. The coefficients between nodes represent the causual relationships of different nodes.

In order to find the method with best performance for constructing epileptic network and validate the effectiveness of our proposed approach, we use a five-node MVAR model to simulate the causal relationship between nodes. The time course of activity assigned on each node is generated by model (11). Simplified MVAR model could be used to simulate epileptic sources (Hosseini et al., 2018). There are both unidirectional and bidirectional connections in model (11). Node *X*_{1} is simulated as epileptic node where seizure starts. Node *X*_{2}, *X*_{3} and *X*_{4} are neighbor nodes., and node *X*_{4} has bidirectional connection with node *X*_{5}, which represents remote normal tissue.

There are different types of measures to construct the brain network. We choose several commonly used network measures based on causality, coherence, or information theory. We also construct the causal network combining with surrogate analysis to determine which method matches the original network best. We compare nine methods for constructing causal network, including PCC, DTF, DTF with surrogate analysis (DTF-SA), PDC, PDC with surrogate analysis (PDC-SA), iCoh, iCoh with surrogate analysis (iCoh-SA), wPLI and KLDIV. We compute the correlation coefficients between the extracted connectivity matrix and the ground truth of the model (11), which are then normalized by the maximum value. According to the value of the correlation coefficient, we choose the method with the best performance to construct brain network.

### 2.5 Multi-unit computational model

The Z6 model can simulate the dynamic process of the interaction between nodes due to information transmission during epileptic seizures, and intuitively describe the state transition of nodes. This model contains a fixed point and a limit-cycle. The noise system controls one of two factors in order to control the trajectory of the system. In the noisy system, the deterministic part at the drift coefficient can be expressed by the following single complex equation (Benjamin et al., 2012):

where *z* is a complex parameter, *z* = *x* + *iy*. *a* and *b* are real numbers (*a* = −1, *b* = 2), *ω* controls the oscillation frequency of the system, *λ* is the possible attractors of the system. The parameter *λ* determines the state of the system, and we choose 0 < *λ* < 1. We can consider *λ* to be the excitability parameter of the system. When *λ* approaches 1, the system is more excitable (Benjamin et al., 2012).

The nodes of the brain network have the characteristic of bidirectional functional connectivity, forming a network of interconnected nodes. We extend the equation to a network model with *N* nodes:

where *G*_{ij} is the normalized information connectivity matrix between nodes, *w*(*t*) represents white noise with a mean of 0.0003 and a standard deviation of 0.05 (Sinha et al., 2017), *α* is the coefficient of noise. *β* equals 0.02 here. In order to achieve the same order of magnitude as the undirected symmetric information-connected matrix (Benjamin et al., 2012), the normalized matrix *G* is then multiplied by *K* = 1,000. The connectivity matrix *G* describes the topology of the epileptic network, and determines the interaction between each node of the system. When iEEG dataset is analyzed, the patient’s epileptic networks are constructed by network measures, which is then used as the matrix G of the computational model. In this case, the dynamic behavior of the model is determined by patient’s specific brain network. The model was solved numerically using a fixed step Euler-Maruyama solver with a step size of 0.05.

During the numeric simulation, the network is driven by random noise. The time for a node to change from a stationary state to an oscillating state is called the escape time (*T*_{es}) (Sinha et al., 2017). *T*_{es} is used as an indicator for predicting seizure (Benjamin et al., 2012). We use the Z6 model to simulate the epileptic brain as a bi-stability state network (Goodfellow and Glendinning, 2013; Sinha et al., 2017). The probability of a node entering epileptic state is inversely proportional to *T*_{es} (Petkov et al., 2014; Sinha et al., 2017). It is also proportional to the stability of the system, and the value of *T*_{es} decreases as the parameter *λ* increases (Benjamin et al., 2012). The parameters *λ* of all nodes are set to *λ*_{0} or *λ*_{1} on non-target or target nodes, respectively. The optimal *λ*_{0} and *λ*_{1} are chosen by grid search of *λ*_{0} − *λ*_{1} pairs. *λ*_{0} is chosen between 0 and 0.5 in step of 0.05, and *λ*_{1} is chosen between 0.5 and 1 in step of 0.05. The Δ*T* is calculated at each *λ*_{0} − *λ*_{1} pair. The optimal *λ*_{0} and *λ*_{1} are found when there is the largest standard deviation of Δ*T* for all pairs. First, we set the *λ* of all nodes to *λ*_{0}, and record the *T*_{es} as *T*_{0}. Then, we change the *λ* of each target node to *λ*_{1}, and record the *T*_{es} as *T*_{1}. The difference between *T*_{0} and *T*_{1} is the change in escape time (Δ*T*),

### 2.6 Localizing targets for neuromodulation

The procedure of localizing targets for neuromodulation is shown in Figure 1. First, the segmented iEEG recordings are used, as shown in Figure 1A. The data are processed to construct patient-specific epileptic network (Figure 1B). The epileptic network is in the form of causal connectivity matrix. The multi-unit neural computational model is constructed based on the epileptic network and the Z6 model (Figure 1C). The number of nodes of the multi-unit model is same as the number of channels in iEEG the recordings. The optimal values of *λ* are determined (Figure 1D). Using the selected parameters, the Δ*T* of each node is then calculated (Figure 1E). The distribution of Δ*T* is plotted on patient’s head model, and the node with the largest value of Δ*T* is selected as optimal target. The optimal target is compared with the resected regions of epilepsy patient, as shown in Figure 1F.

**FIGURE 1**. Procedure of localizing target for neuromodulation in drug-resistant epilepsy. **(A)** The segmented iEEG recordings. **(B)** Patient-specific epileptic network based on iEEG data. **(C)** Multi-unit neural computational model based on the epileptic network and the Z6 model. **(D)** Parameter optimization for *λ*. **(E)** Calculating the change of escape time (Δ*T*) of each node. **(F)** Localizing optimal target with the largest Δ*T* value.

### 2.7 Validation of neuromodulation

We calculated *T*_{es} of epileptic brain networks in all patients with inhibitory modulation on the localized target nodes and non-target nodes. The Wilcoxon rank sum test was used (Wilcoxon, 1945), and *p* < 0.01 was chosen as significance threshold. The proposed approach was then performed to localize candidate targets for neuromodulation in patients without seizure-free outcome.

## 3 Result

### 3.1 Localizing the critical node of MVAR model

The model (11) is shown in Figure 2A. The node *X*_{1} is the main driven force for the model to enter oscillatory state, which is simulated as the epileptogenic node. The normalized correlation coefficients of nine different methods are plotted in Figure 2B. Based on the extracted causal connectivity by DTF-SA method, we localize the target for neuromodulation using our proposed approach. The optimal values of *λ*_{1} and *λ*_{0} is 0.85 and 0.50, respectively. The result of Δ*T* is shown in Figure 2C, and the node *X*_{1} is with the highest value. The effectiveness of external modulation on the node *X*_{1} is shown in Figure 2D. The *T*_{es} of the network decreases significantly while the external excitatory stimuli is applied on node *X*_{1}.

**FIGURE 2**. Validating the proposed approach on the five-node model. **(A)** Five-node causal network. **(B)** Normalized correlation coefficient between the constructed network and the ground truth. **(C)** Δ*T* for each node. **(D)** The *T*_{es} value of the network with and without external modulation.

### 3.2 Localizing targets for neuromodulation using patient data

The localized targets for neuromodulation of 11 patients are listed in Table 1. For patient P1-P3 with seizure-free outcome, the targets are inside the resected regions. In non-seizure-free group, localized targets for patients P4, P10, and P11 are outside the resected regions, and localized targets for patients P5-P9 are inside the surgical resected regions. The mean error distance is 8.4 mm in non-seizure-free group. The values of parameter *λ*_{0} for 11 patients is 0.22 ± 0.08, *λ*_{1} is 0.99 ± 0.01 (mean ± SD).

The localized targets for neuromodulation of patient P3 and P4 are plotted in Figures 3A,C, respectively. Red region indicates the stimulation is effective to suppress seizure in Figures 3A,C. Patient P3 belongs to the seizure-free group. The resected region of this patient was mainly in left lateral frontal cortex. The electrode LFG48 with the highest Δ*T* is selected as the target, and reside in the resected region (red rectangle), as shown in Figure 3A. Patient P4 belongs to the non-seizure-free group. The resected region was mainly in left parietal cortex. The electrode LPG7 with the highest Δ*T* is selected as the target, which is 28.0 mm away from the resected region (red rectangle), as shown in Figure 3C. The distribution of Δ*T* for choosing optimal *λ*_{0}−*λ*_{1} pair is plotted in Figures 3B,D. The *λ*_{0} = 0.30 and *λ*_{1} = 1.00 for both patients.

**FIGURE 3**. The localized targets for neuromodulation of patient P3 and P4. **(A)** The distribution of Δ*T* for patient P3. The white circles represent the contacts of the electrodes. Red and blue indicate the magnitude of the value of Δ*T*. The red color indicates that the value of Δ*T* is large. The blue color indicates that the value of Δ*T* is small. **(B)** The optimal value of *λ* for patient P3. *λ*_{0} = 0.30 and *λ*_{1} = 1.00. The red color indicates that the value of standard deviation (S. D.) of Δ*T* is large. **(C)** The distribution of Δ*T* for patient P4. **(D)** The optimal value of *λ* for patient P4. *λ*_{0} = 0.30 and *λ*_{1} = 1.00.

The Δ*T* distributions for the other 9 patients are plotted in Figure 4. Among them, patients P1 and P2 belonged to the seizure-free group. The surgical field in both patients was in the left temporal lobe. Targets identified by our method that were effective in eliminating epilepsy were located at the surgical field of each of the two patients.

**FIGURE 4**. The distribution of Δ*T* for the other nine patients except P3 and P4. Red region indicates that the effect of modulation is strong.

### 3.3 The effectiveness of neuromodulation

The *T*_{es} is normalized relative to its minimum and maximum values of each patient. The horizontal red bar represents the median in Figure 5. The median of the *T*_{es} is 0.21 and 0.06 for target and non-target nodes, respectively. The values of *T*_{es} on the localized target nodes are significantly longer than the values on non-target nodes (*p* < 0.01). The modulation on the localized target nodes is more effective in suppressing seizure than non-target ones. Outlier data indicate that modulation on some non-target nodes is more (or less) efficient than modulation on other nodes.

**FIGURE 5**. The normalized *T*_{es} when external modulation is applied on the localized target and non-target nodes. The horizontal red bar represents the median, and the blue box represents the rang from first to third quartile, respectively. The horizontal black lines represent the upper and lower limits, and the red plus sign represents outlier data.

We remove the iEEG channels in the resected regions, and reconstruct virtual resected networks for the non-seizure-free patients. In the virtual resected network, the iEEG channels in resected regions are removed. The connectivity matrix G is calculated using other channels, which results in a smaller matrix. The value of *λ*_{0} is set to 0.22 for the non-target of the resected network, and *λ*_{1} = 1. The localized candidate targets for neuromodulation are listed in Table 1. The mean distance between the new targets and the surgical resected regions is 26.9 mm. The localized targets for Patient P4, P10, and P11 are not changed before and after surgical resection. The mean distance between the new localized targets and the resected regions is 29.6 mm for patient P5-P9.

## 4 Discussion

Localizing the effective targets is the key to neuromodulation therapy. The proposed approach identified the node *X*_{1} of the MVAR model as optimal target successfully, which is the designed node to drive the model to oscillatory state. The network measures based on correlation, causal effects, phase lag, and information entropy were compared for reconstructing the network. The DTF-SA method showed highest similarity between the reconstructed network and the ground true. Our results could help other study for choosing network measures. We choose the Z6 model as the network node because of the relative low computational cost. Other neural computation models could also be adopted, such as Epileptor model (Proix et al., 2018). The parameter *λ*_{0} of non-target node is the only parameter need to be determined for the network except for connectivity matrix. Our result show that the mean value of *λ*_{0} for non-target nodes is good in most cases when analyzing the iEEG dataset, and the value *λ*_{1} for target node is set to 1, which leads that our proposed approach is easy to use.

Patient P1, P2, and P3 have undergone surgical resection, and achieved good outcome. The epileptogenic tissues are assumed inside the resected regions. Based on the patient-specific epileptic network, all localized targets for those patients reside in the resected regions. Those results indicate that our approach find the target responsible for seizure generation. The localized targets for patient P5-P9 also reside in the resected regions, which is consistent with judgement of neurologist. The epilepsy is a brain network disease. Resection of brain tissue changes the topology of epileptic network, and seizure may start from other brain location. Applying the external stimulation on a given node will result in influence on whole network. In this context, we believe the neuromodulation measure may lead to better outcomes for those patients without seizure-free outcome. The localized targets for patient P4, P10, and P11 are outside the resected regions based on the epileptic network before surgery, and the mean distance from the surgical resected regions is 22.4 mm. After removing the resected nodes, our approach also localizes the targets on the same electrodes. This result reflects that our proposed network is stable on localizing the targets even with the virtual resection.

When applying inhibitory stimulation on the localized target nodes, it can significantly delay the brain network from entering a state of oscillation relative to the non-target nodes, as shown in Figure 5. This demonstrates the effectiveness of the proposed method for preventing epileptic seizures. Neuromodulation is used as a non-destructive means of brain network regulation, such as TMS, tFUS, which have different spatial resolutions (Davis and Gaitanis, 2020; Lin et al., 2020; Zou et al., 2020). Considering the applicability of our proposed method, the number of targets selected for neuromodulation is 1. Selecting 2 or more targets for neuromodulation will have a more obvious modulation effect, but it is not suitable for neuromodulation methods with low spatial resolution, such as TMS.

Our method induces resting-state brain networks into epilepsy *via* stimulation parameters. This approach differs from current neuromodulation treatments. For example, tFUS suppresses seizure by reducing the excitability of the nervous system (Folloni et al., 2019; Lin et al., 2020; Zou et al., 2020). However, in clinical surgery, the traditional method is to find the epilepsy surgery area by evoking electrical stimulation. Combined with the clinical surgical process, we select parameters that can induce the brain network to enter the epileptic state to determine the target of neuromodulation. This choice ensures the practicality of our method.

Furthermore, this approach has limitations. The iEEG recording is an invasive measure mainly for presurgical evaluation of DRE patient. We have not applied this measure on scalp EEG recording. The low coverage of the intracranial electrodes on the epileptogenic zone may result in the inaccurate of constructing epileptic network, which ultimately leads to poor localization of the targets. Validating our method in the real application will further advance the technology. On the one hand, the Z6 model is a noise-driven computational model that simulates epileptic seizures (Sinha et al., 2017). Due to the randomness and uncertainty of noise, multiple calculations are required to avoid accidental factors. The long computational time is another limitation of our method. On the other hand, we chose the Z6 model to simulate the dynamic process of epileptic seizures. Other neural computational models can also replace the Z6 model, such as the Epileptor model. Therefore it is necessary to choose different and more suitable parameters to simulate the process of neuromodulation.

## 5 Conclusion

The effective neuromodulation therapy is very important for DRE patients with bad surgical outcome. The DTF with surrogate analysis is more suitable for constructing patient’s epileptic network using iEEG recording. Multi-unit computational model can be used to simulate the seizure dynamics, and evaluate the effects of external excitatory and inhibitory stimulation. By using iEEG and computational model, our study provided a new approach to localize the optimal targets for the potential neuromodulation of these patients.

## Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: http://www.ieeg.org.

## Author contributions

CL and YL proposed the approach. YL performed the data analysis and model simulation. YL wrote the first draft of the manuscript. CL revised the manuscript. YL and CL approved the submitted version.

## Funding

This work was supported by the National Natural Science Foundation of China (61771323) and Research Fund of Liaoning Provincial Natural Science (2021-KF-12-11), Program for Liaoning Innovative Talents in University, Research project of Liaoning Provincial Department of Education (LJGD20200012).

## Acknowledgments

We would like to thank Fang Zhou for her help with data analysis.

## Conflict of interest

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.

## Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

Abouelleil M., Deshpande N., Ali R. (2022). Emerging trends in neuromodulation for treatment of drug-resistant epilepsy. *Front. Pain Res.* 3, 839463. doi:10.3389/fpain.2022.839463

Akaike H. (1971). Autoregressive model fitting for control. *Ann. Inst. Stat. Math.* 23, 163–180. doi:10.1007/bf02479221

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

Benjamin O., Fitzgerald T. H., Ashwin P., Tsaneva-Atanasova K., Chowdhury F., Richardson M. P., et al. (2012). A phenomenological model of seizure initiation suggests network structure may explain seizure frequency in idiopathic generalised epilepsy. *J. Math. Neurosci.* 2, 1. doi:10.1186/2190-8567-2-1

Brodie M. J., Barry S. J. E., Bamagous G. A., Norrie J. D., Kwan P. (2012). Patterns of treatment response in newly diagnosed epilepsy. *Neurology* 78, 1548–1554. doi:10.1212/WNL.0b013e3182563b19

Choi H., Sell R. L., Lenert L., Muenning P., Goodman R. R., Gilliam F. G., et al. (2008). Epilepsy surgery for pharmacoresistant temporal lobe epilepsy: A decision analysis. *JAMA* 300, 2497–2505. doi:10.1001/jama.2008.771

Coben R., Mohammad-Rezazadeh I. (2015). Neural connectivity in epilepsy as measured by Granger causality. *Front. Hum. Neurosci.* 9, 194. doi:10.3389/fnhum.2015.00194

Creaser J., Lin C., Ridler T., Brown J. T., D’Souza W., Seneviratne U., et al. (2002). Domino-like transient dynamics at seizure onset in epilepsy. *PLoS Comput. Biol.* 16, 1008206. doi:10.1371/journal.pcbi.1008206

Davis P., Gaitanis J. (2020). Neuromodulation for the treatment of epilepsy: a review of current approaches and future directions. *Clin. Ther.* 42, 1140–1154. doi:10.1016/j.clinthera.2020.05.017

de Tisi J., Bell G. S., Peacock J. L., McEvoy A. W., Harkness W. F. J., Sander J. W., et al. (2011). The long-term outcome of adult epilepsy surgery, patterns of seizure remission, and relapse: a cohort study. *Lancet* 378, 1388–1395. doi:10.1016/S0140-6736(11)60890-8

Dolan K. T., Spano M. L. (2001). Surrogate for nonlinear time series analysis. *Phys. Rev. E Stat. Nonlin. Soft Matter Phys.* 64, 046128. doi:10.1103/PhysRevE.64.046128

Folloni D., Verhagen L., Mars R. B., Fouragnan E., Constans C., Aubry J., et al. (2019). Manipulation of subcortical and deep cortical activity in the primate brain using transcranial focused ultrasound stimulation. *Neuron* 101, 1109–1116.e5. doi:10.1016/j.neuron.2019.01.019

Franaszczuk P. J., Blinowska K. J., Kowalczyk M. (1985). The application of parametric multichannel spectral estimates in the study of electrical brain activity. *Biol. Cybern.* 51, 239–247. doi:10.1007/BF00337149

Goodfellow M., Glendinning P. (2013). Mechanisms of intermittent state transitions in a coupled heterogeneous oscillator model of epilepsy. *J. Math. Neurosci.* 3, 17. doi:10.1186/2190-8567-3-17

Hosseini S. A. H., Sohrabpour A., He B. (2018). Electromagnetic source imaging using simultaneous scalp EEG and intracranial EEG: An emerging tool for interacting with pathological brain networks. *Clin. Neurophysiol.* 129 (1), 168–187. doi:10.1016/j.clinph.2017.10.027

Janszky J., Janszky I., Schulz R., Hoppe M., Behne F., Pannek H. W., et al. (2005). Temporal lobe epilepsy with hippocampal sclerosis: predictors for long-term surgical outcome. *Brain* 128, 395–404. doi:10.1093/brain/awh358

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

Kovac S., Vakharia V. N., Scott C., Diehl B. (2017). Invasive epilepsy surgery evaluation. *Seizure* 44, 125–136. doi:10.1016/j.seizure.2016.10.016

Kullback S., Leibler R. A. (1951). On information and sufficiency. *Ann. Math. Stat.* 22, 79–86. doi:10.1214/aoms/1177729694

Kwan P., Brodie M. J. (2000). Early identification of refractory epilepsy. *N. Engl. J. Med.* 342, 314–319. doi:10.1056/NEJM200002033420503

Lam A. D., Zepeda R., Cole A. J., Cash S. S. (2016). Widespread changes in network activity allow non-invasive detection of mesial temporal lobe seizures. *Brain* 139, 2679–2693. doi:10.1093/brain/aww198

Li X., Wu Y., Wei M., Guo Y., Yu Z., Wang H., et al. (2021). A novel index of functional connectivity: phase lag based on wilcoxon signed rank test. *Cogn. Neurodyn.* 15, 621–636. doi:10.1007/s11571-020-09646-x

Lin Z., Meng L., Zou J., Zhou W., Huang X., Xue S., et al. (2020). Non-invasive ultrasonic neuromodulation of neuronal excitability for treatment of epilepsy. *Theranostics* 10, 5514–5526. doi:10.7150/thno.40520

Lopes da Silva F., Blanes W., Kalitzin S. N., Parra J., Suffczynski P., Velis D. N. (2003). Epilepsies as dynamical diseases of brain systems: basic models of the transition between normal and epileptic activity. *Epilepsia* 44, 72–83. doi:10.1111/j.0013-9580.2003.12005.x

Paldino M. J., Golriz F., Zhang W., Chu Z. D. (2019). Normalization enhances brain network features that predict individual intelligence in children with epilepsy. *PloS One* 14, e0212901. doi:10.1371/journal.pone.0212901

Pascual-Marqui R. D., Biscay R. J., Bosch-Bayard J., Lehmann D., Kochi K., Kinoshita T., et al. (2014). Assessing direct paths of intracortical causal information flow of oscillatory activity with the isolated effective coherence (iCoh). *Front. Hum. Neurosci.* 8, 1, 12. doi:10.3389/fnhum.2014.00448

Proix T., Jirsa V. K., Bartolomei F., Guye M., Truccolo W. (2018). Predicting the spatiotemporal diversity of seizure propagation and termination in human focal epilepsy. *Nat. Commun.* 9, 1088. doi:10.1038/s41467-018-02973-y

Richardson M. P. (2012). Large scale brain models of epilepsy: dynamics meets connectomics. *J. Neurol. Neurosurg. Psychiatry* 83, 1238–1248. doi:10.1136/jnnp-2011-301944

Rincon N., Barr D., Velez-Ruiz N. (2021). Neuromodulation in drug resistant epilepsy. *Aging Dis.* 12, 1070–1080. doi:10.14336/AD.2021.0211

Ryvlin P., Jehi L. E. (2021). Neuromodulation for refractory epilepsy. *Epilepsy Curr.* 22, 11–17. doi:10.1177/15357597211065587

Saggio M. L., Crisp D., Scott J. M., Karoly P., Kuhlmann L., Nakatani M., et al. (2020). A taxonomy of seizure dynamotypes. *Elife* 9, 55632. doi:10.7554/eLife.55632

Schneider T., Neumaier A. (2001). Algorithm 808: ARfit – a matlab package for the estimation of parameters and eigenmodes of multivariate autoregressive models. *ACM Trans. Math. Softw.* 27, 58–65. doi:10.1145/382043.382316

Schulze-Bonhage A. (2017). Brain stimulation as a neuromodulatory epilepsy therapy. *Seizure* 44, 169–175. doi:10.1016/j.seizure.2016.10.026

Sinha N., Dauwels J., Kaiser M., Cash S. S., Westover M. B., Wang Y., et al. (2017). Predicting neurosurgical outcomes in focal epilepsy patients using computational modelling. *Brain* 140, 319–332. doi:10.1093/brain/aww299

Sinha N., Wang Y., da Silva N., Miserocchi A., McEvoy A. W., de Tisi J., et al. (2021). Structural brain network abnormalities and the probability of seizure recurrence after epilepsy surgery. *Neurology* 96, e758–e771. doi:10.1212/WNL.0000000000011315

Sip V., Guye M., Bartolomei F., Jirsa V. (2022). Computational modeling of seizure spread on a cortical surface. *J. Comput. Neurosci.* 50, 17–31. doi:10.1007/s10827-021-00802-8

Sisterson N. D., Kokkinos V. (2020). Neuromodulation of epilepsy networks. *Neurosurg. Clin. N. Am.* 31, 459–470. doi:10.1016/j.nec.2020.03.009

Stephen L. J., Kelly K., Mohanraj R., Brodie M. J. (2006). Pharmacological outcomes in older people with newly diagnosed epilepsy. *Epilepsy Behav.* 8, 434–437. doi:10.1016/j.yebeh.2005.11.007

Stern J. M., Spivak N. M., Becerra S. A., Kuhn T. P., Korb A. S., Kronemyer D., et al. (2021). Safety of focused ultrasound neuromodulation in humans with temporal lobe epilepsy. *Brain Stimul.* 14, 1022–1031. doi:10.1016/j.brs.2021.06.003

Taylor P. N., Han C. E., Schoene-Bake J., Weber B., Kaiser M. (2015). Structural connectivity changes in temporal lobe epilepsy: spatial features contribute more than topological measures. *Neuroimage. Clin.* 8, 322–328. doi:10.1016/j.nicl.2015.02.004

Terry J. R., Benjamin O., Richardson M. P. (2012). Seizure generation: the role of nodes and networks. *Epilepsia* 53, e166–e169. doi:10.1111/j.1528-1167.2012.03560.x

Trinka E., Cock H., Hesdorffer D., Rossetti A. O., Scheffer I. E., Shinnar S., et al. (2015). A definition and classification of status epilepticus–report of the ILAE task force on classification of status epilepticus. *Epilepsia* 56, 1515–1523. doi:10.1111/epi.13121

Tsuboyama M., Kaye H. L., Rotenberg A. (2020). Review of transcranial magnetic stimulation in epilepsy. *Clin. Ther.* 42, 1155–1168. doi:10.1016/j.clinthera.2020.05.016

van Diessen E., Diederen S. J. H., Braun K. P. J., Jansen F. E., Stam C. J. (2013). Functional and structural brain networks in epilepsy: What have we learned? *Epilepsia* 54, 1855–1865. doi:10.1111/epi.12350

Wendling F., Benquet P., Bartolomei F., Jirsa V. (2016). Computational models of epileptiform activity. *J. Neurosci. Methods* 260, 233–251. doi:10.1016/j.jneumeth.2015.03.027

Wieser H. G., Blume W. T., Fish D., Goldensohn E., Hufnagel A., King D., et al. (2001). Proposal for a new classification of outcome with respect to epileptic seizures following epilepsy surgery. *Epilepsia* 42, 282–286. doi:10.1046/j.1528-1157.2001.35100.x

Wilcoxon F. (1945). Individual comparisons by ranking methods. *Biom. Bull.* 1, 80–83. doi:10.2307/3001968

Wilke C., Worrell G., 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

Keywords: neural computational model, neuromodulation, drug-resistant epilepsy, intracranial EEG, optimal target

Citation: Liu Y and Li C (2022) Localizing targets for neuromodulation in drug-resistant epilepsy using intracranial EEG and computational model. *Front. Physiol.* 13:1015838. doi: 10.3389/fphys.2022.1015838

Received: 10 August 2022; Accepted: 10 October 2022;

Published: 20 October 2022.

Edited by:

Rong Liu, Dalian University of Technology, ChinaReviewed by:

Yangsong Zhang, Southwest University of Science and Technology, ChinaPaulo Ricardo Protachevicz, University of São Paulo, Brazil

Copyright © 2022 Liu and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Chunsheng Li, lichunsheng@sut.edu.cn