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

Abstract Epilepsy surgery is the treatment of choice for drug-resistant epilepsy patients, but up to 50% of patients continue to have seizures one year after the resection. In order to aid presurgical planning and predict postsurgical outcome on a patient-by-patient basis, we developed a framework of individualized computational models that combines epidemic spreading with patient-specific connectivity and epileptogeneity maps: the Epidemic Spreading Seizure and Epilepsy Surgery framework (ESSES). ESSES parameters were fitted in a retrospective study (N = 15) to reproduce invasive electroencephalography (iEEG)-recorded seizures. ESSES reproduced the iEEG-recorded seizures, and significantly better so for patients with good (seizure-free, SF) than bad (nonseizure-free, NSF) outcome. We illustrate here the clinical applicability of ESSES with a pseudo-prospective study (N = 34) with a blind setting (to the resection strategy and surgical outcome) that emulated presurgical conditions. By setting the model parameters in the retrospective study, ESSES could be applied also to patients without iEEG data. ESSES could predict the chances of good outcome after any resection by finding patient-specific model-based optimal resection strategies, which we found to be smaller for SF than NSF patients, suggesting an intrinsic difference in the network organization or presurgical evaluation results of NSF patients. The actual surgical plan overlapped more with the model-based optimal resection, and had a larger effect in decreasing modeled seizure propagation, for SF patients than for NSF patients. Overall, ESSES could correctly predict 75% of NSF and 80.8% of SF cases pseudo-prospectively. Our results show that individualised computational models may inform surgical planning by suggesting alternative resections and providing information on the likelihood of a good outcome after a proposed resection. This is the first time that such a model is validated with a fully independent cohort and without the need for iEEG recordings.


Patient data
Basic statistics of the patient groups are summarized in Supp.tables 1 (modeling cohort) and 2 (validation cohort).

Individualized MEG brain networks
Seizure propagation was modeled on the patient-specific brain networks, as derived from MEG data, for both cohorts.One to three eyes-closed resting-state (supine position) MEG recordings of 10 to 15 minutes each were acquired for each patient in a magnetically shielded room (Vacuumschmelze GmbH, Hanau, Germany) during routine presurgical clinical practice (whole-head MEG system (Elekta Neuromag Oy, Helsinki, Finland) with 306 channels consisting of 102 magnetometers and 204 gradiometers).The data were sampled at 1250 Hz, and online filtered with an antialiasing filter at 410 Hz and a high-pass filter of 0.1 Hz.The head's position relative to the MEG sensors was determined using the signals from 4 to 5 head-localization coils that were recorded continuously.The positions of the head-localization coils and the outline of the scalp (roughly 500 points) were measured with a 3D digitizer (Fastrak, Polhemus, Colchester, VT, USA).The temporal extension of Signal Space Separation (tSSS) [1,2] was used to remove artifacts using Maxfilter software (Elekta Neuromag, Oy; version 2.1), and the MEG data were filtered in the broadband (0.5 − 48.0 Hz).For a detailed description and parameter settings see Hillebrand et al. [3].
Pre-operative MRI scans (8 channel phased-array head coil on a 3T whole-body MR scanner, Discovery MR750, GE Healthcare, Milwaukee, Wisconsin, USA) were used for co-registration with the MEG data.Anatomical 3D T1-weighted images were obtained with a fast spoiled gradient-recalled echo sequence, and were interpolated to 1 mm isotropic resolution during reconstruction.The points on the scalp surface were used for co-registration of the MEG scans with the anatomical MRI of the patient through surface-matching.A single sphere was fitted to the outline of the scalp and used as a volume conductor model for the beamforming approach.
Neuronal activity was reconstructed using an atlas-based beamforming approach, modified from Hillebrand et al. [4], to reconstruct the time-series of neuronal activation of the ROI centroids [5].We considered the 246 ROIs of the Brainetome (BNA) atlas [6], the centroids of which were inversely transformed to the co-registered MRI of the  patient.Then, a scalar beamformer (Elekta Neuromag Oy; beamformer; version 2.2.10) was applied to reconstruct each centroid's time-series, as detailed elsewhere [5].The time-series of each centroid were visually inspected for epileptiform activity and artifacts.On average, 59.08 (range: 48 − 83) 16384-sample interictal epochs of sufficient quality were selected for each patient.The epochs were further analyzed in Brainwave (version 0.9.151.5 [7]) and were down-sampled to 312 Hz, and filtered in the broadband (0.5 -48 Hz).

iEEG Propagation Pattern
Patients in the modeling cohort underwent invasive EEG recordings using stereotactic electrode implantation as described also in [8].The number and location of the intracerebral electrodes (Ad-Tech, Medical Instrument Corporation, USA, 10-15 contacts, 1.12 mm electrode diameter, 5 mm intercontact spacing; and DIXIE, 10-19 contacts, 0.8 mm electrode diameter, 2 mm contact length, 1.5 insulator length, 16 − 80.5 insulator spacer length) were planned individually for each patient by the clinical team, and were based on the location of the hypothesized SOZ and seizure propagation pattern.The number of electrodes per patient varied between 9 and 15 (average = 12.1 ± 1.8) and the total number of contacts between 194 and 99 (average = 121 ± 24).The locations of the electrode contact points were obtained for each patient from the post-implantation CT scan (containing the iEEG electrodes) that was co-registered to the preoperative MRI scan using FSL FLIRT (version 4.1.6)12 parameter afine transformation.Each contact point was assigned the location of the nearest BNA ROI centroid.Given that BNA ROIs are in general larger than the separation between contact points, different contact points could be assigned to the same BNA ROI.
An iEEG seizure pattern was derived for each patient by a clinician (ECWvS).First, a representative seizure was chosen for each patient.Then, the onset time of ictal activity was identified for each iEEG channel, and the channels were sorted according to their onset times.In case two channels became active at the same time, they were assigned the same activation order.This activation pattern was then translated into the BNA space, so that the each sampled ROI i was assigned an activation step.If channels with different onset times belonged to the same BNA ROI, the ROI was assigned the earliest of the possible times.This constituted the iEEG seizure pattern.

Seed-probability maps
In order to characterize the seed-probability maps, in Supp.figure 1 we show the fraction of seed-probability accounted for by RA nodes (panels A and E), the fraction of nodes with non-zero seed-probability that belonged to the RA  (panels B and F) and the comparison of the RA seed-probability (panels C and G) and node (panels D and G) fractions between the SF and NSF groups, respectively for the modeling and validation cohort.

SIR Dynamics
The discrete-time SIR model [9,10] was used to simulate the propagation of ictal activity from a set of seed nodes (or brain regions) that were set to be infected (or in the ictal state) at the beginning of the simulation.The SIR dynamics considered here are characterized 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, and due to the limited amount of data available (one iEEG-recorded seizure per patient, and N = 15 patients) we considered here a global recovery rate γ i = γ and spreading probabilities given by the MEG network structure: As discussed in our previous studies [8,11], the iEEG seizure patterns only accounted for the propagation of ictal activity and not the recovery to the healthy state, thus considering local recovery rates is beyond the scope of this study.In this manner we avoided introducing arbitrary time scales (i.e. to account for the duration of the ictal state) into the seizure propagation model.
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 (and only one) new node became infected, and they were iterated N R = 10 4 times per set of (ρ, γ) parameters.

Parameter fit: the ESSES model
The two control parameters of ESSES are the density of connections in the network, ρ, indirectly setting the spreading probabilities, and the global recovery rate γ.ρ and γ were fitted in the retrospective study (modeling cohort) by comparing the spatio-temporal propagation pattern of the ESSES-modeled seizures to the patient's iEEG seizure pattern (constructed as described above), when the resection area was used as the seed for epidemic spreading.For each set of ESSES model parameters, we measured the spatio-temporal seizure propagation profile described by the probability p i (t) that each ROI i becomes infected at step t.Then, mean infection time t i of each ROI i is defined as where T is the maximum integration time.t i was then sub-sampled to the ROIs sampled by the iEEG electrodes, and each of these ROIs was assigned an infection order according to t i .This constituted the ESSES seizure pattern.
The goodness-of-fit C(ρ, γ), quantified how similar the ESSES and iEEG seizure-propagation patterns were, and it was measured following the same procedure as in Ref. [8].It took into account two factors: the weighted correlation between infection orders of ROIs that were infected in both patterns, C w , and the overlap between the infected (or active) and susceptible (or inactive) ROI sets of both patterns, P overlap , i.e.
In order to take into account the stochastic nature of the SIR dynamics, C w was weighted by the fraction of realizations that each ROI i got infected during a modeled seizure, P IR (i).Similarly, P overlap was weighted by P IR (i) as: where N iEEG is the number of ROIs sampled by the iEEG electrodes (on average = 43.6 ± 7.9), and S and H are respectively the sets of active (in the seizure state) and inactive (in the healthy state) ROIs in the iEEG pattern.The total correlation C equals 1 in case of exactly equal activation patterns, 0 in the case of null-overlap or correlation, and −1 in the case of complete anti-correlation of activation times (but equal seizure areas).We note however that C decays from 1 faster than a simple correlation metric, since it takes into consideration not only the activation times, but also the activation areas.
In order to fit ρ and γ, we measured 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 the resection area as the seed of seizure propagation [8,11].For each parameter configuration we performed 10 iterations (each over N R = 10 4 iterations of the SIR dynamics), to obtain average C values and their fluctuation for each patient.We then found the parameter set that maximized C for each patient (individual model fit) and at the group level (population model fit).The population model fit was used for all subsequent analyses in the main text (for both the modeling and validation cohorts).As in our previous study [8], we found that the population model fit led to a better classification between the SF and NSF groups, likely due to a reduced influence of noise.
The results for the individual model fit are summarized in Supp.figure 2 (on average, C = 0.12 with standard deviation std C = 0.07).The model provided a (not significantly) better fit for the SF (C SF = 0.14 ± 0.08) than NSF (C NSF = 0.08 ± 0.02) patient groups (C SF − C NSF = 0.05, ranksum rks= 99, p = 0.18, exact two-sided Wilcoxon ranksum test), as shown in Supp.figure 2A.A ROC classification analysis based on the goodness-of-fit returned a good classification result with an area under the curve (AU C) of AU C = 0.750 (Supp.figure 2B).The confusion matrix corresponding to the optimal point of the ROC curve, according to the Youden criterion (to control for the class imbalance), is shown in 2C.We found an accuracy of 0.80, precision = 0.60, sensitivity = 0.75 and F 1 = 0.67.There were no significant differences in the fit parameters ρ and γ between the groups.
Supp. figure 3 shows the results of the model fit at the group level for the modeling cohort.The details of the statistical comparison of the model fit results between the seizure-free and non-seizure-free groups are shown in Supp.table 3.

Alternative resection strategies
We used an optimization method based on simulated annealing to find alternative resection strategies R of increasing size S(R).The optimization algorithm made use of the link between network structure and spreading dynamics to find resections that led to a minimum seizure propagation by minimizing the efficiency of the seed, E seed .The seed in this case was defined as the set of regions with (significantly) non-zero seed-probability, and the contribution of each ROI was weighted by its seed-probability (see methods section for more details).At the population level, E seed (R) decreased with the size of the resection for all patients (see an exemplary case in Supp.figure 4A-top, and the remaining cases in Supp.figure 5).To diminish differences due to seed extent and initial efficiency, we computed the normalized seed efficiency e R (seed) = E 0 (seed) − E R (seed) (Supp.figure 4B-top for the population average).At the group level, the SF group showed a significantly larger e R (seed) than the NSF group (F (19) = 14.80, p = 10 −32 ), although the effect of increasing the resection size was not significantly different (F (19) = 0.74, p = 0.8).
The actual effect of each resection on seizure propagation in the model was quantified using the SIR dynamics and 300 independent seed realizations.An exemplary case is shown in Supp.figure 4A-bottom, and the group average can be seen in Supp.figure 4B-middle.Seizure propagation depended heavily on the seed realization.A bi-stable regime emerged where, depending on the seed realization, the modelled seizures either propagated macroscopically or died locally.For each resection size we measured the normalized decrease in seizure propagation δIR(R).At the group level, the SF group had a significantly smaller IR (F (19) = 7.86, p = 10 −16 , Supp.figure 4B-middle), but there was no significant difference in the effect of increasing the resection size (F (19) = 0.29, p ≈ 1).
We defined the normalized decrease in seizure propagation to compare between different patients as δIR(R) = (IR 0 − IR R )/IR 0 (Supp.figure 4B-bottom).At a group level, the SF group presented significantly larger δIR(R) for S(R) > 1 (F (19) = 8.57, p = 2 • 10 −18 ), and again there was no significant difference in the effect of increasing the resection size between the groups (F (19) = 0.27, p ≈ 1).
In order to compare the two groups systematically, we defined the optimal resection R op as the one leading to a 90% decrease in seizure propagation, δIR = 0.90, and the disconnecting resection R D as the smallest resection leading to seed disconnection, e R D = 0. We found that the SF group had smaller optimal resections (panel C-top) than the NSF group, but the difference was not significant (diff= −5.73, rks = 81 p = 0.37).The overlap between the optimal and actual resection was also larger, but not significantly so, for the SF group (panel D-top, diff= 0.23, rks = 98 p = 0.20).The disconnecting resection (panel E-top) was smaller, but not significantly so, for the SF than for the NSF group (diff= −6.48, rks = 81 p = 0.37).Either of these three variables could classify the two groups with good AUC results of 0.66, 0.71 and 0.66, respectively (mid subpanels of panels C, D and E).The confusion matrices corresponding to the optimal points (Youden criterion) are also shown in Supp.figure 4C,D,E.The classifications resulted in accuracy = 0.80, precision = 0.67, sensitivity = 0.50 and F 1 = 0.57 for the size of optimal resections, accuracy = 0.67, precision = 0.44, sensitivity = 1.00 and F 1 = 0.62 for the overlap of the optimal resections with the resection area, and accuracy = 0.80, precision = 0.67, sensitivity = 0.50 and F 1 = 0.57 for the size of disconnecting resections.
Supp. Figure 4: Optimization of virtual resections (modeling cohort).A Effect of resections of size S(R) for an exemplary patient (case 3) as shown by the seed efficiency E(seed) R (top) and total seizure propagation IR R after the resection (bottom).Light lines in the top panel correspond to individual realizations of the optimization algorithm, and the solid wide line to the optimal resection of each size.This resection was used to quantify the effect on seizure propagation in the bottom panel.Light lines correspond each to a seed realization, whereas the solid wide line corresponds to the average over seed realizations.Brain plots on the right show E(seed) R (left) and IR R (right) for three exemplary resection sizes (S(R) = 1, 3, S D ), where S D is the size of the disconnecting resection.B Group level analysis of the normalized seed efficiency e(seed), seizure propagation IR R and normalized decrease in seizure propagation δIR(R) after (optimal) resections of size S(R).Blue dashed lines and pink solid lines stand respectively for NSF and SF patients.Wide lines indicate the group averages with uncertainties shown by the shaded areas, as given by the standard deviation.The dotted black line in the bottom panel indicates the 90% threshold used to define the optimal resection.C, D, E Group level comparison of the size of optimal resections (C), its overlap with the actual resection area (D) and the size of disconnecting resections (E).The top panels illustrate the violin plots for the SF and NSF groups, with significance results given by exact two-sided Wilcoxon ranksum tests.Data-points are coded with the same colour and marker type as in Supp.figure 1A.The middle panels show the corresponding ROC classification analyses, where FPR and TPR stand respectively for the false positive and true positive rates.The bottom panels show the confusion matrices corresponding to the optimal classification points (Youden criterion).Supp. Figure 5: Optimization of virtual resection strategies for SF (panel A) and NSF (panel B) patients, showing the seed efficiency E R (seed) after virtual resections of size S(R).The light lines indicate individual realizations of the simulated annealing algorithm, whereas the thick lines indicate the best iteration (i.e. the one leading to the minimum E R (seed).Each patient is shown by a different colour as indicated in the legends, with the same color code as Supp.figure 1A.

Disconnecting resection
The disconnecting resection R D was defined as the smallest resection leading to seed disconnection, i.e. to E R (seed) = 0.It is computationally much faster to obtain than R op given that it does not require the simulation of the SIR spreading dynamics.As shown in Supp.figure 6E, the size of R D is almost identical (correlation coefficient = 0.99) to the size of R op .
We also show in Supp.figure 6 the correlation among the 4 model biomarkers (size of optimal and disconnecting resections S(R op ) and S(R D ), overlap between the optimal and actual resection Ov(RA, R op ) and (model-based) effect of the actual resection δIR(RA).Only S(R op ) and S(R D ) were strongly correlated.Supp. Figure 6: Pearson's correlation coefficient between each pair of model-derived variables: δIR(RA), S(R op ), Ov(RA, R op ), S(R D ).Black (dark) data-points stand for patients in the modeling cohorts, whereas orange (light) markers stand for patients in the validation cohort (thus the combined cohort includes all data-points).The orange (light) lines indicate the linear correlation for the validation cohort, and the brown (dark) ones stand for the combined cohort V + M .The correlation coefficients are indicated in the legends.

Modeling pipeline
The complete modeling pipeline is shown in Supp.figure 7. The study was performed in two steps.First, a retrospective study (modeling cohort) took place to fit the parameters of the different ESSES modules.For this cohort, the MEG and resection-area data had already been processed in a previous study [8].The second step has been described in the main text (figure 2), and involves the validation cohort.Supp. Figure 7: Methodological pipeline.The study was performed in two steps, with first a retrospective study to set-up the model parameters, followed by a pseudo-prospective application of the framework.The elements in this figure are color-coded according to the analysis step: data pre-processing (blue), model fit results (green, only modeling cohort), optimization of alternative resections (pink) and simulation of the resection plan (yellow).The final step consisted of the statistical analyses (group differences, patient classification and outcome prediction).

Simulation of the resection strategy: Modeling cohort
We made use of the population model and the seed-probability maps to simulate the effect of the clinical resection of each patient.The resection was simulated by disconnecting the resection area from the network, and the seed regions were derived probabilistically from the seed-probability maps.The effect of the resection was quantified as the relative decrease in seizure propagation due to the resection δIR(RA) = (IR 0 − IR R (RA)) /IR 0 , where IR 0 and IR R (RA) stand respectively for the propagation before and after the resection (Supp.figure 8A).The seizure propagation and the effect of the resection depended strongly on the seed realization, so the same set of seeds were used before and after the resection to estimate δIR(RA), and we considered 300 independent seed realizations to diminish noise effects.
The effect of the surgery was larger for the SF than for the NSF group (Supp.figure 8B), although the group difference was not significant (δIR(RA, SF ) − δIR(RA, N SF ) = 0.14, rks = 98, p = 0.22).A ROC analysis using δIR(RA) to classify the two groups found an AUC= 0.73 (Supp.figure 8C).The confusion matrix corresponding to the optimal classification point (Youden criterion, Supp.figure 8D) shows an accuracy of 0.73, precision = 0.50, sensitivity = 1.00 and F 1 = 0.67.Again, in this case the classification analysis was able to identify all NSF patients correctly.
Supp. Figure 8: Simulation of the resective surgery (modelling cohort).A Effect of the virtual resections of the resection area in the model.The top panel shows the total seizure propagation IR before and after the resection for each patient (respectively left and right violin plots for each x-tick), for 300 independent realizations of the seed regions.The white markers (respectively circles for SF and triangles for NSF cases) show the average over the seed realizations.NSF cases are also highlighted by red labels on the x-axis.The bottom panel shows the relative decrease in seizure propagation δIR(RA) for each patient, as given by the average (data-points) and standard deviation (error bars).The colorcode is as in Supp.

Figure 1 :
Characterization of seed-probability maps, modeling (top) and validation (bottom) cohorts.A and E Fraction of the total seed-probability in RA nodes, f seed (RA), for each patient.Filled markers indicate that the seed-probability in the RA nodes was larger than expected by chance (this occurred for all patients in the modeling cohort).B and F Fraction of nodes with non-zero seed-probability that belonged to the RA, n seed (RA).In A, B, E and F, SF (NSF) patients are indicated by circular (triangular) markers.The colorcode indicates the patient, also given by the x-axis.Filled markers indicate that there were more RA nodes with non-zero seed-probability than expected by chance.C,D and G,H Comparison of the RA seed-probability (C,G) and node (D,H) fractions between the SF and NSF groups using a two-sided Wilcoxon ranksum test (panel C: diff=0.035,p = 0.85, ranksum rks = 90; panel D: diff=0.037,p = 0.92, rks = 89, panel G: diff=0.17,p = 4 • 10 −5 , rks = 542; panel H: diff=0.081,p = 0.16, rks = 489.5).

Figure 2 :
Seizure propagation model: model fit results.Group comparison of the goodness-of-fit using the individual model fit.A Violin plots of the goodness-of-fit distributions for each group.The significance analysis was performed with an exact two-sided Wilcoxon rank test.B ROC (receiving-operator-characteristic) curve analyses, where TPR and FPR respectively indicate the true positive (NSF patient classified as NSF) and false positive (SF patient classified as NSF) rates.C Confusion matrix corresponding to the optimal operating point of the ROC curve in panel B (shown by a black asterisk).Cohort SF NSF diff rks p C Modeling 0.08 ± 0.08 −0.02 ± 0.08 0.11 72 0.04 Supp.Table 3: Summary of statistical comparisons: difference in goodness-of-fit between SF and NSF groups (modeling cohort).Supp. Figure 3: ESSES parameter fitting (modeling cohort).A C std (ρ, γ) map displaying the standard deviation of the model fit C for the modeling cohort.B Violin plots of the goodness-of-fit distributions for the SF and NSF groups.C ROC (receiving-operator-characteristic) curve analysis, where TPR and FPR respectively indicate the true positive (NSF patient classified as NSF) and false positive (SF patient classified as NSF) rates.The optimal classification point (Youden criterion) is shown by a black asterisk, and the corresponding confusion matrix is shown in panel D. The confusion matrix indicates the number of SF and NSF cases that were correctly (diagonal elements) and incorrectly (off-diagonal elements) classified.
figure 1A.B, C, D Virtual resections of the resection area had a (not significantly) larger effect for the SF than NSF group, as shown by the violin plots of δIR(RA) in panel B. A ROC classification analysis of the two groups based on δIR(RA) showed an AUC= 0.73 (panel C).TPR and FPR respectively stand for the true positive (NSF patients classified as NSF) and false positive (SF patients classified as NSF) rates.The confusion matrix corresponding to the optimal operating point (Youden criterion) is shown in panel D. For this threshold we found accuracy = 0.73, precision = 0.50, sensitivity = 1.00, and F 1 = 0.67.

Table 1 :
Case Sex Resection Area S RA Patient data (modellng cohort).