Individualized epidemic spreading models predict epilepsy surgery outcomes: A pseudo-prospective study

Surgical resection is often the most effective treatment to achieve seizure control for patients with drug-resistant focal epilepsy. The surgery requires the generation of an hypothesis of the epileptogenic zone (EZ) by means of extensive presurgical evaluations, and its subsequent removal or disconnection during surgery (Lüders, Najm, Nair, Widdess-Walsh, & Bingman, 2006). Despite extensive investigations, there has only been a slight improvement in prognosis over the past two decades (Baxendale et al., 2019; Jehi et al., 2015), and between 30% to 50% of the patients who undergo surgery continue to have seizures 1 year later, depending on etiology and location of the EZ (Englot et al., 2015). A key conceptual change in recent years is the notion of epileptogenic networks, which takes into account the complex interplay between different brain regions in promoting and inhibiting seizure generation and propagation (Bartolomei et al., 2017; Kramer & Cash, 2012; van Diessen, Diederen, Braun, Jansen, & Stam, 2013). As a consequence, the effect of a given surgery is to be measured against the whole epileptogenic network: a small resection involving heavily connected regions may have widespread effects, but it may also be compensated for by the remaining network (Hebbink, Meijer, Huiskamp, van Gils, & Leijten, 2017; Nissen et al., 2018). This perspective aligns with the commonly accepted view that large-scale brain organization can be regarded as an emerging phenomenon taking place on a complex network, which has spurred numerous data- and model-based studies (Seguin, Jedynak, et al., 2023; Seguin, Sporns, & Zalesky, 2023). Several network-based studies have found group-level differences between seizure-free and nonseizure-free patients (da Silva et al., 2020; Nissen et al., 2018; Taylor et al., 2018), with removal of pathological hub (i.e., central) regions typically associated with seizure-freedom (Nissen et al., 2017). These results highlight the need to consider patient-specific connectivity (van den Heuvel & Sporns, 2019) in order to tailor the surgery to each patient (Gerster et al., 2021).

A data-driven manner to study the relation between individual brain networks and surgical outcomes involves computational models of seizure dynamics, which allow us to simulate seizure propagation in silico. Different resection strategies can be tested on the computational model before the actual surgery (Goodfellow et al., 2016; Hutchings et al., 2015; Jirsa et al., 2017; Laiou et al., 2019; Lopes et al., 2017; Nissen et al., 2021; Olmi, Petkoski, Guye, Bartolomei, & Jirsa, 2019; Proix, Bartolomei, Chauvel, Bernard, & Jirsa, 2014; Sinha et al., 2017; Taylor, Kaiser, & Dauwels, 2014). The models can be fitted to patient-specific data of brain structure and seizure dynamics, allowing us to tailor the resection strategy for each patient. Within this perspective, previous studies have obtained remarkable success at a group level: Sinha et al. (2017) found that the removal of regions identified as epileptogenic according to an EEG-brain network dynamical model predicted surgical outcome with 81.3% accuracy. Proix, Bartolomei, Guye, and Jirsa (2017), using a seizure model known as the epileptor (Jirsa, Stacey, Quilichini, Ivanov, & Bernard, 2014) based on MRI (magnetic resonance imaging) connectivity, found significant differences in the overlap between the model-based propagation zone and the area sampled by iEEG between patients with good (Engel class I) and bad (Engel class III) outcomes at the group level. Subsequent studies also found a better match between the modeled and clinically observed epileptogenic regions for seizure-free than nonseizure-free patients (Makhalova et al., 2022; Vattikonda et al., 2021). Similarly, Sip et al. (2021) simulated patient-specific resection strategies by means of virtual resections, and found that virtual resections in their model correlated with surgical outcome, with larger effects found for patients with good outcome (Engel classes I and II). In an independent study, Goodfellow et al. (2016) also found significant differences in the model prediction between Engel class I and class IV patients, using an electrocorticogram modeling framework.

Following the same rationale, we developed a computational model of seizure propagation and epilepsy surgery based on epidemic spreading dynamics and patient-specific MEG brain connectivity (Millán et al., 2022), to which we refer here as the Epidemic Spreading Seizure and Epilepsy Surgery model (ESSES). Epidemic models describe the spread of an infectious agent through a network. Epidemic processes on fixed networks have a rich mathematical history (Pastor-Satorras, Castellano, Van Mieghem, & Vespignani, 2015) with a plethora of models that can be exploited for epilepsy surgery optimization (Millán et al., 2022; Nissen et al., 2021). Although such models ignore the underlying biophysical processes that lead to seizure generation and propagation, they describe the basic rules that govern spreading processes. In previous studies (Millán et al., 2022, 2023), we found that epidemic spreading models could reproduce stereotypical patterns of seizure propagation as recorded via invasive electroencephalography (iEEG) recordings. Moreover, once fitted with patient-specific data, ESSES could identify alternative resection strategies, either of smaller size or at a different location than the actual surgery (Millán et al., 2022; Nissen et al., 2021). In a more recent study (Millán et al., 2023), we showed that the goodness of fit of ESSES seizures to those recorded via iEEG predicted surgical outcome—with an area under the curve of 88.6%—indicating that ESSES not only reproduces the basic aspects of seizure propagation, but it also captures the differences, either in the location of the resection area relative to the EZ, or intrinsically in the iEEG or MEG data, between patients with good and bad outcome. Importantly, ESSES’s global parameters were defined at the population level, and the model was individualized for each patient via patient-specific MEG networks, which characterized the local spreading probabilities. As a consequence, ESSES can be extended to patients without iEEG recordings, in contrast to previous modeling studies, which typically required the existence of patient-specific iEEG data to individualize the model for each patient (Bernabei et al., 2023; Gunnarsdottir et al., 2022; Makhalova et al., 2022; Proix et al., 2017; Runfola, Sheheitli, Bartolomei, Wang, & Jirsa, 2023; Sinha et al., 2017; Y. Wang et al., 2023). IEEG allows for a highly resolved description of seizure dynamics, but its spatial sampling is sparse and it is highly invasive. Consequently, it is only part of the presurgical evaluation in a selection of patients.

Here we performed a pseudo-prospective blind study (34-patient validation cohort) to validate the clinical applicability of ESSES to (a) identify model-based optimal resection strategies and (b) predict the likelihood of a good outcome after a proposed resection strategy, on a patient-by-patient basis. In order to emulate the clinical presurgical conditions, the research team was blind to the patients’ postsurgical data, namely the resection area and surgical outcome, during ESSES’s analyses, and the multimodal presurgical information available for eachpatient was integrated into ESSES. ESSES can identify resection strategies that perform optimally in the model, that is, by minimizing modeled seizure propagation, for a given resection size. We refer to these resections as optimal resections, in agreement with previous works (An, Bartolomei, Guye, & Jirsa, 2019; Millán et al., 2022; Nissen et al., 2021; Sinha et al., 2017). ESSES can also simulate the effect of a given resection in silico. Within this setup, we tested three hypotheses: (a) seizure-free (SF) patients would have smaller optimal resections than nonseizure-free (NSF) patients, (b) SF patients would have a larger overlap between optimal and planned (clinical) resections, and (c) the planned resection would have a larger effect (in ESSES) for SF than for NSF patients. We found that these three ESSES biomarkers, namely the size of the optimal resection, their overlap with the planned resection, and the effect of the planned resection on ESSES seizures, provided estimates of the likelihood of a good outcome after the surgery, as well as suggesting alternative resection strategies that performed optimally in the model. We envisage that the implementation of a modeling scheme such as ESSES in clinical practice may inform the planning of epilepsy surgery. Different surgical plans can be tested with ESSES for each patient, such that strategies that lead to a large decrease of propagation in the model are more likely to lead to seizure freedom. ESSES may also suggest optimal (alternative) resection strategies, for cases where ESSES predicts a bad outcome with the planned resection. Optimal strategies can then lead to new surgical plans, the effect of which can then be tested in ESSES again.

Here we validated the clinical applicability of ESSES to (a) identify optimal resection strategies that may improve surgical outcomes and (b) provide estimates of the probability of postsurgical seizure freedom, given a surgical plan. The key goal of ESSES is to identify surgical candidates who would have a bad outcome (NSF patients) so that the surgical plan can be adjusted. This study combined a retrospective analysis on a modeling cohort (N = 15) that was used to set the model hyperparameters (following our previous retrospective study (Millán et al., 2023) on this same cohort), and a pseudo-prospective study on a validation cohort (N = 34) to validate ESSES findings and to emulate its clinical application in a blind setup that mimics the clinical presurgical conditions. The researchers were blind to the performed surgery and surgical outcome during the application of ESSES to the validation cohort.

The study was performed as follows:

1. Seizure model: definition and fitting (modeling cohort). An SIR-type of epidemic spreading process modeled seizure propagation over patient-specific brain connectivity. IEEG data from the modeling cohort was used to fit the global parameters of the spreading model so that ESSES-modeled seizures matched those recorded via iEEG, as shown in Figure 1A.

2. Individualized ESSES framework: patient-specific models. ESSES was individualized for each patient: patient-specific MEG brain connectivity defined the network on which ESSES computed seizure propagation. Multimodal patient-specific data, available from presurgical evaluations, defined the seed regions (i.e., the seizure onset regions) based on epileptogenicity or seed probability maps.

3. Alternative resection strategies (aim A). ESSES incorporates an optimization algorithm to determine model-based optimal resection strategies for each patient. These acted as a benchmark against which the planned resection for each patient could be tested. These resections were optimal in the model in the sense that they minimized modeled seizure propagation.

4. Simulation of the planned resection strategy (aim B). The resection plan for each patient was simulated in ESSES with a virtual resection that emulated the actual surgical resection, and the subsequent decrease in seizure propagation was measured.

5. Statistical analyses (aim B). We compared ESSES’s predictions (steps 3 and 4) between patients with good and bad outcome. We defined the NSF class as the positive class for classification and prediction testing.

We modeled seizure propagation by a susceptible-infected-recovered (SIR) epidemic process, as illustrated in Figure 1. The S-I-R states account respectively for the healthy (pre-ictal), ictal and healthy (post-ictal) states, coupled with patient-specific brain connectivity (derived from MEG data) to define the local spreading probabilities. The SIR model describes the spread of an infection from an initial set of infected nodes, the seed regions, to the other nodes in the network, and the recovery of the infected nodes, without reinfections (Barrat, Barthélemy, & Vespignani, 2008; Pastor-Satorras et al., 2015). Here we confined ourselves to one of the simplest compartmental SIR models, using a discrete-time setting where the spreading probability from node i to node j corresponded to the coupling strength w_ij on the patient-specific brain network and where the recovery probability γ was set to be equal for all nodes. The brain network was initially thresholded (by setting the weakest links to zero) at different densities ρ indicating the fraction of nonzero links remaining in the network after thresholding (see Methods section and Supporting Information section 5).

The two control parameters of ESSES are thus the global recovery probability γ and the network density ρ. We followed the inference method presented in our previous study (Millán et al., 2023) to fit the model parameters to iEEG-recorded seizures of the modeling cohort. We note that the modeling framework as presented here differs slightly from the one in Millán et al. (2023), which included an extra parameter to set the global spreading rate. The details of the model fit can be found in the Methods section, and the fit results are reported in the Supporting Information (Supporting Information section 5.2, see also Supporting Information Figures S2 and S3).

The degree of similarity between the ESSES and iEEG seizures was measured with the goodness of fit C(ρ, γ) (Equation 1). The resulting diagram resembled a familiar phase transition (Figure 1B), with an interface of high goodness of fit (yellow regions) corresponding to a roughly constant spreading-to-recovery ratio ρ/γ = const, in agreement with other studies (Moosavi, Jirsa, & Truccolo, 2022). The maximum goodness of fit is indicated by a blue square in Figure 1B and sets the working point of ESSES for the remaining analyses. At this working point, the SF group presented a significantly better fit than the NSF group (p = 0.04, see Supporting Information Table S3 and Supporting Information Figure S3B for details).

A ROC classification analysis indicated a good classification (AUC = 0.79, see Supporting Information Table S4 and Supporting Information Figure S3C) between the SF and NSF groups. At the optimal classification point (Youden criterion, Supporting Information Figure S3D), all NSF patients were correctly identified. The high sensitivity suggests that all patients identified as SF by ESSES could proceed to surgery with high expectations (100% in this group) of a good outcome. On the contrary, patients identified as NSF should be examined further (e.g., by performing further presurgical evaluations or considering other resection plans) as they had a 57% chance of bad outcome with the proposed surgery (to be compared with a 26% chance of bad outcome expected simply from the relative group sizes).

A key ingredient of ESSES is the definition of the epileptogenic or seed regions. Here we defined epileptogenicity or seed probability maps SP_i, indicating the probability that each brain region i gave rise to a seizure. The seed probability maps integrated patient-specific multimodal presurgical information (encoded in the local patient database (Castor Electronic Data Capture, n.d.)) in a quantitative and systematic manner that was adapted for each patient to include the data from the presurgical evaluations that they had undergone (see Methods section and Supporting Information section 4 for details). The resulting seed probability maps for two representative patients (modeling cohort) are illustrated in Figures 3B and D together with the corresponding resection areas (panels A, C). The seed probability maps show wider spatial patterns than the resection areas, and may involve several lobes in both hemispheres. The resection areas for the two cases shown here were contained within the most likely seed regions. In general, the resection areas had a larger seed probability than expected by chance for all patients. We did not find significant differences in the overlap between the resection areas and the seed probability maps between SF and NSF patients (see Supporting Information Figure S1).

ESSES can derive individualized alternative resection strategies—which minimize modeled seizure propagation—via an optimization algorithm based on simulated annealing (Millán et al., 2022; Nissen et al., 2021). The optimization algorithm parameters were set on the modeling cohort data (see Methods section for the algorithm details and Supporting Information section 5.3 and Supporting Information Figures S4 and S5 for the modeling cohort results), and the algorithm was then applied to the validation cohort in a blind setting.

The optimization algorithm searched for resections R of increasing size S(R) that minimized the seed efficiency E_R(seed), that is, the average distance (on the network) from the seed nodes to the other network nodes. This procedure exploits the link between epidemic spreading dynamics and network structure, such that spreading to a region is strongly influenced by its distance to the seed (Pastor-Satorras et al., 2015). In Figure 4A we show the normalized seed efficiency e_R(seed), which is normalized to the seed efficiency in the unresected network so as to diminish differences due to seed extent and initial efficiency. e_R(seed) decreased with the size of the resection for all patients. At the group level, the SF group showed a significantly smaller e_R(seed) than the NSF group (repeated measures ANOVA test, F(19) = 37.95, p < 10⁻⁸⁹) for all considered seed sizes except S(R) = 1. Moreover, the effect of increasing the resection size on e_R(seed) was larger for the SF than for the NSF group (F(19) = 3.78, p < 10⁻⁶).

The actual effect of a resection R on modeled seizure propagation was quantified by measuring the normalized decrease in seizure propagation due to the resection, δIR(R) (Figure 4B), again relative to propagation on the unresected network. Seizure propagation depended heavily on the seed realization such that a bistable regime emerged in which ESSES seizures either propagated macroscopically or died locally (an exemplary case is shown in Supporting Information Figure S4). Thus, results reported here were averaged over 300 independent realizations of the seed regions and SIR dynamics. At the group level, the SF group presented a larger decrease in seizure propagation (F(19) = 25.88, p < 10⁻⁶⁵) and a larger effect of increasing the resection size (F(19) = 2.90, p = 4 · 10⁻⁵). There were large differences in the dependence of δIR(R) on the resection size between different patients. Whereas in the majority of the cases δIR(R) increased roughly exponentially with S(R), for several patients there was an abrupt (discontinuous) jump at a given resection size.

We defined the optimal resection R_op as the one leading to a 90% decrease in seizure propagation, δIR(R_op) = 0.90. The SF group had significantly smaller optimal resections, and these presented a significantly larger overlap with the actual resection strategy Ov(R_op, RA) (see panels C and F of Figure 4 and Table 1), than the NSF group. We found good classification results using either of these variables to classify between the SF and NSF groups (AUC = 0.71, 0.69 respectively for S(R_op) and Ov(R_op, RA); see Figures 4D and G). Both variables led to very similar classification results at the optimal classification point (Youden criterion), correctly identifying 6/8 NSF cases (panels E and H). The classification results for the validation cohort are summarized in Table 2 (see Supporting Information Table S4 for the modeling cohort results).

In summary, these results indicate that the planned resection strategy (accounted for here by the resection area) presented a larger overlap with the optimal resection for patients with good outcome. In particular, 90.0% of SF and 42.9% of NSF patients were correctly classified by Ov(RA, R_op). Remarkably, ESSES could also distinguish between SF and NSF patients without taking into account the information of the surgical plan. In fact, up to 90.4% of SF and 46% of NSF patients were correctly identified by S(R_op) (in relation to only a 76.5% SF-chance and 23.5% NSF-chance according to the group ratios). As this analysis did not depend on the planned resection strategy, a bad prognosis would be indicative of the need to perform a more exhaustive presurgical evaluation, and potentially imply an unavoidable nonseizure-free outcome after any surgery.

Finally, we note that almost equivalent results may be obtained by considering the disconnecting resection, that is, the smallest resection leading to disconnection of the seed, instead of the optimal resection (see Supporting Information section 5.4 and Supporting Information Figure S6). This is due to the strong link between network topology and emergent SIR dynamics, a result that can be used to speed up computations considerably, by using a purely network-based analysis of the effect of different resection strategies.

We simulated the effect of the planned surgery in ESSES for each patient by performing virtual resections of the resection area, which was considered as a proxy for the surgical plan here. We report here on the results for the validation cohort (Figure 5); results for the modeling cohort can be found in the Supporting Information (Supporting Information section 7, Supporting Information Figure S8). As in previous sections, all modeling details had already been set during the modeling step. The effect of the resection strategy on (modeled) seizure propagation, δIR(RA), was significantly larger for the SF than the NSF group (Figure 5B, Table 1). A ROC classification analysis revealed a good classification between the two groups (AUC = 0.78, Figure 5C), and at the optimal point (Youden criterion, black asterisk in panel C) the majority of the patients were correctly identified (Figure 5D, Table 2). In particular, there was a 91.3% chance that a patient classified as SF had a good outcome, and a 54.5% chance that a patient classified as NSF had a bad outcome, compared to a 76.5% and 23.5% chance based on the relative group sizes.

The classification analyses in the previous sections were informed by each patient’s surgical outcome. In a prospective setting the outcome for the patient is not yet known, and thus cannot be used to build the classification model. In order to emulate a true prospective setting, we performed a prediction analysis based on leave-one-out cross-validation. That is, in order to predict the outcome of each patient of the validation cohort, a prediction model was built using data from the remaining 33 cases. Results from this analysis are shown in Figure 6, with the statistical details reported in Table 3. The prediction results were slightly worse than the classification ones (previous sections), particularly for the NSF class where there was a 12.5% reduction in the group size. In any case, respectively, 4, 5, and 5 NSF cases and 19, 18, and 21 SF cases were correctly identified by each ESSES biomarker (Figure 6A). Moreover, 75% of NSF cases (6/8) and only 19.2% (5/26) of SF cases were identified by two or more biomarkers as NSF (Figure 6B). For this cohort, if ESSES predicted a good outcome with at least two markers, there was a 80.8% chance of seizure freedom after the surgery (compared to a 76.5% expectancy of surgery success according to the group rates). Conversely, if the model predicted a bad outcome, then there was a 75% chance that the surgery would fail (compared to a 23.5% expectancy of surgery failure according to the group rates). In clinical practice, a good ESSES prediction could then be interpreted as a large (80.8%) chance of seizure freedom after the surgery and thus support the decision to proceed with surgery. On the contrary, a bad ESSES prediction would indicate a 76.5% chance that the surgery would fail. This may be suggestive of the need of more presurgical evaluations or a different resection strategy, and eventually indicate a low probability of complete seizure freedom after the surgery.

Finally, in order to test whether the information provided by the three biomarkers could be combined to improve the prediction results, we performed a machine learning analysis using an adaptive boosting algorithm with random undersampling and leave-one-out cross-validation (Figures 7A and B). The input variables for the classification algorithm were δIR(RA), S(R_op), and Ov(RA, R_op). We found that even though the accuracy of the model was good (0.71), the machine learning model was biased towards the majority class (SF), with only 35% of NSF cases correctly identified (precision = 0.37, sensitivity = 0.35) and a poor result for F1 = 0.36, even though the considered algorithm (RUSboost) was designed to correct for class imbalance. However, the minority class in our case contained only eight cases, likely preventing the model from being able to generalize. In order to address this issue, we created a combined cohort (N = 49) pooling together the patients from the modeling and validation cohorts (Figures 7C and D). The combined cohort had 12 NSF cases (50% increase), and the new model was able to identify the majority of NSF cases correctly (72% of SF cases and 63% of NSF cases). Even though the accuracy of the model (0.70) did not improve, the remaining measures, which are less affected by class imbalance, did (precision = 0.42, sensitivity = 0.63, F1 = 0.51). Overall, the machine learning model was not able to improve upon the results found using the individual variables (see Table 3), and indeed the prediction was predominantly based only on a single biomarker, namely the effect of the planned resection on the modeled seizures, δIR(RA). Due to the small sample size, we could not determine whether this was due to intrinsic model limitations, suboptimal hyperparameters, or simply a too small group size (particularly of the minority class). Our setup (leave-one-out cross-validation combined with random undersampling) was designed to minimize the effects of the small sample size, but could not avoid them fully.

Personalized models of brain dynamics can aid the treatment of patients with neurological disorders. In this study we presented ESSES (Epidemic Spreading Seizure and Epilepsy Surgery model), a framework to aid epilepsy surgery planning on a patient-by-patient basis. ESSES defines individualized seizure propagation models that integrate multimodal presurgical data, and can propose alternative resection strategies and provide confidence bounds for the probability of success of a given strategy. The implementation of ESSES in clinical practice may thus eventually improve the chances of achieving a good postsurgical outcome.

In this study we proposed a combined setting such that ESSES’ parameters were fitted in a retrospective study (N = 15) using iEEG data of ictal activity, in analogy with previous studies (Goodfellow et al., 2016; Jirsa et al., 2017; Kini et al., 2019; Makhalova et al., 2022; Moosavi et al., 2022; H. E. Wang et al., 2023). We validated that ESSES captured the main aspects of seizure propagation and was able to reproduce the iEEG-recorded seizures, in agreement with our previous studies (Millán et al., 2022, 2023). Remarkably, the goodness of fit of ESSES-modeled seizures to iEEG data could identify patients with a bad outcome with AUC = 0.79, 100% sensitivity, and 57% precision. Such information may be integrated in the presurgical evaluation of the patients for whom iEEG data is available: different resection strategies may be tested as the origin of the ESSES-modeled seizures (Millán et al., 2023), with a low goodness of fit being indicative of a low chance of seizure freedom. In particular, a bad prediction by the model would indicate (in this cohort) a 57% chance of a bad outcome (to be compared with only a 26.7% NSF rate in this cohort). Conversely, all patients identified as SF by the model could proceed to surgery with high expectations (100% in this group) of good outcome.

The novel aspect of this study consisted of a subsequent pseudo-prospective study with an independent cohort and in a blind setting. Importantly, we did not require the presence of iEEG data in the pseudo-prospective study, and instead the multimodal presurgical information available for each patient was integrated into seed probability maps. In this manner ESSES can be adapted to the information available for each patient, in a quantitative and systematic manner. IEEG data is highly invasive and burdensome for the patient, and thus not always part of the presurgical evaluation. For instance, only 19 of the 34 patients of the validation cohort had undergone it. Thus, by not requiring iEEG data ESSES can be applied to a much larger patient population than traditional approaches (Goodfellow et al., 2016; Jirsa et al., 2017; Kini et al., 2019), with the expected wider impact.

ESSES may be applied prospectively as follows. First of all, ESSES may suggest optimal resection strategies, in analogy with previous studies (An et al., 2019; Laiou et al., 2019; Millán et al., 2022; Nissen et al., 2021), with the advantage that all multimodal presurgical information available for each patient is integrated into ESSES, instead of considering only one source used for network reconstruction. We note that these resections are optimal within the framework of the model, and this does therefore not guarantee optimal clinical outcome. Nevertheless, we have found that these virtual resections have good predictive value of surgical outcome. The optimal resection strategy, defined here as the smallest resection leading to a 90% decrease in (modeled) seizure propagation, can be used as a first indicator of the chances of seizure freedom after any surgery. In our pseudo-prospective predictive framework (emulating the presurgical conditions) the size of this resection could predict 50% of patients with bad outcome (Table 3), whereas the relative NSF rate in this group was 23.5%. This result is independent of the resection strategy, and it is completely characterized by the presurgical information available for each patient. Thus, a bad prognosis could indicate that either the presurgical information available is not of sufficient quality, or that the patient is unlikely to be seizure-free with any resection strategy.

ESSES can also provide information about the prognosis after a particular resection by (i) comparing it to the optimal ESSES resection strategy and (ii) quantifying its effect on seizure propagation in the patient-specific ESSES model. Here we found that resections with a larger overlap with the optimal virtual resection were more likely to lead to seizure freedom, in agreement with previous studies (Goodfellow et al., 2016; Kini et al., 2019; Makhalova et al., 2022). Similarly, resections leading to a larger decrease in seizure propagation in ESSES were associated with a larger probability of seizure freedom after the resection, in agreement with other modeling (Goodfellow et al., 2016; Kini et al., 2019) and network-based (Bartolomei et al., 2017; Lopes et al., 2017; Nissen et al., 2017) studies. Here we considered only the planned resection strategy, which was approximated here by the resection area, since this information could be derived in a systematic manner, and this setup allowed us to validate ESSES findings. In a presurgical setting, different strategies could be tested to measure the probability of seizure freedom after each one. In particular, we found that when combining the information from the three model-based biomarkers (namely, the size of the optimal resection, its overlap with the planned resections, and the effect of the planned resection on modeled seizure propagation) we could predict pseudo-prospectively 81% and 75% of SF and NSF cases (see Table 3), whereas the relative group ratios were 76.5% and 23.5%, respectively. Clinically, this implies that if a good prognosis is found by at least two biomarkers, then there is a 91.3% (true negative rate, 21 cases were SF of the 23 predicted by the model) chance that the patient will be seizure-free, and the patient can proceed with the surgery with the knowledge that they will likely have a good outcome. Conversely, a bad prognosis by at least two biomarkers indicates a 55% chance of bad outcome, and may be interpreted as an ESSES suggestion to perform more presurgical testing or consider alternative resection strategies. Importantly, epilepsy surgery may still improve the quality of life of the patient even when complete seizure freedom can not be achieved. Thus, moderate a priori chance of a bad outcome is not necessarily a contraindication for surgery, but it is important in the presurgical counseling of the patients.

Our findings here did not depend on the presence of iEEG data, and even when iEEG data were available we only included a low-resolution description of them. IEEG data does provide the most detailed information of epileptogenic activity, and is it often the most valuable tool to identify the epileptogenic zone or predict surgical outcome for patients with complicated ethiology (Bernabei et al., 2023; Gunnarsdottir et al., 2022; Makhalova et al., 2022; Proix et al., 2017; Runfola et al., 2023; Sinha et al., 2017; Y. Wang et al., 2023). In fact, for the modeling cohort we found the best classification results when using the goodness of fit of ESSES-predicted seizure propagation patterns to the iEEG seizures, in agreement with previous studies (Makhalova et al., 2022). IEEG imaging, however, is burdensome to the patient, has risk of complications, and has limited spatial coverage. A first prediction of surgical outcome could thus be performed with ESSES when the results of noninvasive testing have been obtained, and an iEEG study might be avoided if the model already predicts a good outcome with the existing data.

In summary, we showed here that ESSES could identify patients with good outcome presurgically based on (i) the smaller size of the optimal ESSES resection strategies, (ii) a larger overlap of the planned resection strategy with the optimal ESSES resection, and (iii) a larger effect of the planned resection strategy on decreasing (modeled) seizure propagation. Our findings here indicate that ESSES could be generalized to other patient populations (as we did with the validation cohort), with the only requirement of a patient-specific brain network, and can incorporate multimodal information from the existing presurgical evaluation, in particular, without requiring the presence of iEEG data. The ESSES-based biomarkers identified here could be taken into account during presurgical planning to evaluate the need for more testing, or may lead to the decision to forgo the surgery, if a bad outcome is predicted. This extra information may be particularly valuable for patients with complicated ethiology (e.g., discordant information from different modalities, variable seizure propagation patterns, multiple seizure onset zones), for whom the discussion of whether or not to perform the surgery is challenging.

ESSES consists of different interconnected elements, namely (i) the underlying network structure, (ii) the seizure propagation model (and parameter fitting), (iii) the seizure onset zone model, (iv) the virtual resection model, and (v) the virtual resection optimization algorithm. Each of these different elements was designed to model a particular aspect of epilepsy surgery in a synergistic manner. For instance, the emergent properties of the seizure propagation model (the SIR model) led the design of the virtual resection optimization algorithm. At the same time, the modular organization of the framework allows for the independent improvement or modification of each of the modules. In fact, different modules were developed and analyzed in detail in our previous studies. For instance, the virtual resection algorithm model was initially designed in Nissen et al. (2021) and improved in Millán et al. (2022), whereas the seizure propagation and parameter fitting model as used here was mainly defined in Millán et al. (2023). Below we discuss the main modeling considerations and results for each ESSES module.

As the underlying network structure, we considered MEG-derived whole-brain networks as a proxy for structural connectivity, following our previous works (Millán et al., 2022, 2023), and in contrast with other works (An et al., 2019; Jirsa et al., 2017; Nissen et al., 2021; Sip et al., 2021). MEG provides highly temporally resolved information with good spatial resolution and uniform coverage. Our previous studies showed that MEG networks based on the amplitude envelope correlation (AEC) can integrate information from both short-range structural connections (by not correcting for volume conduction) and long-range functional coupling. Thus, AEC-MEG networks can be used as a cost-effective proxy for structural connectivity (Millán et al., 2022) with much lower computational cost than DWI (diffusion-weighted imaging) networks, while also being more sensitive to long-range connections, in particular interhemispheric ones, that may often be missed by DWI (Chen et al., 2015). It would be an interesting question for future studies to discriminate whether structural or functional connections drive seizure propagation, in analogy to recent studies on the spreading of abnormal proteins associated with Alzheimer’s disease (Schoonhoven et al., 2023).

The MEG networks were thresholded at different levels to prune out spurious connections, following previous studies (Millán et al., 2022, 2023; Nissen et al., 2021; Schoonhoven et al., 2023). This requires the use of an arbitrary threshold, which we fitted to the iEEG data. In all cases we considered sparse networks (the maximum density considered was 0.35), and the operating point of ESSES was set at a very low density (0.03). This small density prevented weak or negative correlations from being included in the thresholded network. The proposed thresholding method can become a limitation if denser networks, including more connections, are considered.

ESSES was based on a simple epidemic spreading model, the SIR model. Epidemic spreading models, such as the SIR or SIS (susceptible-infected-susceptible) models, describe the basic aspects of spreading phenomena on networked systems (Pastor-Satorras et al., 2015), and have been used to describe other neurophysiological processes before, such as the spreading of pathological proteins on brain networks (Peraza et al., 2019; Schoonhoven et al., 2023) or the relation between brain structure and function (Stam et al., 2016). Epidemic spreading models have been extensively studied on different network substrates (Pastor-Satorras et al., 2015) and are supported by a well-grounded mathematical and computational framework that we can use to our advantage in the context ofepilepsy surgery. For instance, from an epidemic spreading perspective, it is to be expected that hub removal plays a major role in the decrease of seizure propagation, as found experimentally (Lopes et al., 2017; Nissen et al., 2017), with the spreading threshold heavily influenced by the existence of hubs (Pastor-Satorras et al., 2015). This theoretical background guided the design of an efficient virtual resection optimization algorithm, such that the decrease in seizure propagation after a virtual resection could be approximated by the decrease of centrality of the seed regions.

As we showed here and in previous works, epidemic spreading models can also reproduce the fundamental aspects of seizure propagation at the whole-brain level in epilepsy patients (Millán et al., 2022, 2023). As ESSES’s working point we chose here the values of the global parameters that led to the maximum average goodness of fit of the modeling cohort (Figure 1). Importantly, ESSES was still individualized for each patient by means of the patient-specific brain connectivity, setting the local spreading probabilities, and the patient-specific seed regions (based on the seed probability maps built with multimodal presurgical information). As we showed in our previous study (Millán et al., 2023) and in the Supporting Information here (Supporting Information section 5.2), by not individualizing the global model parameters (namely ρ and γ) for each patient, we were able to reduce noise effects by integrating together ictal data from different patients. Moreover, this formulation allowed us to generalize ESSES to patients for whom iEEG seizure propagation patterns were not available.

Our findings in this study indicated that the iEEG seizure propagation patterns were significantly better explained by ESSES for SF patients, and in fact all NSF cases could be identified by a bad ESSES fit, and 73% of the SF cases by a good fit. There are several possible explanations for these findings. Given that the epidemic seed was based on the resection area for each patient in this part of the analyses, a simple explanation is that the resection strategy might have been better for SF patients given the existing information. However, the difference could also arise from the iEEG data: the sampling may have been inadequate for NSF patients (Sip et al., 2021), or these may have presented seizure dynamotypes (Saggio et al., 2020) that were not well explained by the considered epidemic spreading model (SIR model). The fact that the optimization of virtual resections analysis—which did not depend on the clinical resection area—also found differences between the SF and NSF groups points towards an intrinsic difference between the presurgical data of the two groups, and not only to a suboptimal surgical strategy for the NSF group.

The next ingredient of ESSES was the definition of the seizure onset zone in the model, that is, the set of brain regions from which seizures originate. In this study we presented a method to combine the multimodal presurgical information available for each patient into seed probability maps. This setup thus emulated the clinical situation prior to the surgery, where a surgical strategy has been devised based on the information that is available from the presurgical evaluation. It would also allow for flexibility in the clinical application of ESSES: if more evaluations become available these could be readily integrated into the seed probability map to update ESSES’s results.

The final key ingredients of ESSES were the simulation and optimization of resection strategies. Here we considered a node-based resection such that the resected nodes were disconnected from the network. This approach, however, does not take into account possible widespread effects or plasticity mechanisms, which could also be included into the model (Demuru et al., 2020). The virtual resection optimization algorithm was originally validated in our previous studies (Millán et al., 2022; Nissen et al., 2021). Given that optimizing virtual resections is highly computationally demanding, the algorithm took advantage of the mathematical link between network structure and SIR dynamics to reduce the dynamics-based optimization problem (i.e., finding the resection leading to a minimum seizure propagation) into a network optimization problem (i.e., finding the resection leading to a minimum seed efficiency). This was also motivated by our previous finding that the effect of a resection on the model depended strongly on the centrality of the seed regions after the resection (Millán et al., 2022; Nissen et al., 2021). In particular, Nissen et al. (2021) found that removing connections to the network hubs was the most efficient way to decrease seizure propagation, whereas Millán et al. (2022) verified a strong correlation between a decrease in closeness centrality of the seed and a decrease in seizure propagation following a virtual resection. The effect of a resection on seizure propagation is also influenced by other network and model properties, and as a consequence the optimal network-based and SIR-based resections may differ slightly (Millán et al., 2022). However, the intrinsic noise in the seed definition, in the seed probability maps, and in the actual origin and propagation patterns of iEEG-recorded seizures created variability in the clinical data that absorbed the differences between the network-based and SIR-based optimal resections (which we previously found to be small anyway (Millán et al., 2022)).

The virtual resection optimization algorithm considered here imposed no conditions on the location of the resected regions, nor did it force that the resection strategy was made up of only one set of adjacent regions. Conditions on the resection strategies could be imposed, such as preserving eloquent cortex or forbidding bi-hemispheric resections (An et al., 2019; Laiou et al., 2019). This would limit the dimensionality of the space of possible resection strategies and simplify the computations. However, by not imposing any conditions here we derived an optimal ESSES resection against which other, perhaps clinically more realistic, strategies could be tested (by, e.g., measuring their overlap as we did here).

There are inherent limitations in the modeling of virtual resections, as the findings cannot be directly tested and we often rely on retrospective data. Here we have attempted to simulate how an epilepsy surgery model could be used in the clinic, that is, prospectively, by considering only the presurgical information that is typically available to the clinical team. However, the optimal resections suggested by ESSES can still not be tested in practice, and in fact can only be considered optimal within the context of the model. Only long-term testing of the framework in the clinic can truly validate the use of computational models in epilepsy surgery.

ESSES is an abstraction of seizure dynamics that does not aim to reproduce the detailed biophysiological processes involved in seizure generation and propagation, but aims to focus only on the most relevant features of seizure propagation (Millán et al., 2022, 2023; Nissen et al., 2021; Sip et al., 2021). In order to validate ESSES as a framework to simulate seizures, we compared the modeled seizures with those recorded via iEEG. This required, however, a simplified representation of the iEEG data. In particular, as there was no intrinsic timescale in the SIR model, and to avoid introducing an arbitrary one, we reduced the iEEG data to a pattern that describes the activation order of the sampled ROIs. Furthermore, even if ESSES provides a good representation of the iEEG seizures, extrapolating these results to the simulation of the effect of a resection is not trivial. Moreover, our virtual resection technique assumed that the effect of a surgery could be approximated simply by removing or disconnecting the resected regions, whereas in practice widespread effects and compensation mechanisms are expected (Demuru et al., 2020). Here we validated ESSES’ results against postsurgical outcome, but seizure freedom is not a perfect gold standard either. For instance, in cases with a good outcome a smaller resection could potentially also have led to seizure freedom (Millán et al., 2022; Nissen et al., 2021).

All modeling frameworks are affected by the need to (sometimes arbitrarily) choose modeling parameters, which go from the data reduction process to the choices of thresholds and metrics for the final analyses. Here we considered well-established data preprocessing techniques (Hillebrand et al., 2016). ESSES was validated in previous studies (Millán et al., 2022, 2023; Nissen et al., 2021), and importantly we found that the results held for an independent cohort, and that modeling details (such as the simulation algorithm for the SIR model) did not affect the main results (Millán et al., 2023). A simple model to simulate seizure propagation (the SIR model) also reduced the number of modeling parameters so that the findings could be more easily generalized. Some arbitrary choices were still needed, such as the definition of the 90% threshold to select the optimal resection strategy. However, we validated that similar results were obtained when another resection (the disconnecting resection) was considered.

The seed probability maps were based on an existing low-resolution database (Castor Electronic Data Capture, n.d.). Seed regions were consequently widespread over the network. This also led to a large variability in the results of different simulations for each patient (see, for instance, Figures 4A and B and 5A), as these depended strongly on the seed realization. In order to improve the resolution of the model and minimize noise, the data from each modality could be integrated directly into the model, skipping the 34-region description in the database.

Finally, a limitation of this study is the small size of the nonseizure-free group, with only four cases in the modeling cohort and eight in the validation cohort. This small size limited the classification and prediction analyses, and prevented us from building a more sophisticated machine learning model based on our analysis. With the proposed leave-one-out cross-validation method, combined with random undersampling and a small input space (only three data points per patient), we attempted to overcome these limitations, but we were not able to improve upon the simpler ROC-basedprediction results. Future studies involving more than one center have the potential to at least diminish this limitation.

Individualized computational models of seizure propagation and epilepsy surgery based on patient-specific brain connectivity can reproduce individual iEEG seizure propagation patterns and aid epilepsy surgery planning by proposing alternative resection strategies and providing estimates on the likelihood of seizure freedom after the surgery. Here we presented the ESSES framework for seizure propagation and epilepsy surgery. ESSES combines SIR epidemic spreading dynamics over patient-specific MEG brain connectivity with a virtual resection framework. We defined a method to derive patient-specific regional epileptogenicity maps from the presurgical evaluations of the patients in a systematic and quantitative manner, and integrated them into ESSES. We performed a pseudo-prospective study emulating the use of ESSES in clinical practice, prior to surgery. In the pseudo-prospective analyses we did not require the presence of iEEG data, demonstrating that the model could be applied to larger patient populations. We found that the goodness of fit of ESSES to the iEEG seizures (in a retrospective study), the effect of the planned resection strategy, as well as the size of ESSES optimal resections and their overlap with the planned resection, predicted surgical outcome with 0.68–0.76 AUC and 0.50–0.63 sensitivity to identify nonseizure-free patients. Our results thus prescribe the use of ESSES during the presurgical evaluation to evaluate the need for further presurgical testing on a case-by-case basis or, conversely, support the decision to proceed with surgery in the case of a good outcome prediction. For cases where a bad outcome is predicted, the surgical plan may be altered to include ESSES’s results.

The general design of the study is detailed in Figure 2 and Supporting Information Figure S7. Namely, we first set the hyperparameters of ESSES using a modeling cohort (N = 15) for which seizure propagation patterns derived from iEEG recordings were available. Then, ESSES was fitted with multimodal patient-specific data (in the form of seed probability maps) and it was used to (a) identify optimal resection strategies for each patient and (b) predict the chance of a good outcome after a given resection. Then, ESSES was applied to a validation cohort (N = 34) in a pseudo-prospective analysis with a blind setting to emulate the presurgical conditions. That is, during the application of ESSES to determine optimal resection strategies, the researchers were blind to the actual clinical resection and surgical outcome of each patient. This data was subsequently de-blinded in two stages. First, the resection areas were obtained to be used as a proxy for the surgical plan of each patient to (a) compare them with ESSES’s optimal resection strategy and (b) simulate the effect of the surgical plan in ESSES. Finally, we de-blinded the 1-year surgical outcome to enable a statistical validation of the results.

We included two patient groups in this study, the modeling cohort for the model definition (retrospective study) and the validation cohort for the pseudo-prospective validation. All patients had undergone resective surgery for epilepsy at the Amsterdam University Medical Center, location VUmc, between 2013 and 2019. All patients had received an MEG recording, and underwent pre- and postsurgical magnetic resonance imaging (MRI). All patients gave written informed consent and the study was performed in accordance with the Declaration of Helsinki and approved by the VUmc Medical Ethics Committee. The excluding criterion was the existence of a prior brain surgery.

Both patient groups were heterogeneous with temporal and extratemporal resection locations and different etiology (see Supporting Information Tables S1 and S2 for details). Surgical outcome was classified according to the Engel classification at least 1 year after the surgery (Engel, 1993). Patients with Engel class 1A were labeled as seizure-free (SF), and patients with any other class were labeled as nonseizure-free (NSF). The modeling cohort consisted of 15 patients (4 NSF, 11 females) who had also undergone an iEEG (invasive electroencephalography) study, including postimplantation CT-scans. This same cohort was already included in Millán et al. (2023) and partially in Millán et al. (2022). The validation cohort consisted of 34 patients (8 NSF, 13 females). No extra requirements (other than the presence of an MEG recording of sufficient quality) were placed. In order to maintain the pseudo-prospective setting, the research team was blind to the resection area and outcome of the validation cohort patients. In order to perform the final analyses, for which this information was needed, the data was coded to avoid identification. For two cases of the validation cohort (cases 2 and 9) the data of surgical outcome was de-blinded together with the data of the resection area as the research team became aware of a subsequent resective surgery (indicative of a bad outcome of the first surgery).

Seizure propagation was modeled on the patient-specific brain networks, as derived from MEG data, for both cohorts (see Supporting Information Figure S7). For each patient, a 10 to 15 minutes eyes-closed resting-state (supine position) MEG recording was used to derive broadband (0.5–48.0 Hz) MEG functional connectivity. All instrumental and methodological details were equal to our previous studies (Millán et al., 2022, 2023) and are detailed in the Supporting Information (Supporting Information section 2). Functional networks were generated considering each of the 246 ROIs of the Brainnetome (BNA) atlas (Fan et al., 2016) as nodes. The elements w_ij of the connectivity matrix, indicating the strength of the connection between ROIs i and j, were estimated by the AEC (Amplitude Envelop Correlation) (Brookes et al., 2011; Bruns, Eckhorn, Jokeit, & Ebner, 2000; Colclough et al., 2016; Hipp, Hawellek, Corbetta, Siegel, & Engel, 2012), without including a correction for volume conduction. The uncorrected AEC maintains information about the structural connections, which are mainly determined by the distance between each ROI pair, by not correcting for volume conduction. We validated the relationship between AEC-MEG and structural networks in a previous study (Millán et al., 2022) by comparing them with a well-validated model for structural connectivity: the exponential distance rule (EDR) network. Based on animal studies, the EDR specifies that the weights of structural connections in the brain, w_ij, decay exponentially with the distance between the ROIs d_ij (Ercsey-Ravasz et al., 2013; Gămănuţ et al., 2018; Theodoni et al., 2022), that is, w_ij ∝ exp(−αd_ij). Recent studies have corroborated this behavior also in human structural connectivity (Deco & Kringelbach, 2020; Deco et al., 2021; Roberts, Perry, Roberts, Mitchell, & Breakspear, 2017), although the EDR cannot capture all details of white matter connectivity, as this is not isotropic (Betzel & Bassett, 2018; Jbabdi, Sotiropoulos, Haber, Van Essen, & Behrens, 2015; Markov et al., 2013), and includes long-range connections that are missed by the EDR (Roberts et al., 2016). However, the EDR is enough to capture the overall scaling of structural connections with the distance as observed in the human structural connectome. In Millán et al. (2022) we validated that AEC-MEG networks were strongly correlated (R² = 0.50) with the corresponding EDR networks, therefore showing that AEC-MEG reproduces at least partially the overall organization of structural connectivity. Moreover, AEC-MEG networks also include long-range connections that may promote seizure propagation, but that may be missing from structural (i.e., DWI) networks (Jones, Knösche, & Turner, 2013; Reveley et al., 2015). Thus, uncorrected AEC-MEG networks are a convenient way to construct a network that resembles a structural network and includes long-range connections.

AEC values were rescaled between 0 (perfect anticorrelation) and 1 (perfect correlation), with 0.5 indicating no coupling (Briels et al., 2020). Functional networks were thresholded at different network densities ρ indicating the fraction of links remaining in the network. We note that the networks were thresholded but not binearized, so that w_ij could take values between 0 and 1. The density thresholds were chosen to be logarithmically distributed between 0.01 to 0.35. The weakest nonzero link included in the network had an average weight of 0.54 (range: 0.52–0.56) for ρ = 0.35. At ESSES’s operating point (best model fit) the density was ρ = 0.03, and the weakest nonzero weight was 0.71 (range: 0.67–0.76).

The resection area (RA) was determined from the 3-month postoperative MRI. For the modeling cohort the resection areas were obtained as part of two previous studies (Millán et al., 2022, 2023). For the validation cohort, to maintain a completely blind setting for the first analysis (Optimization of alternative resections), the resection areas were obtained during a second preprocessing step, as described in Figure 2. Cases 9 and 20 of the validation cohort underwent the postoperative MRI on a different MRI scanner at their resection center, respectively, 1 day and 3 weeks after the surgery. Case 9 also lacked a 3-month postoperative MRI; an MRI from 2 years after the surgery was used instead.

The postresection MRIs were coregistered to the preoperative MRI using FSP FLIRT (version 4.1.6) 12 parameter affine transformation. The resection area was then visually identified and assigned to the corresponding BNA ROIs, namely those for which the centroid had been removed during surgery.

Patients in the modeling cohort underwent invasive EEG recordings using stereotactic electrode implantation as described in Millán et al. (2023). One characteristic iEEG-recorded seizure from each patient was used to derive a seizure propagation pattern in terms on the BNA ROIs, the iEEG seizure pattern, as described in Millán et al. (2023) and in the Supporting Information section 3.

ESSES was based on our previous studies (Millán et al., 2022, 2023; Nissen et al., 2021) where we showed that simple epidemic spreading models could reproduce the spatiotemporal seizure propagation patterns derived from invasive EEG recordings, and that they could be used to simulate the effect of different resection strategies in silico. ESSES was based on a well-known epidemic spreading model: the susceptible-infected-recovered (SIR) model (Pastor-Satorras et al., 2015), which was simulated on the patient-specific MEG brain network. The SIR model simulated the propagation of ictal activity from a set of seed regions that were set to be infected at the beginning of the simulation to the remaining nodes in the network, and the subsequent recovery of infected nodes. The SIR dynamics were defined by two parameters: the probability β_ij that each infected node i propagates the infection to a neighbor j (S → I), and the probability γ_i that each infected node i recovers (I → R). For simplicity, we considered here a global recovery probability γ_i = γ, and spreading probabilities given by the MEG network connectivity: β_ij = w_ij. Thus, the spreading rate was determined by the density of connections in network ρ. The two control parameters of ESSES are thus the network density ρ, and the recovery probability γ. Depending on the network structure, the epidemics can show different spatiotemporal spreading profiles described by the probability p_i(t) that each ROI i becomes infected at step t.

We estimated C for a range of values ρ and γ logarithmically distributed (between 0.01 and 0.35 for ρ and between 0.01 and 1.00 for γ), considering N_R = 10⁴ iterations of the SIR dynamics 10 times in order to determine average C values and their fluctuation for each patient. We then found the parameter set that maximized C for each patient (see Supporting Information section 5.2 and Supporting Information Figure S2) and at the group level (Figure 1A). The model parameters that lead to the best fit at the population level defined the ESSES model and were carried over to the pseudo-prospective analyses. Importantly, even though the SIR global parameters were set equal for all patients, ESSES was individualized for each patient by means of their patient-specific MEG brain connectivity, which defined the spreading probabilities, and their patient-specific seed probability map, which defined the seed regions.

The SIR dynamics was simulated by an adaptive Monte Carlo method (the BKL algorithm) in Matlab in discrete time, such that at each time step one new node became infected. N_R = 10⁴ iterations of the dynamics were run for each model configuration in all analyses.

We built seed probability maps indicating the probability that each ROI started a seizure, for each patient of both cohorts. This is a key difference with our previous studies, where the seed regions were either derived from the resection area (Millán et al., 2022, 2023; Nissen et al., 2021), which can only be known after the surgery, or from the iEEG data (Millán et al., 2022, 2023). Here we defined a framework to integrate data from the different presurgical evaluations that were available for each patient, which was encoded in an existing database (Castor EDC, Ciwit B.V., Amsterdam (Castor Electronic Data Capture, n.d.)).

To compute the seed-probability maps, we considered the information available from 6 presurgical modalities: (i) presence of ictal activity in EEG, (ii) MRI lesions, (iii) MEG abnormalities, (iv) PET lesions, (v) SPECT abnormalities, and (vi) iEEG recordings of ictal activity. All patients had undergone an EEG, MRI, and MEG study, but not all of them presented PET, SPECT, or iEEG data. The presence (1) or absence (0) of data of each modality was encoded in a variable D_m = 0, 1, m = 1, 2, …, 6, for each patient.

The database included information at the level of 34 regions, consisting of six frontal regions (fronto-orbital, frontal-basal, frontal-parasagitaal, frontal-periventricular, frontal-lateral, frontal-operculum), six temporal regions (hippocampus, amygdala, uncus, anterior-neocortical, posterior-neocortical, gyrus-parahippocampalis), two insular regions (anterior and posterior insula), one central, one parietal, and one occipital region, for each hemisphere. The temporal and frontal lobes are the most often involved in EZ and resection strategies, and thus are described in more detail in the database.

A clinician (ECWvS) defined a unique projection of the regions in the database on to the BNA ROIs. In most cases the database regions corresponded to well-defined gyri that are also well described in the BNA documentation. A table describing the projection is included as Supporting Information. We projected the abnormality map A_i from the low-resolution description into the BNA atlas to obtain the seed probability maps SP_i, with i = 1, 2, …, 246. Given that the description provided by the database was broad and homogeneous (i.e., the considered ROIs are much larger than the BNA ROIs), and that cooccurrence of abnormalities in different modalities is a strong indicator of the epileptogenic zone, we included a rescaling factor R to produce more focal seed probability maps: SP_i = (A_j)^R, where j is the region in the database corresponding to the BNA ROI i. We found that for R > 2 the results did not depend strongly on R, and report here for R = 3.

In the second virtual resection study, we simulated the effect of the planned resection for each patient, to measure its effectiveness in reducing seizure propagation. The resection area was used as a proxy for the resection strategy (Figure 2: Simulation of the resection plan), since it could be derived in a systematic manner from the data.

For all virtual resection analyses the seed regions were derived from the patient-specific seed probability maps, and the underlying network was given from the patient-specific MEG network as before. In order to obtain precise results, the effect of each resection was averaged over 300 independent realizations of the seed regions from the seed-probability maps. As described in Figure 2, for the validation cohort we first performed the optimization of alternative resections in a blind setting. Then the resection areas were de-blinded and used as a proxy of the planed resection strategy to (i) quantify the overlap of ESSES’s optimal resections with the resection strategy and (ii) measure the effect of the planed resection in decreasing (modeled) seizure propagation. Finally the 1-year postoperative outcome was also de-blinded and used for the statistical analyses.

The weighted correlation coefficient was used to determine the correlation between the iEEG and ESSES seizure propagation patterns for the modeling cohort. In all analyses, for comparisons between SF and NSF patients, we used a two-sided Wilcoxon ranksum test. Significance thresholds for statistical comparisons were set at p < 0.05.

We performed receiver-operating characteristic (ROC) curve analyses to study the patient classification based on (i) the goodness of fit of the model (modeling cohort), (ii) the size of optimal and disconnecting resections (modeling and validation cohorts), (iii) the overlap between optimal resections and the planed resection (modeling and validation cohorts), and (iv) the effect of the planed resection on modeled seizure propagation (modeling and validation cohorts). A positive result was defined as bad outcome (nonseizure-free, NSF) classification.

In order to account for the noise in the SIR model, the spreading dynamics were averaged over 10⁴ iterations of the SIR dynamics to derive each ESSES seizure pattern. The model fit analyses were repeated 10 times to obtain averaged values. For the virtual resection analyses we performed 300 independent realizations of the seed regions and SIR dynamics. Each seed realization was used to measure seizure propagation in the original (before any resections) network and after the selected resection of each size. For the optimization of resections analysis we also ran the simulated annealing algorithm 10 times for each resection size and selected the iteration that led to the minimal seed efficiency.

For the classification analyses we report the accuracy = (TP + TN)/(TP + FP + FN + TN), precision = TP/(TP + FP), sensitivity = TP/(TP + FN), F1 statistic (harmonic mean between precision and sensitivity) = 2TP/(2TP + FP + FN), and area under the curve (AUC). For the prediction analyses, we built a predictive model for each patient using the data from the remaining patients, in a leave-one-out cross-validation-type setting. The predictive model compounded the prediction results from these N = 34 models. We measured its accuracy, precision, sensitivity and F1 statistic.

In the final analysis of the study we performed a predictive machine learning analysis based on the AdaboostM1 algorithm (Matlab 2018) combined with random undersampling. AdaBoost is an adaptive boosting machine learning algorithm in which the weights of misclassified instances are adjusted iteratively to improve the model. By combining adaptive boosting with random undersampling of the majority class (SF group), the classification algorithm effectively addresses class imbalance and reduces bias to the majority class and overfitting risks (AdaboostM1 - Matlab 2018, n.d.; Friedman, Hastie, & Tibshirani, 2000).

For each patient, three variables were considered as input for the prediction analysis: the size of the optimal resection S(R_op), its overlap with the resection area Ov(R_op, RA), and the effect of the resection strategy on modeled seizure propagation δIR(RA). The goal of the machine learning algorithm was to predict surgical outcome. Due to the small cohort size, we performed a leave-one-out-cross-validation procedure, such that N_pat different training sets were created, each leaving out one patient, which was then used to test the prediction model. The training sets were formed by randomly undersampling the majority class (SF) to the size of the minority (NSF) class. The small cohort size also prevented us from including a validation set and performing parameter tuning. Thus, we used default hyperparameters of AdaboostM1 (see AdaboostM1 - Matlab 2018, n.d. for details): the number of learners in each model was set equal tothe group size minus one, the learning rate was set to 1.0 (default) and results were averaged over 10 iterations of the undersampling and AdaboostM1 procedures for each classification model. The machine learning analysis was performed twice: first considering only the patients in the validation cohort (N_pat = 34), and secondly considering all patients (combined cohort, N_pat = 49).

The data used for this manuscript are not publicly available because the patients did not consent for the sharing of their clinically obtained data. Requests to access to the datasets should be directed to the corresponding author. All user-developed codes are publicly available on Github: https://github.com/anapmillan/ESSES.

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00361.

Ana P. Millán: Conceptualization; Data curation; Formal analysis; Investigation; Methodology; Project administration; Software; Validation; Visualization; Writing – original draft; Writing – review & editing. Elisabeth C. W. van Straaten: Conceptualization; Data curation; Funding acquisition; Investigation; Methodology; Project administration; Resources; Supervision; Writing – review & editing. Cornelis J. Stam: Conceptualization; Investigation; Methodology; Project administration; Resources; Supervision; Writing – review & editing. Ida A. Nissen: Conceptualization; Funding acquisition; Software; Writing – review & editing. Sander Idema: Conceptualization; Funding acquisition; Writing – review & editing. Piet Van Mieghem: Conceptualization; Funding acquisition; Supervision; Writing – review & editing. Arjan Hillebrand: Conceptualization; Resources; Software; Supervision; Writing – review & editing.

Ana P. Millán, ZonMW (https://dx.doi.org/10.13039/501100001826), Award ID: 95105006. Ida A. Nissen, ZonMW (https://dx.doi.org/10.13039/501100001826), Award ID: 95105006. Ana P. Millán, Epilepsiefonds (https://dx.doi.org/10.13039/501100006117), Award ID: 95105006. Ida A. Nissen, Epilepsiefonds (https://dx.doi.org/10.13039/501100006117), Award ID: 95105006. Piet Van Mieghem, H2020 European Research Council (https://dx.doi.org/10.13039/100010663), Award ID: 101019718. Ana P. Millán, Ministerio de Ciencia, Innovación y Universidades (https://dx.doi.org/10.13039/100014440), Award ID: RYC2021-031241-I. Ana P. Millán, Ministerio de Ciencia e Innovación (https://dx.doi.org/10.13039/501100004837), Award ID: PID2020-113681GB-I00.

Epileptogenic zone:

Region(s) in the brain that needs to be removed to stop the occurrence of epileptic seizures.

iEEG:

Invasive electroencephalography, consisting of the recording of local brain activity through the generated electrical fields via depth electrodes that are placed inside the brain.

Virtual resection:

Simulation of a resection strategy in a computational model by removing or disconnecting the regions included in the resection.

MEG:

Magnetoencephalography, noninvasive functional recording that measures the magnetic field generated by brain activity.

Resection area:

Set of brain regions removed during epilepsy surgery.

SIR model:

Susceptible-infected-recovered model, epidemiological model that describes the spreading of an agent over a network, accounting for the infection and recovery processes.

Seizure onset zone:

Brain region (or regions) from which seizure are generated.