Epidemic models characterize seizure propagation and the effects of epilepsy surgery in individualized brain networks based on MEG and invasive EEG recordings

Epilepsy surgery is the treatment of choice for drug-resistant epilepsy patients. However, seizure-freedom is currently achieved in only 2/3 of the patients after surgery. In this study we have developed an individualized computational model based on MEG brain networks to explore seizure propagation and the efficacy of different virtual resections. Eventually, the goal is to obtain individualized models to optimize resection strategy and outcome. We have modelled seizure propagation as an epidemic process using the susceptible-infected (SI) model on individual brain networks derived from presurgical MEG. We included 10 patients who had received epilepsy surgery and for whom the surgery outcome at least one year after surgery was known. The model parameters were tuned in in order to reproduce the patient-specific seizure propagation patterns as recorded with invasive EEG. We defined a personalized search algorithm that combined structural and dynamical information to find resections that maximally decreased seizure propagation for a given resection size. The optimal resection for each patient was defined as the smallest resection leading to at least a 90% reduction in seizure propagation. The individualized model reproduced the basic aspects of seizure propagation for 9 out of 10 patients when using the resection area as the origin of epidemic spreading, and for 10 out of 10 patients with an alternative definition of the seed region. We found that, for 7 patients, the optimal resection was smaller than the resection area, and for 4 patients we also found that a resection smaller than the resection area could lead to a 100% decrease in propagation. Moreover, for two cases these alternative resections included nodes outside the resection area. Epidemic spreading models fitted with patient specific data can capture the fundamental aspects of clinically observed seizure propagation, and can be used to test virtual resections in silico. Combined with optimization algorithms, smaller or alternative resection strategies, that are individually targeted for each patient, can be determined with the ultimate goal to improve surgery outcome. MEG-based networks can provide a good approximation of structural connectivity for computational models of seizure propagation, and facilitate their clinical use.


Optimization of virtual resections
For a given resection size S R (number of resected nodes), the optimal virtual resection R * (S R ) is the one that creates the largest reduction in seizure propagation, ∆ * R (S R ). Because characterizing all possible virtual resections is a combinatorial problem that soon looses tractability, alternative approaches are necessary to find the optimal solution. Therefore, in order to alleviate the exploration of the parameter space, we defined a random search method that combined structural and dynamical information, and made use of a simulated annealing algorithm. The method considers a surrogate topological measure for the SI dynamics, which is used to fasten up the exploration of the space of virtual resections, and the solutions are subsequently refined using the actual SI dynamics. The optimization method is described in detail below, and a schematic view is given in figure S2.

Surrogate model
In the SI model, the probability that the infection has reached a certain node i at a given time t decays in an approximate linear manner with d seed i [7], where d seed i is the average distance from the considered node to the seed (see figure S2). That is, Figure S1: Validation of AEC-MEG networks as surrogates for structural connectivity. a) Average Correlation between AEC-MEG and EDR networks for a range of exponent α values. The shaded areas indicate the standard deviation between the patients. The maximum correlation found is C = 0.71 for α = 0.052. b) Illustration of the correlation analysis for the optimal α. To calculate C, the matrices are linearized by stacking one column after another, and the standard Pearson's Correlation coefficient is computed. c) Average MEG-AEC network for the patient cohort. d) Average EDR network for the patient cohort corresponding to the optimal α value.
where d ij is the topological distance between the pair of nodes (i, j) and S seed is the number of nodes in the seed region. d ij is defined as where the sum is over the edges in the shortest path between nodes i and j, w k is the weight of edge k and w 0 is the minimum non-zero weight, introduced for normalization purposes [7]. The average distance to the seed, , where the sum is over nodes that do not belong to the seed and N = 246 is the total number of nodes in the network (i.e. the number of ROIs), can be used as a topological proxy for the propagation of the SI dynamics. The topological measure that quantifies the effect of a given virtual resection is thus the increase in the average distance to the seed:

Optimization method
In order to find optimal resection strategies, and considering the strong non-linearities present in the model, we have made use of a probabilistic optimization algorithm, the Simulated Annealing (SA) algorithm [8], to optimize the effect of virtual resections, for a given number of resected nodes. In order to do so, we employed a four-step method to find, for each resection size S up to the size of the resection area, S RA , the resection leading to the maximum decrease in seizure propagation (see figure S2): 1. Firstly, we considered the surrogate model, and used the SA algorithm to maximize ∆R T . The SA algorithm is probabilistic, and was therefore run 10 times, yielding a set of surrogate solutions R k T , with k = 1, ..., 10, each characterized by its topological effect ∆R k T , as shown in figure S2a. Figure S3 shows an exemplary case of the SA optimization search.
2. The topological metric ∆R T is only an approximation of the actual effect each resection R has on the propagation dynamics. Consequently, for each surrogate solution R k T , its actual effect was then determined by estimating I k R T (t 0 ) (figure S2b), i.e. seizure propagation on the network once the virtual resection has taken place.
3. The optimal topological resection of size S, R * T (S), was defined as the one leading to the maximum decrease in epidemic spreading (figure S2c).
4. Finally, we made use of the SA algorithm again: the optimal topological resection R * T (S) was used as the initial condition for the search. This time, the algorithm was used to minimize I R (t 0 ) directly (see figure  S2d). This yielded the final optimal resection R * (S), with propagation I * R (S). This method was iterated from S = 1 up to the size of the resection area, S = S RA . Overall, this method allows for a faster exploration of the parameter space (using the surrogate model), and a posterior refinement of the solution using the seizure propagation model.

SA Algorithm
We used a MATLAB implementation of the Simulated Annealing (SA) algorithm (J. Vandekerckhove, general simulated annealing algorithm, version 1.0.0.0, MATLAB Central File Exchange), adapted for our model (see figure  S3). Starting out with a given selection of nodes for the virtual resection, at each step of the SA algorithm this set is updated by changing n of the nodes. n follows a power-law distribution to allow for non-linear exploration of the node-space. This SA algorithm minimizes a cost function, which was defined as the inverse of d seed R for the surrogate case (step 1 in figure S2), and as the fraction of infected nodes at t 0 , I R (t 0 ) for the final step (step 4 in figure S2).
The other parameters of the algorithm were set to default for step 1: initial temperature = 1, stopping temperature = 10 −8 , annealing function f (T ) = 0.8T , maximum number of consecutive rejections = 10 3 , maximum number of tries within one temperature = 300, maximum number of successes within one temperature = 20. The algorithm was run 10 times with different (random) initial configurations, and the one with the largest effect was selected. For step 4, some parameters were changed to allow for a faster approach to the final solution: stopping temperature Figure S2: Schematic representation of the Virtual Resection Optimization Method. We used a four step optimization approach to find the virtual resection that, for a given resection size, minimized seizure propagation. Firstly (step 1), we used the SA algorithm to find resections that maximized the surrogate measure ∆R T (eq. S3). The optimization algorithm was iterated 10 times, yielding a set of 10 possible virtual resections, R k T (S), k = 1, ..., 10, for each resection size. This is depicted in panel (a), where we show d seed k as found by the SA algorithm for VRs of different sizes. In most cases all solutions had the same d seed k and overlap in the figure. For this case, the seed became completely disconnected for S = 5, which is 1 node less than the size of the RA. Next (step 2), we calculated the actual effect of each resection R k T (S), I k (t 0 ), as depicted in panel (b), where we show I k (t 0 ) for those same resections (error bars indicating the standard deviation among 10 iterations of 10 4 runs of the SI dynamics overlapped with data points and are not visible). Panel (c) illustrates the relationship between the surrogate and actual effects of each resection, i.e. between the topological measure d seed k and dynamical one I k (t 0 ), for each VR. I(t 0 ) decayed roughly proportionally with d seed until it saturated and equaled 0 for complete seed isolation (in which case d seed diverged, not shown here). In step 3, the resection leading to the maximum actual effect for each resection size, R * T (S), was selected as shown in panel (d). Finally, for each resection size separately, this resection was used as the initial condition for the SA algorithm to minimize seizure propagation, as measured by I R (t 0 ) (panel (e)). This led to the final optimal resection, R * (S), characterized by the propagation I * R (S), for each resection size. In this case, we found that the propagation was stopped for S = 5.
= 10 −6 , annealing function f (T ) = 0.6T . In this case the algorithm was only run once for each resection size, using the optimal solution from the surrogate SA method as the seed.
Due to the random nature of the SI dynamics, different realizations with the same seed can lead to different outcomes. The SA algorithm could thus artificially minimize I(t 0 ) by selecting realizations that, by chance, have a SA -Simulated Annealing Figure S3: Surrogate Simulated Annealing. Here we illustrate the Simulated Annealing (SA) method, for an exemplary VR of size S = 4, for patient 10. The SA algorithm performs a random search, starting from a random virtual resection, so that at each step a new resection is calculated by switching some of the selected nodes. The new resection is then accepted with a probability that depends on the difference of effect between the resections, so that resections that increase d seed have a higher probability of being accepted. This process is iterated until a stable solution is found. Due to its probabilistic nature, different realizations of the algorithm can lead to different outcomes, of which the optimal one is selected. In this case, iterations 1, 2, 3, 5 reached the same resection (target nodes 19, 30, 33, 17) that maximally increased d seed .  Table S1: Model optimization. For each patient the maximum correlation between the clinical and model seizure patterns, and the corresponding optimal network connectivity, for the broad-band (BB) networks, is given, both for the RA and SOZ seeds. The seizure outcome is also indicated. The final three rows indicate the average correlation and density obtained for the whole population (all), seizure free patients (SF) and not seizure free patients (NSF). * Indicates a significant correlation (p < 0.05 ).

RA
smaller propagation at t 0 . Averaging over a large number of realizations (10 4 in the current case) minimizes this problem, but it can still occur. To avoid this artificial minimization, after accepting a new virtual resection R the propagation I R (t 0 ) was calculated again with a new simulation of the SI dynamics. In this manner the realizations that were used to select a resection were not carried over in subsequent iterations. This step was not necessary in case of the surrogate model since the calculation of d seed R is exact.
3 Reproduction of seizure propagation patterns: Results.

Optimal correlation for BB networks.
In table S1 we report the optimal correlation found for each patient using the patient-specific BB-connectivity network, and the corresponding value of the mean connectivity.

Results for α−band networks.
We report here on the results of the seizure propagation model fit and validation analysis for the α-band networks. As described in the methods section of the main text, we measured the correlation C between the SI and SEEG propagation patterns for different values of the mean network connectivity κ. The resulting C(κ) curves are shown in figure S4, in panels (a) and (b) for the RA and SOZ seeds, respectively. As it was the case for the BB networks, we found a bimodal dependence of C(κ) for most patients, so that there is a maximum for low density values and another one for large densities. As before, we restricted our considerations to the first maximum, which leads to networks with only the fundamental connections that are needed to reproduce seizure propagation. Then, we found the connectivity value κ max that led to the maximum correlation, independently for each patient. We found that the model significantly reproduced the seizure propagation patterns for 9 and 10 out of 10 patients for the RA and SOZ seeds, respectively. The maximum correlation C max found for each patient is indicated in figure S4 with a triangle, and the corresponding values are reported in table S2 and figure S5.
We also report in figure S5 the optimal correlation found when using the average α-band network and the fully connected network. The average values for each group are reported on the main text. Finally, in figure S4c we show the average model as given by the median correlation of C(κ) for all patients, for both the RA and SOZ seeds. Figure S4: Reproduction of seizure propagation. Correlation between the modelled and clinically observed seizure propagation pattern as a function of network density, for each patient, using the RA (panel a) and SOZ seeds (panel b), for the alpha-band networks. Panel (c) indicates the median curves for each case, as indicated by the legend. Circles in panels (a) and (b) denote significant correlations (p < 0.05), and the optimal correlation for each case is marked with a triangle.  Table S2: Model optimization. For each patient, the maximum correlation between the clinical and modelled seizures patterns, and the corresponding optimal network connectivity, for the α−band networks, is given for the RA and SOZ seeds. The seizure outcome is also indicated. The final three rows indicate the average correlation and density obtained for the whole population (all), seizure free patients (SF) and not seizure free patients (NSF). * Indicates a significant correlation (p < 0.05 ).   Figure S6 reports the dependence of the normalized size of the 90% and 100% (panel a) and of the normalized effect of the 1-node resection (panel b) on the number of links of the resection area, E RA . E RA combines the effect of the size and connectivity of the resection area. We found that s 90 and s 100 correlated positively and significantly with E RA (r(8) = 0.81, p = 0.005 and r(8) = 0.61, p = 0.07, respectively for s 90 and s 100 ) although the correlation was only significant for s 90 . We analyzed whether this effect was led by the outlier case with largest E RA by removing it from the correlation analysis. The results (r(7) = 0.85, p = 0.004 and r(7) = 0.64, p = 0.06, respectively for s 90 and s 100 ) indicate that this was not the case case. i R1 correlated negatively and significantly with E RA (r(8) = −0.71, p = 0.02), and the effect was also not led by the outlier case (r(7) = −0.65, p = 0.05). Thus, larger or more densely connected epileptogenic regions would require larger minimal resections, as expected.