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

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):

$\sum _{j=0}^p{A}_{j}x\left(t-j\right)=w\left(t\right),(1)$

where x(t) = [x₁(t), x₂(t), …, x_N(t)] is the vector of EEG N-channel process, p is the order of the model. A₀ is identity matrix, A₁, A₂, … , A_p are the N × N matrices of model coefficients, w(t) = [w₁(t), w₂(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-s}^T$, where x^T is transposed vector of x. We get following equation (Kaminski and Blinowska, 1991):

$R\left(-s\right)+A_{1}R\left(1-s\right)+\cdots +A_{p}R\left(p-s\right)=0,(2)$

where $R(s)=E(x(t),x_{t-s}^T)$ is the covariance matrix with lag 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

$X\left(z\right)=H\left(z\right)W\left(z\right).(3)$

where H(z) is the transfer function. Set z⁻¹ = 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:

$r_{ij}^2\left(f\right)=\frac{{\left|H_{ij}\left(f\right)\right|}^2}{{\sum }_{r=1}^n{\left|H_{ir}\left(f\right)\right|}^2}.(4)$

Note that the value of $r_{ij}^2(f)$ is between 0 and 1.

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:

$P_{a,b}=\frac{\mathrm{cov}\left(a,b\right)}{{\sigma }_a{\sigma }_b},(5)$

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 ${\bar{H}}_{ij}(f)$ of PDC is calculated by the following equation (Baccala and Sameshima, 2001):

${\bar{H}}_{ij}\left(f\right)=I-{\sum }_{j=1}^p{A}_j{e}^{-i2\pi f\mathrm{\Delta }t},(6)$

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):

$PLI=|\mathit{sgn}\left(\sin \left(\mathrm{\Delta }\theta \left(\mathrm{t}\right)\right)\right)|,(7)$

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

$wPLI=\frac{|{\mathrm{A}}_1{\mathrm{A}}_2\sin \left(\mathrm{\Delta }\theta \left(\mathrm{t}\right)\right)|}{\left|{\mathrm{A}}_1{\mathrm{A}}_2\sin \left(\mathrm{\Delta }\theta \left(\mathrm{t}\right)\right)\right|},(8)$

where A₁ and A₂ are the corresponding amplitudes of the s₁ and s₂, 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.

$W_x\left(n,f\right)=\frac{|X\left(n,f\right){|}^2}{{{\sum }_{n,f}|X\left(n,f\right)|}^2},(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:

$D_{KL}\left(W_x,W_y\right)=\sum _{n,f}{W}_x\left(n,f\right)\log \frac{W_x\left(n,f\right)}{W_y\left(n,f\right)}.(10)$

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),

$\left\{\begin{array}{c}{X}_1\left(n\right)=0.95\sqrt{2}{X}_1\left(n-1\right)-0.9025X_1\left(n-2\right)+w_1\left(n\right),\hfill \\ X_2\left(n\right)=0.5X_1\left(n-2\right)+w_2\left(n\right),\hfill \\ X_3\left(n\right)=-0.4X_1\left(n-3\right)+w_3\left(n\right),\hfill \\ X_4\left(n\right)=-0.5X_1\left(n-2\right)+0.25\sqrt{2}{X}_4\left(n-1\right)+0.25\sqrt{2}{X}_5\left(n-1\right)+w_4\left(n\right),\hfill \\ X_5\left(n\right)=-0.25\sqrt{2}{X}_4\left(n-1\right)+0.25\sqrt{2}{X}_5\left(n-1\right)+w_5\left(n\right),\hfill \end{array}\right.(11)$

where X_i, i = 1, 2, … , 5, represents the ith 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₁ is simulated as epileptic node where seizure starts. Node X₂, X₃ and X₄ are neighbor nodes., and node X₄ has bidirectional connection with node X₅, 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):

$\frac{\mathrm{d}\mathrm{z}}{\mathrm{d}\mathrm{t}}=f\left(z\right)\equiv \left(a|\mathrm{z}{|}^4+\mathrm{b}|\mathrm{z}{|}^2+\lambda -1+i\omega \right)z,(12)$

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:

$\frac{\mathrm{d}{\mathrm{z}}_i\left(t\right)}{\mathrm{d}\mathrm{t}}=f\left(z_i\right)+\beta \sum _{j\ne i}^{N}K{G}_{ij}\left(z_j-z_i\right)+\alpha w_i\left(t\right),(13)$

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 λ₀ or λ₁ on non-target or target nodes, respectively. The optimal λ₀ and λ₁ are chosen by grid search of λ₀ − λ₁ pairs. λ₀ is chosen between 0 and 0.5 in step of 0.05, and λ₁ is chosen between 0.5 and 1 in step of 0.05. The ΔT is calculated at each λ₀ − λ₁ pair. The optimal λ₀ and λ₁ are found when there is the largest standard deviation of ΔT for all pairs. First, we set the λ of all nodes to λ₀, and record the T_es as T₀. Then, we change the λ of each target node to λ₁, and record the T_es as T₁. The difference between T₀ and T₁ is the change in escape time (ΔT), $\mathrm{\Delta }T=\left|T_0-T_1\right|$, represents the effectiveness of the neuromodulation applied on a given node.

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

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₁ 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 λ₁ and λ₀ is 0.85 and 0.50, respectively. The result of ΔT is shown in Figure 2C, and the node X₁ is with the highest value. The effectiveness of external modulation on the node X₁ is shown in Figure 2D. The T_es of the network decreases significantly while the external excitatory stimuli is applied on node X₁.

FIGURE 2

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 λ₀ for 11 patients is 0.22 ± 0.08, λ₁ is 0.99 ± 0.01 (mean ± SD).

TABLE 1

TABLE 1. Patient information, surgical results, optimal parameters, and target location.

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 λ₀−λ₁ pair is plotted in Figures 3B,D. The λ₀ = 0.30 and λ₁ = 1.00 for both patients.

FIGURE 3

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.30 and λ₁ = 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.30 and λ₁ = 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

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

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 λ₀ is set to 0.22 for the non-target of the resected network, and λ₁ = 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₁ 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 λ₀ 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 λ₀ for non-target nodes is good in most cases when analyzing the iEEG dataset, and the value λ₁ 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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou J., Meng L., Lin Z., Qiao Y., Tie C., Wang Y., et al. (2020). Ultrasound neuromodulation inhibits seizures in acute epileptic monkeys. iScience 23, 101066. doi:10.1016/j.isci.2020.101066

PubMed Abstract | CrossRef Full Text | Google Scholar