PatientsWe used two datasets: (1) Clinical data with 6 patients with 12 seizures from the neurosurgery department of Shanxi Provincial People’s Hospital; (2) HUP data with 10 patients with 20 seizures from the public dataset (Hospital of the University of Pennsylvania).For the Clinical data (Table S1), the study protocols complied with the Code of Ethics of the Declaration of Helsinki and approved by Shanxi Provincial People’s Hospital. All the participants provided written informed consent as approved by the institutional review board. All patients underwent a one-year postoperative follow-up. The sEEG electrode implantation area is based on non-invasive discovery. Electrode placement was planned by a panel of epileptologists according to the clinical protocol with a separate reference electrode placed in the location of seizure onset (Table S2). Each electrode had a diameter of 0.8 mm and consists of 5–18 contacts with a length of 2 mm and a spacing of 1.5 mm at 1000 Hz sampling rate. Under general anesthesia, electrodes were implanted using a Robotic Surgical assistant (ROSA) system.For the HUP data (Table S3), cortical surface electrode configurations, determined by a multidisciplinary team of neurologists and neurosurgeons, consisted of linear and 2D arrays (2.3-mm diameter with 10-mm intercontact spacing). Signals were recorded and digitized at 500 Hz sampling rate and preprocessed to eliminate line noise. Details are available on the public web site (https://openneuro.org/).Preprocessing with the Brainstorm tool29: First we removed the bad channels, then performed ICA to remove any artifacts, bandpass filtered at 0.16–97 Hz, and down-sample from 1000 Hz to 500 Hz. By reducing the sampling rate, important information can be retained while reducing the amount of data, making subsequent analysis and processing more efficient. The bad channels were mainly determined according to sEEG observations and the doctor’s record during the operation. The ICA-based artifact correction approach was applied to separate and remove artifacts in the SEEG data through linear decomposition.Simulation dataResection schemes could be obtained using real data, but the current gold standard for judging their correctness is still EEG analysis by a trained clinician, which is time-consuming, subjective, and might miss important features. We used the “Epileptor” model, proposed by Jirsa30, to generated simulated data to evaluate the accuracy of the resection plan. This model has been applied to explain many dynamic pathways of epileptic seizures31 and demonstrate changes in brain network functional connectivity32. The model is carried out here using the following five equations:$$\:{\dot{x}}_{1,i}={y}_{1,i}-{f}_{1}\left({x}_{1,i},{x}_{2,i}\right)-z+{I}_{1}$$
(1)
$$\:{\dot{y}}_{1,i}=1-5{x}_{1,i}^{2}-{y}_{1,i\:}\:$$
(2)
$$\:{\dot{z}}_{i}=\frac{1}{{\tau\:}_{0}}\left(4\left({x}_{1,i}-{x}_{0}\right)-{z}_{i}-\sum\:_{j=1}^{N}K{C}_{ij}\left({x}_{1,j}-{x}_{1,i}\right)\right)\:\:$$
(3)
$$\:{\dot{x}}_{2,i}=-{y}_{2,i}+{x}_{2,i}-{x}_{2,i}^{3}+{I}_{2}+0.002g\left({x}_{1,i}\right)-0.3\left({z}_{i}-3.5\right)\:\:\:$$
(4)
$$\:{\dot{y}}_{2,i}=\frac{1}{{\tau\:}_{2}}\left(-{y}_{2,i}+{f}_{2}\left({x}_{2,i}\right)\right)\:\:\:\:\:\:\:$$
(5)
where$$\:{f}_{1}\left({x}_{1,i},{x}_{2,i}\right)=\left\{\begin{array}{c}{x}_{1,i}^{3}-3{x}_{1,i}^{2}\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:if{x}_{1,i}<0\\\:\left({x}_{2,i}-0.6{\left({z}_{i}-4\right)}^{2}\right){x}_{1,i}\:\:\:\:\:\:\:\:\:\:if{x}_{1,i}\ge\:0\end{array}\right.\:\:\:\:\:\:\:\:\:\:$$
(6)
$$\:{f}_{2}\left({x}_{2,i}\right)=\left\{\begin{array}{c}0\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:if{x}_{2,i}<-0.25\\\:6\left({x}_{2,i}+0.25\right)\:\:\:\:\:\:\:\:\:\:\:if{x}_{2,i}\ge\:-0.25\end{array}\right.\:\:\:\:\:$$
(7)
$$\:g\left({x}_{1,i}\right)={\int\:}_{-{t}_{0}}^{t}{exp}^{-\gamma\:\left(t-\tau\:\right)}{x}_{1,i}\left(\tau\:\right)dt\:\:\:\:\:\:\:\:\:\:\:$$
(8)
where τ0 = 2857, τ2 = 10, I1 = 3.1, I2 = 0.45, and γ = 0.01. The parameter x0 indicates the excitability of the brain region (see Appendix S1 for details).We used the “Epileptor” model to generate five sets of simulated data to explore the general laws of virtual resection. The data were generated from 588 nodes and 56 nodes affected by seizures. For the five simulations, different parameters were set respectively for the EZ, early propagation and late propagation depth electrodes during each simulation (Table S4). The simulations were generalizability: spread within the same depth electrodes; spread to the contralateral brain; spread between different depth electrodes.Construct high order efficient networkConstruct efficient networkWe construct a directed brain network using transfer entropy (TE) to estimate the flow of information from one time series to another33. TE could also be used as a measure of causality, because it considers the information transfer between variables without assuming that there is a specific relationship between variables, especially for nonlinear systems34,35. The TE method is described in greater detail in Appendix S2:$$\:TE\left(Y\:\to\:X\:\:\right)=\sum\:_{{x}_{i+1},{x}_{i}^{\left(k\right)},{y}_{i}^{\left(l\right)}}p({x}_{i+1},{x}_{i}^{\left(k\right)},{y}_{i}^{\left(l\right)})log\frac{p\left({x}_{i+1}\right|{x}_{i}^{\left(k\right)},{y}_{i}^{\left(l\right)})}{p\left({x}_{i+1}\right|{x}_{i}^{\left(k\right)})}\:\:\:\:\:\:\:\:\:\:\:\:\:$$
(9)
Perform Fisher transformation on the obtained effective matrix. Next, we apply FDR (Benjamin and Hochberg) correction at the 0.001 level to control the false positive rate. Finally, we binarize the matrix to obtain directed graphs, which will serve as input data for high-order TE (HTE) analysis.Capture propagation modeIn order to find the specific seizure propagation pathways, we proposed HTE to describe the propagation of effective connections in brain networks, which can identify characteristic regions connected to specific seed brain regions at different link-step levels.All HTE analyses were conducted at the individual level of each participant. In our framework, step refers to the number of links (edges) connecting an electrode to the path of the target electrode. Considering that epilepsy patients reach a multimodal integrated network in the third stage, only three linking steps were proposed, which could capture three neighbors. This is because we found that the third-order network of epilepsy patients basically includes the EZ region and PZ region27. By using the HTE method, the propagation path of epilepsy could be accurately captured. On this basis, surgical schemes could be formulated, but evaluation criteria are lacking.$$\:{HTE}^{K}=\left\{\begin{array}{c}{HTE}^{1}*{HTE}^{K-1},\:\:\:\:\:k\ge\:2\\\:HTE,\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:k=1\:\:\end{array}\right.\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:$$
(10)
Virtual resection evaluationOn these basis, we first measure network synchronization by calculating the Laplacian matrix of each adjacency matrix36,37. The Laplacian matrix could be interpreted as a measure of how easily information spreads among electrodes in a network38:\(\:D\) is the outdegree matrix of \(\:HTE\).Next, we calculate the asymmetric synchronization as the ratio of the third smallest eigenvalue to the largest eigenvalue of the Laplacian matrix, which quantifies the stability of the synchronization state39:$$\:R={\lambda\:}_{N}/{\lambda\:}_{2}$$
(12)
Where \(\:{\lambda\:}_{N}\) is the largest eigenvalue of the Laplacian matrix and \(\:{\lambda\:}_{2}\)is the second smallest eigenvalue of the Laplacian matrix.To simulate the effect of excision surgery, we used the method of virtual cortical excision, which quantifies the contribution of control centrality as synchronicity40. The calculation process of control centrality: step 1, HTE captured the propagation path; Step 2, gradually removed contacts from the propagation path; Step 3, re-calculate the asymmetric synchronicity of the remaining contact composition area to evaluate its stability. Control centrality can be calculated after remaining each contact or the entire region of interest (composed of many contacts). In this study, based on the propagation path, we cut it off within each layer according to the outdegree (from left to right, from top to bottom), and calculated the control centrality at the contact level. The epilepsy network analysis is shown in Fig. 1, including 2.4 and 2.5.Fig. 1A pipeline for epilepsy network analysis. The sEEG pipeline uses transfer entropy to construct a directed connection matrix, and a higher-order matrix containing a propagation network is obtained by step process. The removed contacts are determined from the higher-order matrix and the control centrality after removing area is calculated.
Virtual resection evaluation based on sEEG propagation network for drug-resistant epilepsy
