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

2.1 Seizure propagation as an epidemic spreading process

We modeled seizure propagation by a Susceptible-Infected-Recovered (SIR) epidemic process, where 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 re-infections [31, 37]. Here we confined 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 non-zero links remaining in the network after thresholding (see Methods sections 5.2 and 5.5, and Supp. section 5).

The two control parameters of the ESSES model are thus the global recovery probability γ and the network density ρ. Using the iEEG data for the modeling cohort patients, we fit γ and ρ to reproduce the iEEG-recorded seizure propagation patterns with the ESSES-modeled seizures, setting the resection area was used as the epidemic seed (see 5.5 and Supp. section 5.2 for more details) [30, 32]. The degree of similarity between the ESSES and iEEG seizures was measured with the goodness-of-fit C(ρ, γ) as defined in eq. 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 [38]. Most individual best fits (data-points in figure 1B) fell within this region, although there was large variability among the individual fit maps (in fact, we found low signal to noise ratios of approx. 1/2, see figure 1C).

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 (, C_SF − C_NSF = 0.11, rank-sum rks = 72, p = 0.04), as shown in figure 1D. A ROC classification analysis indicated a good classification (AUC = 0.79) between the SF and NSF groups based (figure 1E). At the optimal classification point (Youden criterion, figure 1F), all NSF patients were correctly identified in this case (accuracy = 0.80, precision = 0.57, sensitivity = 1.00 and F 1 = 0.71). 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).

2.2 Presurgical hypothesis of the seed regions

A key ingredient of ESSES is the definition of the epileptogenic or seed regions. Here we defined epileptogenicity or seed-propbability maps SP_i indicating the probability that each brain region i gave rise to a seizure The seed-probability maps integrated the patient-specific multimodal presurgical information (encoded in the local patient database [39]) 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 5.6 and Supp. section 4 for details). The resulting seed-probability maps for two representative patients from the modeling cohort are illustrated in figure 3B,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 in 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 map between SF and NSF patients (see Sup. Fig. S1 and S2).

• Download figure
• Open in new tab

Figure 3:

Seed-probability maps. Resection areas (left) and seed-probability maps (right) as derived from the database with presurgical information for two representative cases from the modeling cohort: patient 3 (SF, top) and 6 (NSF, bottom).

2.3 Optimal resection strategies

ESSES can derive individualized optimal (alternative) resection strategies via an optimization algorithm based on simulated annealing. The optimization algorithm was first designed using the modeling cohort data (see Methods section 5.7), and was then applied to the validation cohort in a blind setting. Here we report on the results for the validation cohort, corresponding analyses for the modeling cohort can be found in the supplementary information (Supp. section 5.3 and Supp. figures 3 and 4).

The optimization algorithm searched for resections R of increasing size S(R) that minimized the seed efficiency E_R(seed), exploiting the link between epidemic spreading dynamics and network structure. To diminish differences due to seed extent and initial efficiency, we defined the normalized seed efficiency e_R(seed) = E_R(seed) − E₀(seed), where E₀(seed) is the initial seed efficiency (figure 4A). e_R(seed) decreased with the size of the resection for all patients. At the group level, the SF group showed a significantly larger 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 the NSF group (F (19) = 3.78, p < 10⁻⁶).

• Download figure
• Open in new tab

Figure 4:

Optimal (alternative) resection strategies (validation cohort). Effect of optimal virtual resections of size S(R) on A normalized seed efficiency e_R(seed), B seizure propagation after the resection IR(R), and C normalized decrease in seizure propagation δIR(R). Blue dashed lines stand for NSF patients, and pink solid lines for SF patients. Thin lines show individual patients, and darker wide lines the group averages, with shaded areas indicating the standard deviations. D,E Group level comparison of the size of optimal resections S(R_op) (D) and their overlap with the resection area Ov(R_op, R_D) (E). The left panels show the distribution of values of each patient group, with significance results obtained with exact two-sided Wilcoxon ranksum tests. The center panels show the corresponding ROC classification analyses, where TPR and FPR stand respectively for the true positive and false positive rates. Finally, the right panels show the confusion matrices corresponding to the optimal point (Youden criterion, black asterisks in the middle panels) of the ROC curves.

The expected effect of a resection R on seizure propagation was quantified by comparing seizure propagation IR before (IR(0)) and after (IR(R)) the resection. Seizure propagation depended heavily on the seed realization such that a bi-stable regime emerged where, depending on the seed realization, the ESSES seizures either propagated macroscopically or died locally (an exemplary case is shown in Supp. figure 3). Thus, results reported here were averaged over 300 independent realizations of the seed regions and SIR dynamics. The SF group presented a significantly larger decrease in seizure propagation (F (19) = 16.65, p < 10⁻⁴³), than the NSF group. In this case the effect of increasing the resection size was not significantly different between the two groups (F (19) = 1.43, p = 0.11).

We defined the normalized decrease in seizure propagation δIR(R) = (IR₀ − IR_R)/IR₀ to compare the effect of the alternative resections between different patients (figure 4C). 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 δIR(R) as a function of 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 which may be small (spanning only a few nodes) or large. In order to quantify these differences, 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 (figure 4D-left, diff=−4.34, p = 0.03, rks = 411.5), and these presented a significantly larger overlap with the actual resection strategy Ov(R_op, RA) (figure 4E-left, diff= 0.20, p = 0.03, rks = 495.5), that 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), figure 4D,E). Both variables led to very similar classification results at the optimal classification point (Youden criterion), correctly identifying 6/8 NSF cases. In particular, we found accuracy = 0.74, precision = 0.46, sensitivity = 0.75, F 1 = 0.57 for S(R_op), and accuracy = 0.71, precision = 0.43, sensitivity = 0.75, F 1 = 0.55 for Ov(R_op, RA). The classification results for both the validation and modeling cohorts are summarized in table 1.

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, there was only a 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 according would be indicative of the need to perform a more exhaustive presurgical evaluation, and potentially imply an unavoidable non-seizure free outcome after any surgery.

Finally, we note that almost equivalent results may be obtained by considering the disconnecting resection, i.e. the smallest resection leading to disconnection of the seed, instead of the optimal resection (see Supp. section 5.4 and Supp. figure 5). 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.

2.4 Simulation of the surgical plan

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 supplementary information (Supp. section 6, Supp. figure 6). As in the previous section, 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, δIR(RA, SF) − δIR(RA, NSF) = 0.26, rks = 513, p = 0.02). 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 patients were correctly identified (accuracy = 0.79, precision = 0.55, sensitivity = 0.75, F 1 = 0.63, figure 5D). 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.

• Download figure
• Open in new tab

Figure 5:

Simulation of the planned resection strategy (validation cohort). A The top panel shows seizure propagation IR before (left point cloud for each patient) and after (right point clouds) the resection, for 300 iterations of the seed regions, for each patient. The bottom panel shows the average relative decrease in seizure propagation δIR(RA), with errorbars given by the standard deviation over seed iterations. B Comparison of the relative decrease in seizure propagation δIR(RA) between the SF and NSF groups. Each point corresponds to one patient. C ROC curve of the group classification based on δIR(RA). TPR and FPR indicate respectively the true positive (NSF cases classified as NSF) and false positive (SF cases classified as NSF) rates. D Classification results for the optimal point (black asterisk in panel B) of the ROC curve according to the Youden criterion.

2.5 Prediction of surgical outcome

The classification analyses in the previous sections were informed by each patient’s surgical outcome. In a prospec-tive 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 crossvalidation. 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. The results from these analyses were slightly worse than in the previous section, 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 7B). 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 probabily of complete seizure freedom after the surgery.

• Download figure
• Open in new tab

Figure 6:

Prediction of surgical outcome: pseudo-prospective study. A Prediction results using each of the three model-based biomarkers of surgical outcome: the size of optimal resections S(R_op), the overlap between optimal resections and the resection area, Ov(R_op, RA), and the decrease in seizure propagation due to simulation of the planned resection strategy, δIR(RA). NSF (SF) cases are shown by black (white) rectangles. The bottom row shows the fraction of models (0 − 3 out of 3) with a positive (i.e. NSF) classification, for each patient. Surgical outcome is shown in the top row. NSF cases are highlighted by a red arrow and by red labels. B Relative number of cases identified as NSF by n models, n = 0, 1, 2, 3, respectively for the SF (blue, left-side bars, N = 26) and NSF (red, right-side triangles, N = 8) groups.

• Download figure
• Open in new tab

Figure 7:

Prediction of surgical outcome using a machine learning algorithm (RUSBoost) with leave-one-out cross validation. As input variables we used the normalized decrease in seizure propagation after virtual resection of the RA, δIR(RA), the size of optimal resections S(R_op) and the overlap of optimal and clinical resections Ov(R_op, RA). Panels A,B show the confusion matrix and predictor importance for the validation cohort (N = 34, 8 NSF), and panels C,D are for the combined cohort (N = 49, 12 NSF).

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 (figure 7A,B). The input variables for the classification algorithm where δ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 F 1 = 0.36, even though the considered algorithm (RUSboost) was designed to correct for class imbalance. However, the minority class in our case contained only 8 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 (figure 7C,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, F 1 = 0.51). Overall, the machine learning model was not able to improve upon the results found using the individual variables (see table 2). 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).

3.1 ESSES modeling framework

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 model. 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 Ref. [24] and improved in Ref. [30], whereas the seizure propagation and parameter fitting model as used here was mainly defined in Ref. [32]. 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 [30, 32], and in contrast with other works [21, 24, 29, 42]. 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 [30] with much lower computational cost than DTI networks, whilst also being more sensitive to long range connections, in particular inter-hemispheric ones, that may often be missed by DTI [43].

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 [31], and have been used to describe other neuro-physiological processes before, such as the spreading of pathological proteins on brain networks [44, 45] or the relation between brain structure and function [46]. Epidemic spreading models have been extensively studied on different network substrates [31] and are supported by a well-grounded mathematical and computational framework that we can use to our advantage in the context of epilepsy 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 [12, 20], with the spreading threshold heavily influenced by the existence of hubs [31]. 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 [30, 32]. As the working point of the ESSES model 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 multi-modal presurgical information). As we showed in our previous study [32] and in the supplementary information here (Supp. 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 by 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 [29], or these may have presented seizure dynamotypes [47] 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 and intrinsic difference between the presurgical data of the two groups, and not only to a sub-optimal 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. This set-up thus emulated the clinical situation prior to the surgery, where a surgical strategy has been devised based on the available information 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 [48]. The virtual resection optimization algorithm was originally validated in our previous studies [24, 30]. 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 [24, 30]. In particular, Nissen et al. [24] found that removing connections to the network hubs was the most efficient way to decrease seizure propagation, whereas Millán et al. [30] 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 [30]. However, the intrinsic noise in the seed definition, in the seed-probability maps, and in the actual origin and propagation patterns of iEEG-recorded seizure 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 [30]). 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 [23, 42]. 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).

3.2 Modeling considerations and limitations

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, i.e. 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. 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 bio-physiological processes involved in seizure generation and propagation, but aims to focus only on the most relevant features of seizure propagation [24, 29, 30, 32]. In order to validate ESSES as a framework to simulate seizures, we compared the modeled seizures with those recorded via iEEG. This required however of a simplified representation of the iEEG data. In particular, as there was no intrinsic time-scale 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 [48]. 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 [24, 30].

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 [49]. ESSES was validated in previous studies [24, 30, 32], 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 [32]. 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 [39]. Seed regions were conse-quently widespread over the network. This also led to large variability in the analyses for each patient (see for instance figures 4A-C and 5A), as the results of each simulation 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 database 34-region description.

Finally, it is a limitation of this study the small size of the non-seizure-free group, with only 4 cases in the modeling cohort and 8 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. Future studies involving more that one center have the potential to at least diminish this limitation.

5.1 Patient groups

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 post-surgical 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 Supp. tables 1 and 2 for details). Surgical outcome was classified according to the Engel classification at least one year after the surgery [50]. Patients with Engel class 1A were labelled as seizure-free (SF), and patients with any other class were labelled as non-seizure-free (NSF). The modeling cohort consisted of 15 patients (4 NSF, 11 females) who had also undergone an iEEG (invasive electroencephalography) study, including post-implantation CT-scans. This same cohort was already included in Ref. [32], and partially in [30]. 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).

5.2 Individualized Brain Networks

Seizure propagation was modeled on the patient-specific brain networks, as derived from MEG data, for both cohorts (see figure 2). 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 [30, 32] and are detailed in the supplementary information (Supp. section 2). Functional networks were generated considering each of the 246 ROIs of the Brainnetome (BNA) atlas [51] 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) [52–55], 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, it. We validated the relationship between AEC-MEG and structural networks in a previous study [30] 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 [56–58], i.e. w_ij ∝ exp(−αd_ij). Recent studies have corroborated that the EDR reproduces human DTI data well [59–61], and in Ref. [30] we validated that AEC-MEG networks were strongly correlated (R² = 0.50) with the corresponding EDR networks, therefore showing that AEC-MEG encodes structural information. Moreover, AEC-MEG networks also include long-range connections that may promote seizure propagation, but that may be missing from structural (i.e. DTI) networks [62, 63]. 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 re-scaled between 0 (perfect anti-correlation) and 1 (perfect correlation), with 0.5 indicating no coupling [64]. 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.

5.3 Resection Area

The resection area (RA) was determined from the three-month post-operative MRI. For the modeling cohort the resection areas were obtained as part of two previous studies [30, 32]. 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 pre-processing step, as described in figure 2. Cases 9 and 20 of the validation cohort underwent the post-operative MRI on a different MRI scanner at their resection center, respectively one day and three 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 post-resection MRIs were co-registered to the pre-operative 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.

5.4 iEEG Seizure Propagation Pattern

Patients in the modeling cohort underwent invasive EEG recordings using stereo-tactic electrode implantation as described in [32]. 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 [32] and in the Supp. section 3.

5.5 Seizure Propagation Model

ESSES was based on our previous studies [24, 30, 32] where we showed that simple epidemic spreading models could reproduce the spatio-temporal 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 [31], 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 neighbour 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 the ESSES model are thus the network density ρ, and the recovery probability γ. Depending on the network structure, the epidemics can show different spatio-temporal spreading profiles described by the probability p_i(t) that each ROI i becomes infected at step t.

ρ and γ were fitted to the iEEG seizure-propagation patterns at the group level. The resection area was set as the seed of epidemic spreading, and an ESSES seizure propagation pattern was built that described the set of infected and non-infected ROIs during the SIR-simulated seizures, as well as the order in which infected ROIs became infected. In order to take into account the stochastic nature of the SIR dynamics, the participation of each ROI was weighted by the fraction of realizations in which it was involved in the simulated seizure (since different ROIs became infected in different realizations). The goodness-of-fit of the model, C(ρ, γ) [32], quantified how similar the ESSES and iEEG patterns were. It took into account two factors: the weighted correlation between activation orders of ROIs that were active in both patterns, C_w, and the overlap between the active and inactive ROI sets of both patterns, P_overlap, i.e. The details of this definition can be found in Supp. section 5.2.

We estimated C for a range of values ρ ∈ {0.025, 0.05, 0.10, 0.20, 0.30} and γ ∈ 10⁻⁴, 10⁻³, 10⁻², 10⁻¹, 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 Supp. section 5.2 and Supp. figure 2) 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 algoritihm) 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.

5.6 Presurgical hypothesis of the seed regions

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 [24, 30, 32], which can only be known after the surgery, or from the iEEG data [30, 32]. 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 [39]).

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 6 frontal regions (fronto-orbital, frontal-basal, frontal-parasagitaal, frontal-periventricular, frontal-lateral, frontal-operculum), 6 temporal regions (hippocampus, amygdala, uncus, anterior-neocortical, posterior-neocortical, gyrus-parahippocampalis), 2 insular regions (anterior and posterior insula), 1 central, 1 parietal and 1 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.

For each region i and modality m, the database indicates the presence (1) or absence (0) of abnormalities, from which we derived binary abnormality maps a_i,m = 0, 1. The overall abnormality map A_i was obtained by aggregating over all modalities available for each patient. Not all modalities are equally relevant to establish the probability that a region is involved in epileptogenic activity: EEG is the least focal, whereas iEEG provides the most localized information, and its results also integrate information from the other modalities (as these affect where the iEEG electrodes are placed). In order to gauge these differences, we weighted each modality m by a relevance factor ω_m, with ω = 1 for EEG, 2 for MRI, MEG, PET and ISPECT, and 4 for iEEG. Thus, the overall abnormality map was defined as where the normalization factor n is defined as

A clinician (ECWvS) defined a unique projection of the regions in the database unto 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 supplementary material. 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 co-occurrence of abnormalities in different modalities is a strong indicator of the epileptogenic zone, we included a re-scaling 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.

5.7 Virtual Resections

We conducted virtual resections of sets of nodes R by disconnecting them from the network, by setting to 0 all their connections. The effect of each resection was characterized by the normalized decrease in seizure propagation δIR(R) in the resected network (R) with respect to the original (0) one: where IR is the fraction of nodes that became infected at any point during the modeled seizure, namely, That is, IR takes into account all nodes that became infected during the simulated seizure, regardless of whether they eventually recovered or not. This characterizes the total extent of the simulated seizure.

We performed two virtual resection studies, as detailed in figure 2. Firstly, we performed an Optimization of alternative resections analysis. We derived optimal virtual resections R of increasing sizes S(R) (defined as the number of resected nodes) with an optimization algorithm based on simulated annealing [65] and derived in our previous studies [24, 30]. The optimization method took advantage of the relationship between SIR spreading and network structure to use a structural metric –the seed efficiency– as a proxy for the actual effect of the resection on seizure propagation δIR(R). Thus, for each resection size S(R), the simulated annealing algorithm searched for the resection R that minimized the seed efficiency E_R(seed) [37, 66, 67]. All nodes were considered as possible targets of the resection. To compare between different patients we defined the normalized seed efficiency where E₀(seed) is the seed efficiency in the original (un-resected) network. The actual effect of each resection was quantified by the seizure propagation level after the resection, IR(R), and the normalized decrease in seizure propagation δIR(R). We defined the optimal ESSES resection R_op, as the smallest resection leading to (at least) a 90% decrease in (modeled) seizure propagation. This resection was characterized by its size S(R_op) and overlap with the resection area Ov(RA, R_op). We also defined the disconnecting resection R_D as the smallest resection that lead to seed disconnection (see Supp. section 5.4 and Supp. figure 5).

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 the right side of figure 2, for the validation cohort we first performed the Optimization of alternative resections (figure 2) 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 the ESSES optimal resections with the resection strategy and ii) measure the effect of the planed resection in decreasing (modeled) seizure propagation. Finally the one-year postoperative outcome was also de-blinded and used for the statistical analyses.

5.8 Statistics

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 (non-seizure-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= (T P + T N)/(T P + FP + FN + T N), precision = T P/(T P + FP), sensitivity= T P/(T P + FN), F 1 statistic (harmonic mean between precision and sensitivity) = 2T P/(2T P + 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 crossvalidation-type setting. The predictive model compounded the prediction results from these N = 34 models. We measured its accuracy, precision, sensitivity and F 1 statistic.

Finally, for the Machine Learning analyses we used the Adaboost1 algorithm (Matlab 2018) using leave-one-out-cross-validation and default hyperparameters. The numbers of learners in each model was equal to the group size minus one (33 for the validation cohort and 48 for the combined cohort), and results were averaged over 10 iterations of the Adaboost1 algorithm for each classification model.

5.9 Data availability

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 data-sets should be directed to the corresponding author. All user-developed codes are publicly available on Github https://github.com/anapmillan/computer_model_for_epilepsy.