# Effects of the spatial resolution of the Virtual Epileptic Patient on the identification of epileptogenic networks

Jean-Didier Lemaréchal, Paul Triebkorn, Anirudh Nihalani Vattikonda, Meysam Hashemi, Marmaduke Woodman, Maxime Guye, Fabrice Bartolomei, Huifang E Wang, Viktor Jirsa

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Jean-Didier Lemaréchal, Paul Triebkorn, Anirudh Nihalani Vattikonda, Meysam Hashemi, Marmaduke Woodman, et al.. Effects of the spatial resolution of the Virtual Epileptic Patient on the identification of epileptogenic networks. Imaging Neuroscience, 2024, 2, pp.1 - 22. 10.1162/imag_a_00153. hal-04636707

HAL Id: hal-04636707

https://hal.science/hal-04636707v1

Submitted on 5 Jul 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Effects of the spatial resolution of the Virtual Epileptic Patient on the identification of epileptogenic networks

Jean-Didier Lemaréchal<sup>a</sup>, Paul Triebkorn<sup>a</sup>, Anirudh Nihalani Vattikonda<sup>a</sup>, Meysam Hashemi<sup>a</sup>, Marmaduke Woodman<sup>a</sup>, Maxime Guye<sup>b,c</sup>, Fabrice Bartolomei<sup>a,d</sup>, Huifang E. Wang<sup>a,*</sup>, Viktor Jirsa<sup>a,*</sup>

Corresponding Author: Jean-Didier Lemaréchal (jean-didier.lemarechal@sorbonne-universite.fr)

## ABSTRACT

Digital twins play an increasing role in clinical decision making. This study evaluates a digital brain twin approach in presurgical evaluation, the Virtual Epileptic Patient (VEP), which estimates the epileptogenic zone in patients with drug-resistant epilepsy. We built the personalized digital brain twins of 14 patients and a series of synthetic dataset by considering different spatial configurations of the epileptogenic and/or propagation zone networks (EZN and PZN, respectively). Brain source signals were simulated with a high spatial resolution neural field model (NFM) composed of 81942 nodes, embedding both long-range (between brain regions) and short-range (within brain regions) coupling. Brain signals were then projected to stereotactic electroencephalographic (SEEG) contacts with an accurate forward solution. An inversion procedure based on a low spatial resolution neural mass model (NMM) composed of 162 nodes was applied to estimate the excitability of each region in each simulation. The ensuing estimated EZN/PZN was compared to the simulated ground truth by means of classification metrics. Overall, we observed correct but degraded performance when using an NMM to estimate the EZN from data simulated with an NFM, which was significant for the simplest spatial configurations. We quantified the reduced performance and demonstrated that the oversimplification of the forward problem is its principal cause. We showed that the absence of local coupling in the NMM affects the inversion process by an overestimation of the excitability, representing a significant clinical impact when using this procedure in the context of presurgical planning. In conclusion, this study highlighted the importance to shift from an NMM towards a full NFM modeling approach for the estimation of EZN, with a particularly relevant need when considering the most complex clinical cases.

**Keywords:** Virtual Epileptic Patient (VEP), epileptogenic network, stereotactic electroencephalography, neural field model, neural mass model, inference

## 1. INTRODUCTION

### 1.1. Introduction of the VEP

Epilepsy is a common neurological disorder characterized by recurrent seizures (Fisher et al., 2014; Noebels et al.,

2012). Epileptic seizures very often involve networks of connected brain regions (Bartolomei et al., 2017; Besson et al., 2017) and present complex spatiotemporal brain dynamics (da Silva et al., 2003; Kalitzin et al., 2019; Saggio et al., 2020). For patients with pharmaco-resistant

<figure><img src="image_3.png" /></figure>

<footnote><sup>a</sup>Aix-Marseille Université, Institut National de la Santé et de la Recherche Médicale, Institut de Neurosciences des Systèmes (INS) UMR1106, Marseille, France</footnote>

<footnote><sup>b</sup>Aix-Marseille Université, CNRS, CRMBM, Marseille, France</footnote>

<footnote><sup>c</sup>APHM, Timone University Hospital, CEMEREM, Marseille, France</footnote>

<footnote><sup>d</sup>APHM, Epileptology and Clinical Neurophysiology Department, Timone Hospital, Marseille, France</footnote>

<footnote>*Huifang E. Wang and Viktor Jirsa contributed equally to the work</footnote>
focal epilepsy, resective surgery is one invasive solution to become seizure free. The goal of this strategy is to remove a sufficient and minimal set of epileptogenic regions which are responsible for seizure initiation (Jehi, 2018; Lüders et al., 2006; Rosenow & Lüders, 2001). Stereotactic electroencephalographic (SEEG) is one of the tools used during presurgical planning in order to identify those early recruited regions (Bartolomei et al., 2018; Bernabei et al., 2023; Cardinale et al., 2019; Mercier et al., 2022). This procedure consists in implanting a series of electrodes in the brain of one patient. The location of these electrodes is chosen to test and validate a clinical hypothesis derived from other diagnostic techniques (imaging, EEG, clinical tests). Brain activity is recorded continuously at the milliseconds time scale along the contacts of the depth electrodes separated from a few millimeters. These very high spatiotemporal resolved signals may record different kinds of typical epileptic activity (interictal epileptiform discharges, spontaneous or induced seizures), which are then used to infer the location of the targeted regions for the surgery (Lagarde et al., 2019). To better estimate the localization of the seizure onset zone based on SEEG recordings, many data-driven approaches based on the analysis of functional connectivity have been developed to identify and characterize the dynamics of interictal and ictal brain activity, as well as transitions between them (Wendling et al., 2010), such as the epileptogenicity index (Bartolomei et al., 2008). To complement and extend this approach, a series of computational models have been developed in order to capture specific dynamic features (Wendling et al., 2016). These models, either biophysical or more phenomenological, usually reduce these complex and non-linear dynamics to a set of coupled differential equations which are solved numerically to form solutions. Importantly, they represent an efficient tool to integrate multimodal data (structural, functional, pathophysiological) of one patient to obtain a personalized model (V. Jirsa et al., 2023). Moreover, when considering empirical data of one specific patient, the parameters of these models can be estimated based on inference techniques, such that the predictions of the model best fit the observations of the patient.

In this context, the Virtual Epileptic Patient (VEP) has been introduced as a personalized generative whole brain model of epilepsy spread (V. K. Jirsa et al., 2017; Wang et al., 2023). The VEP integrates the individual structural data of the patient (brain geometry, structural connectivity, location of SEEG electrodes) with a phenomenological model of epileptic activity (V. K. Jirsa et al., 2014). In particular, it implements a straightforward and efficient way of parameterizing the virtual epileptogenic network (Bartolomei et al., 2017) of the

patient in terms of an epileptogenic zone network (EZN), a propagation zone network (PZN), and non-involved zones (NIZ). The distinction between these three networks follows the actual concept of a hierarchical organization taking place during a focal seizure within the brain of a patient with epilepsy. The EZN gathers brain regions with high epileptogenicity which are able to trigger seizures. After the initiation of a seizure, abnormal ictal activity may propagate to regions with lower epileptogenicity, the PZN. Finally, brain regions which do not show any significant change of their activity during a seizure are considered as NIZ. One major application of the VEP is the estimation of the EZN based on empirical SEEG recordings of spontaneous or stimulation-induced seizures (model inversion). Once the EZN has been accurately identified, the VEP can also reliably predict the spatiotemporal pattern of propagation (Proix et al., 2017; Wang et al., 2023). Notably, this specific VEP approach falls within a broader context of computational modeling methods (Bernabei et al., 2023; Lytton, 2008) that have been developed to estimate EZN and predict effect of therapy (Goodfellow et al., 2016; Kini et al., 2019; Olmi et al., 2019; Sinha et al., 2017). Given that currently, surgery does not systematically cure the patients (Cardinale et al., 2019), these studies show how modeling can bring significant contribution in planning therapeutic strategies tailored to the specificity of each patient and that collective efforts are still needed for a better clinical translation of these computational techniques (Litt, 2022).

## 1.2. Distinction between NMM and NFM

In this study, we focus the attention on two essential attributes of a VEP: the spatial model and the neuronal model.

The spatial model is either a Neural Mass Model (NMM) or a Neural Field Model (NFM) (Deco et al., 2008). In brief, the key difference between NMM and NFM is their spatial resolution. On the one hand, NMM is a brain model specified as a finite number of regions corresponding to macroscopic parts of the brain (a few $cm^2$ or $cm^3$) and usually defined by a specific parcellation. In this study, we considered the VEP atlas (Wang et al., 2021) comprising 162 regions (73 cortical and 8 subcortical regions per hemisphere). These regions are represented as point sources interacting with each other through long-range connections (order of cm). On the other hand, NFM overcomes the spatial limitation of NMM by considering the brain as a spatially continuous network (V. K. Jirsa, 2009). This has at least two significant implications for NFM. First, in addition to a long-range (cm) coupling between distant macroscopic brain regions, that is, global connectivity, a short-range (mm) coupling, that is, local
connectivity, can be modeled. This property is crucial for modeling seizure spread along cortical or subcortical surfaces (Liou et al., 2020; Martinet et al., 2017; Schlafly et al., 2022). And second, the projection of the simulated brain activity to the SEEG signals recorded at the electrode level, that is, the forward solution, can be modeled with better accuracy. Indeed, it allows to model brain sources as electrical dipoles and to consider their relative position and orientation with respect to the SEEG contacts, which is essential when computing electromagnetic propagation of electrical activity. In practical implementations, NFM is typically discretized on brain surfaces rather than being discretized point-wise, as is done in the NMM approach. In this study, we use 81924 vertices for the cortical surface. Necessarily, such very high spatial resolution of NFM, compared to the 162 nodes of NMM, comes with a significant additional computational cost.

In the present context, the neuronal model used at each node of the network is either the full Epileptor with 5 state variables (5D) in its original implementation (V. K. Jirsa et al., 2014) or its two-dimensional (2D) reduction (Proix et al., 2014).

## 1.3. Development of the VEP

Since the introduction of the VEP (Proix et al., 2017; V. K. Jirsa et al., 2017), many studies have developed, improved, and validated the approach. From a methodological point of view, the VEP has been specified as a fully Bayesian probabilistic generative model (Hashemi et al., 2020), in which the importance of informative pri-ori knowledge has been highlighted (Hashemi et al., 2021). VEP has been validated with synthetic data and applied to empirical SEEG recordings for the estimation of the EZN with different inference techniques: maximum a posteriori (Vattikonda et al., 2021), Hamiltonian Monte Carlo (Jha et al., 2022), and simulation-based inference (Hashemi et al., 2023). From a clinical perspective, the performance of the VEP has been recently assessed retrospectively in a cohort of 53 patients by comparing the predictions of the VEP with the predictions of clinical experts, with respect to the surgical outcome of the patients (Makhalova et al., 2022; Wang et al., 2023). One essential conclusion was that the concordance of the predictions was better for seizure-free compared to non-seizure-free patients. In addition, it is worth mentioning that the VEP methodology we are considering in the present study to identify EZN is currently at the heart of an ongoing large prospective multi-centric clinical trial (EPINOV, NCT03643016; V. Jirsa et al., 2023), with more than 300 patients, in order to evaluate its abilities to help plan surgery strategies.

## 1.4. Objective of the study

All the previously mentioned studies have used NMM to estimate the EZN from empirical observations. Indeed, considering the high spatial resolution of NFM, using NFM for model inversion is still a very difficult task (Vattikonda et al., 2023). For validation, that is, comparing the parameters estimated from synthetic data to the parameters used to simulate the same synthetic data (ground truth), all these studies have simulated synthetic observations with NMM. In other words, these studies have used only NMM for both simulation and estimation.

As emphasized previously, the repertoire of spatio-temporal dynamics generated with NFM is significantly more extended than that of NMM. Indeed, some specific patterns observed in empirical data, such as traveling waves, cannot be explained with NMM.

To push further the validation process of NMM for the estimation of the EZN, and demonstrate the relevance of its usage in a clinical context, it is therefore necessary to simulate synthetic observations with high spatial resolution NFM. This procedure has already been proposed in Wang et al. (2023) and the authors demonstrated its proof of concept for one dataset. In the present study, we carried out an extensive and systematic evaluation of this approach in order to address two fundamental questions: what is the performance of the neural mass-based VEP pipeline when the ground truth simulated with NFM is known? And under which conditions would it succeed or fail?

## 1.5. Scheme of the study

In this study, we evaluated the performances of VEP with 14 patients with epilepsy. For each patient, we created a series of spatial configurations (ground truth), each defined with different EZN and PZN. We also varied the regions involved in these networks in terms of their size and their number. To disentangle the effects of using a simplified forward solution from the effects of absence of local coupling, two modeling features that an NMM inversion is intrinsically lacking, we simulated data in different conditions: with an NMM (control condition), with an NFM with local coupling equal to 0 (to evaluate the effects of forward model only), and with an NFM with local coupling equal to 2 and 4 (to evaluate the effects of local coupling only). Then, for each simulation, the EZN was systematically estimated with an NMM equipped with the reduced 2D Epileptor. Predictions were evaluated at the individual and group level in terms of goodness of fit and classification metrics. Finally, we separately demonstrated the effects of the two key limitations of the NMM-based estimation, namely the simplified forward solution and the absence of local propagation, with two specific examples.
# 2. MATERIALS AND METHODS

## 2.1. Selection of patients and data acquisition

We selected fourteen patients suffering from drug-resistant focal epilepsy who underwent a standard presurgical protocol at La Timone Hospital in Marseille. Informed written consent was obtained for all patients, and the study was approved by the local ethics committee (Comité de Protection des Personnes sud Méditerranée 1).

The evaluation included non-invasive T1-weighted (T1W) MRI (magnetization prepared rapid acquisition gradient echo sequence, either with repetition time = 1.9 s, echo time = 2.21 ms (1 patient, id010 in Supplementary Table 1) or repetition time = 2.3 s, echo time = 2.98 ms (13 patients), voxel size 1.0 × 1.0 × 1.0 mm) and diffusion-weighted (DW) MRI (diffusion tensor imaging-MR sequence, either with angular gradient set of 64 directions, repetition time = 10.7 s, echo time = 95 ms, voxel size 1.95 × 1.95 × 2.0 mm, b-weighting of 1000 s/mm² (2 patients, id010 and id022 in Supplementary Table 1), or with angular gradient set of 200 directions, repetition time = 3 s, echo time = 88 ms, voxel size 2.0 × 2.0 × 2.0 mm, b-weighting of 1800 s/mm² (12 patients)). Images were acquired using a Siemens Magnetom Verio 3T scanner.

SEEG electrodes (10–18 contacts 2 mm long and separated by 1.5 mm for all patients) were implanted individually for each patient as part of the clinical routine and according to the hypotheses about the EZN. A post-implantation computed tomography (CT) scan was performed to obtain the location of the implanted electrodes.

These patients were selected so as a majority of regions from each of the four cortical lobes of the brain (frontal, parietal, occipital, and temporal) could be represented. We considered that it was unnecessary to include more patients, especially if the clinical hypothesis about the localization of their EZN (and the related implantation) was similar to the ones already present in the study.

## 2.2. General workflow

This section summarizes the pipeline used in the present study (Fig. 1). Specific details regarding the modeling, the creation of the synthetic dataset, and the estimation of the EZN with NMM and its evaluation are provided in subsequent sections.

First, a common procedure was applied to build an individualized VEP for each patient (Fig. 1A). A brain segmentation was computed from the T1W MRI with the FreeSurfer software package (http://surfer.nmr.mgh.harvard.edu) and parcellated with respect to the VEP atlas (Wang et al., 2021). We used the ico6 spatial resolution, providing a cortical surface composed of 81924 vertices (~1.5 mm

between the vertices). The structural connectome was derived from the DW MRI with the MRtrix software package (https://www.mrtrix.org) and the structural connectivity matrix (SC), corresponding to long-range (or global) coupling, was computed by counting the streamlines connecting the parcels of the VEP atlas (Vattikonda et al., 2021; Wang et al., 2023). The location of each contact of the SEEG electrodes was obtained from the CT scan using GARDEL (Medina Villalon et al., 2018) and coregistered to the MRI space. Together with the brain surfaces and the SC, these elements constitute the anatomical structure of an individualized VEP and constrain the spatiotemporal dynamics of the neuronal activity generated within the brain of the patient and measured with the SEEG electrodes. Second, we built a series of {EZN, PZN} spatial configurations based on individual clinical hypotheses. Clinical hypotheses (Supplementary Table 1) contain the list of regions suspected to be part of EZN or PZN. They were provided by the clinicians of La Timone Hospital in Marseille based on the presurgical evaluation of each patient, including the computation of an epileptogenicity index for each brain region (Bartolomei et al., 2008).

Then, we simulated each configuration in 4 different conditions (Fig. 1B): with an NMM and with an NFM for increasing strengths of the local coupling (Ic: 0, 2, and 4. The 4 conditions were simulated using the full 5D Epilepitor and labeled respectively Ic_none, Ic_0, Ic_2, and Ic_4. Increasing the strength of local coupling increases the speed of propagation of brain activity (Proix et al., 2018). It makes regions seize earlier (Fig. 3B, C and D, first column: the onset time of contact GC6-7 in simulated data), but may also change the number of recruited regions. This range of local coupling strengths also needs to maintain a balance with the other parameters, such as structural connectivity or global coupling (Spiegler & Jirsa, 2013). Indeed, for higher values of local coupling strength, the propagation of activity gets dominated by local coupling (fast propagation) and regions in PZN seize immediately after regions in EZN: it becomes almost impossible to distinguish EZN from PZN. Finally, simulated brain signals were multiplied by a gain matrix to form the measurements at the SEEG contacts, that is, the forward solution. Then, for each synthetic SEEG dataset (the observations), we used a probabilistic inversion scheme based on an NMM equipped with the reduced 2D Epilepitor (Fig. 1C). This automatic procedure consists in finding the most plausible values (maximum a posteriori, Section 2.5.2) for each of the free parameters of the NMM of the patient such that the predictions of the model best fit the observations.

The quality of the estimation was evaluated by two ways (Fig. 1D). We computed a goodness of fit to evaluate how
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** General workflow of the study. (A) Structural neuroimaging data of the patient are used to build a personalized brain network model: the Virtual Epileptic Patient (VEP). In particular, the T1-MRI, the DW-MRI, and the CT scan specify the brain space, the long-range connectome, and the position of the SEEG contacts respectively, which are essential attributes of a personalized VEP. (B) Using a specific spatial configuration of the Epileptogenic and Propagation Zone Networks (EZN and PZN respectively) and a high spatial resolution Neural Field Model (NFM) composed of 81924 cortical and subcortical nodes, synthetic data are simulated at the brain level for all the nodes (top panel: EZN and PZN nodes are shown on the first and second row respectively) and at the SEEG level (bottom panel). (C) Data features are extracted from the SEEG time series, and a low spatial resolution Neural Mass Model (NMM) composed of 162 nodes is used to predict the epileptogenicity (or excitability) of each node. (D) The estimation of the EZN and PZN is evaluated by comparing the ground truth excitabilities used to simulate the data, with the predicted ones (NIZ: Non-Involved Zones).</figcaption></figure>

well the SEEG signals predicted by the statistical model fitted the observations. Based on the estimated parameters of the model, we also derived an estimation of the {EZN, PZN} that we compared to the ground truth {EZN, PZN} used to simulate the data by means of classification metrics (precision, recall, F1 score).

## 2.3. Modeling brain activity

### 2.3.1. Neural field model

In this section, we describe the NFM brain network model used to simulate synthetic data for the conditions lc_0, lc_2 and lc_4.
**2.3.1.1. Simulating brain activity.** We used a mixed approach to represent the brain space, namely an NFM for the cortical structures and an NMM for the subcortical structures.

One reason to justify the importance of using NFM when modeling SEEG signals is the possibility to compute an accurate forward solution when the orientation of each source (with respect to the SEEG electrodes) is well defined. This is indeed the case in cortical structures where neuronal populations are oriented perpendicular to the surface and each node (vertex) can express its own orientation. However, there is so far no consensus in the literature about such a systematic orientation property in the subcortical structures, and therefore, no proper way

to implement an accurate forward solution. Because of this limitation, subcortical structures were modeled with NMM and their associated forward solution (Section 2.3.2). Nevertheless, the fact that in this study, we restricted the location of EZN and PZN to cortical structures mitigates the effects of this simplification.

In this context, the brain space was defined with a high spatial resolution for the cortical surface (81924 nodes) and a low spatial resolution for the subcortical structures (18 nodes).

The full five-dimensional Epileptor neuronal model was assigned to each node and activity of node$_i$ located at x$_i$, with $1 \le i \le 81924 + 18$, was given by the following equations (Proix et al., 2018):

$$ \dot{u}_{1,i} = u_{2,1} - f_1(u_{1,i}, q_{1,i}) - v_i + I_1 + LC_{scaling}\gamma_{11}w_1 * S(u_{1,i}, \theta_{11}) + GC_{scaling}\sum_{j=1}^{162} C_{i,j}S(\tilde{u}_{1,j}, \theta_{GC}) $$

$$ \dot{u}_{2,i} = 1 - 5u_{1,i}^2 - u_{2,i} $$

$$ \dot{v}_i = \frac{1}{\tau_0}(4(u_{1,i} - u_0(x_i)) - v_i) $$

$$ \dot{q}_{1,i} = -q_{2,i} + q_{1,i} - q_{1,i}^3 + I_2 + 0.002g(u_{1,i}) - 0.3(v_i - 3.5) + LC_{scaling}\gamma_{22}w_2 * S(q_{1,i}, \theta_{22}) $$

$$ \dot{q}_{2,i} = \frac{1}{\tau_2}(-q_{2,i} + f_2(q_{1,i})) $$

with

$$ g(u_{1,i}) = \int_{t_0}^{t} e^{-(t-s)/\tau_{12}} (a_{12}u_{1,i} + LC_{scaling}\gamma_{12}w_{12} * S(u_{1,i}, \theta_{12}))ds $$

$$ f_1(u_{1,i}, q_{1,i}) = \begin{cases} u_{1,i}^3 - 3u_{1,i}^2 & \text{if } u_{1,i} < 0 \\ (q_{1,i} - 0.6(v_i - 4)^2)u_{1,i} & \text{if } u_{1,i} \ge 0 \end{cases} $$

$$ f_2(q_{1,i}) = \begin{cases} 0 & \text{if } q_{1,i} < -0.25 \\ 6(q_{1,i} + 0.25) & \text{if } q_{1,i} \ge -0.25 \end{cases} $$

with $I_1 = 3.1$, $I_2 = 0.8$, $\gamma_{11} = 0.17$, $\gamma_{22} = 0.32$, $\gamma_{12} = 0.03$, $\theta_{11} = -1$, $\theta_{GC} = -1$, $\theta_{22} = -0.5$, $\theta_{12} = -1$, $\tau_0 = 2857/2$, $\tau_2 = 10$, $\tau_{12} = 100$, $a_{12} = 3$. Here, all parameters and variables have arbitrary units. Briefly, this phenomenological model couples together three populations, $u_i$, $v_i$ and $q_i$, evolving at different time scales. With a slow time scale, the population $v_i$, usually referred as the slow permittivity variable, controls when the system enters and leaves the ictal state, which crucially depends on the fixed excitability $u_0(x_i)$ of node$_i$, that is, its aptitude to autonomously generate seizure without any input from other nodes of the network. The populations with fast and intermediate time scales, $u_i$ and $q_i$, account respectively for the presence of (low-voltage) fast oscillations and spike-wave discharges

during the ictal state. Default values of the parameters have been chosen to best mimic characteristic features of spontaneous seizures recorded in patients with epilepsy (and particularly the presence of fast oscillations and spike-wave discharges; V. K. Jirsa et al., 2014), as well as the transitions between interictal and ictal activity (direct current shift at seizure onset, slowing down of activity at seizure offset). The observable activity of one Epileptor node$_i$ is the quantity $q_{1,i} - u_{1,i}$ which represents the cooccurrence of spike-wave discharges with fast oscillations.

$LC_{scaling}$ and $GC_{scaling}$ are the coefficients used to scale local and global coupling at the whole network level, respectively. The local coupling terms $w * S(u, \theta)$ correspond to the spatial convolution at position $x_i$ of a
Laplacian kernel w (cutoff: 10 mm) with the local firing rate implemented with the Heaviside step function: S(u,θ) = 0 if u < θ and 1 otherwise. C<sub>i,j</sub> is the connection strength between the region node<sub>i</sub> is belonging to and region j and ũ<sub>1,j</sub> is the averaged activity of region j across all nodes belonging to the region. For more details regarding the general mechanic of the neuronal model, the technical implementation of the coupled oscillators u<sub>i</sub> and q<sub>i</sub>, the value taken by each parameter, or the use of the model in the context of NFM, we refer the interested reader to V. K. Jirsa et al. (2014) and Proix et al. (2018).

For subcortical nodes modeled as point sources NMM, the equations read the same except that the local connectivity is not considered (γ<sub>11</sub> = γ<sub>22</sub> = γ<sub>12</sub> = 0).

Whole-brain dynamics was integrated with a fourth-order Runge-Kutta deterministic integration scheme (step size: 0.1) for 20000 samples. In this study, we chose the millisecond as the unit of time and each simulation was therefore 2000 ms long. Initial conditions were set so that the simulated network could start in the interictal state and far away from seizure onset.

2.3.1.1.1. Dynamics of the Epileptor model. Linear stability analysis of the Epileptor (a complete bifurcation analysis is proposed in El Houssaini et al., 2020) has shown that depending on its excitability u₀, the system could exhibit different stability regimes, as described in previous studies (Proix et al., 2014; V. K. Jirsa et al., 2017) and used in previous simulation and inference investigations (Hashemi et al., 2020, 2021; Jha et al., 2022; Vattikonda et al., 2021). For an excitability u₀ lower than a critical value -2.05, and without any external input from other nodes of the network (no coupling), an isolated Epileptor has a single stable fixed point and does not seize autonomously. When the excitability equals the critical value, a saddle-node bifurcation occurs and an unstable fixed point appears. Then, seizures happen spontaneously and the node is part of the Epileptogenic Zone Network (EZN). When excitability is lower but still sufficiently close to the critical value, a relatively small external input from the network may destabilize the system and trigger a seizure in the recruited node: the node is part of the Propagation Zone Network (PZN). Otherwise, if the excitability is lower and far enough from the critical value, the external input cannot trigger a seizure and the node is in the Non-Involved Zone (NIZ).

Importantly, this means that one single threshold to distinguish PZN nodes from NIZ nodes across the whole network does not exist in general. Instead, this threshold depends on the intrinsic properties of the node under consideration, the configuration of the network activity, and the coupling parameters (SC, GC<sub>scaling</sub> and LC<sub>scaling</sub>).

2.3.1.2. Mapping brain activity to SEEG measurements. For cortical nodes, the gain from node<sub>i</sub> to SEEG contact k located respectively at positions x<sub>i</sub> (vertex i) and r<sub>k</sub> was computed by considering the brain as an infinite homogeneous volume conductor with the following equation (Sarvas, 1987):

$$g_{k,i} = \frac{a_i}{4\pi\sigma} \mathbf{Q}^* \frac{\mathbf{r}_k - \mathbf{x}_i}{|\mathbf{r}_k - \mathbf{x}_i|^3}$$

where a<sub>i</sub> is the area of vertex i (one-third of the sum of surrounding triangles areas) and counterbalances the (relatively small) spatial inhomogeneities due to the surface discretization, σ is the homogeneous electrical conductivity, Q is the unitary dipolar moment (the unit vector normal to the surface at vertex i), and Q*(r<sub>k</sub> - x<sub>i</sub>) is the dot product of Q and the vector (r<sub>k</sub> - x<sub>i</sub>), going from vertex i to SEEG contact k. For subcortical nodes, the gain was computed using the NMM equation (Section 2.3.2).

SEEG measurements were computed with the following equation:

$$SEEG_k(t) = \sum_{i=1}^{81924+18} g_{k,i} (q_{1,i}(t) - u_{1,i}(t))$$

where q<sub>1,i</sub> - u<sub>1,i</sub> is the instantaneous observable activity of Epileptor node<sub>i</sub> (Section 2.3.1.1).

## 2.3.2. Neural mass model

In this section, we describe the NMM brain network model used to simulate synthetic data for the condition lc_none.

2.3.2.1. Simulating brain activity. The point source network model was composed of 162 nodes (cortical: 146, subcortical: 18). The equations read the same as the equations of NFM (Section 2.3.1.1) except that 1 ≤ i ≤ 162, no local connectivity is considered (LC<sub>scaling</sub> = 0) and the averaged activity ũ<sub>1,j</sub> of region j is the activity of its single node u<sub>1,j</sub>: ũ<sub>1,j</sub> = u<sub>1,j</sub>.

2.3.2.2. Mapping brain activity to SEEG measurements. The gain from node<sub>i</sub> to SEEG contact k located respectively at positions x<sub>i</sub> (vertex i) and r<sub>k</sub> was computed with the following equation:

$$g_{k,i} = \frac{1}{4\pi\sigma} \sum_{j=1}^{N_i} \frac{a_j}{|\mathbf{r}_k - \mathbf{x}_j|^2}$$

where σ is the homogeneous electrical conductivity, and j identifies one of the N<sub>i</sub> vertices belonging to region i and a<sub>j</sub> is its area (Section 2.3.1.2).
Unlike NFM nodes, NMM regions cannot be properly represented by electrical dipoles with well-defined fixed position and orientation, due to their macroscopic sizes (Section 1.2). Instead, NMM systematically assumes orientations which maximize the gain values. Indeed, in the previous NFM gain equation (Section 2.3.1.2), if we assume that the dipolar moment **Q** of node *i* is systematically pointing into the direction of SEEG contact *k* (thus completely ignoring the curvature of the cortex and the normal to the surface at vertex *i*), then gain $g_{k,i}$ is maximal, dot product $\mathbf{Q} \cdot (\mathbf{r}_k - \mathbf{x}_i)$ reads $|\mathbf{r}_k - \mathbf{x}_i|$, and we obtain the NMM gain equation.

SEEG measurements were computed with the following equation:

$$SEEG_k(t) = \sum_{i=1}^{162} g_{k,i} (q_{1,i}(t) - u_{1,i}(t))$$

## 2.4. Configurations

For each of the 14 patients, we defined a series of {EZN/ PZN} spatial configurations following the same procedure. Each configuration defined the set of regions belonging to the EZN and to the PZN, based on the clinical hypothesis related to their specific EZN (Supplementary Table 1). Here, regions refer to the parcels of the VEP anatomical atlas (Wang et al., 2021), as stated in Section 1.2. We introduced 3 types of spatial configurations with increasing level of complexity. For the first type, only one region was included in the EZN. This resulted in a total of 62 configurations (corresponding to the 62 anatomical regions listed in Supplementary Table 1). In addition, to test the sensitivity of the estimation of the EZN with respect to the size of the region, we also considered a second series of configurations for which the size of each region was reduced to half. For this purpose, we computed the barycenter of each region. Starting from this cortical vertex, we used a region growing approach, successively adding surrounding vertices until the surface of the growing region reached half size of the surface of the full region. This procedure was followed by visual inspection. For the second type, one configuration was created for each patient in which all the regions from the hypothesis were included in the EZN (11 configurations, with at least 2 regions in the EZN, average number of regions in the EZN is 5.3 and maximum is 8). For the third type, one configuration was created for each patient for which the EZN contained one region and the PZN contained at least one region. Regions from the PZN were adjacent to the region of the EZN (11 configurations: 8 configurations with one region in PZN and 3 configurations with 2 regions in PZN). For each configuration, an excitability $u_0(x_i)$ was set at each node $x_i$ of the network to reflect the epileptogenicity

of the node, that is, its ability to trigger a seizure. In this study, considering the dynamics of the Epileptor model (Section 2.3.1.1), we set the following excitability values: -1.6 for the EZN (the node can seize in an autonomous way), -2.1 for the PZN (the node needs input from the network to be recruited and to seize), and -3.6 for the NIZ (the node does not seize). These values have been determined from the analysis of the Epileptor model dynamics (Section 2.3.1.1).

Each spatial configuration was then simulated in 4 conditions: with an NMM and with an NFM using 3 different scalings of the local coupling $LC_{scaling}$ (0, 2, and 4). These 4 conditions were simulated with the full five-dimensional Epileptor and in the following, they are identified as lc_none, lc_0, lc_2, and lc_4. The first condition lc_none corresponds to the best-case scenario, that is, a control condition in which an NMM is used for both the simulation and the estimation of the EZN. To dissociate the effects of the forward solution from the effects of the local coupling, we introduced the 2nd condition lc_0 which differs from NMM only by the way the forward solution is computed. Indeed, without local coupling, all the nodes belonging to the same region share the same dynamics (because of common initial conditions and a deterministic integration) which is also the same as the corresponding node from the 1st condition. And 3rd and 4th conditions added the effects of the local coupling.

For all configurations with empty PZN, the global coupling $GC_{scaling}$ was set to 1 whereas for configurations with non-empty PZN, the global coupling was adjusted such that the propagation was sufficient to trigger a seizure in the PZN. Importantly, for each simulation, we checked that parameters, and especially the excitabilities, were consistent with the expected behavior (EZN, PZN, or NIZ) of each node (Section 2.3.1.1). Specifically, we verified that without any coupling, only EZN nodes seized and with coupling, only EZN and PZN nodes seized. Table 1 summarizes the ensemble of simulated dataset, and Figure 2 illustrates one simulation with an NFM using EZN = {Right-T2-Anterior}, PZN = {Right-T2-posterior}, $GC_{scaling} = 6$ and $LC_{scaling} = 2$.

## 2.5. Estimation of the EZN

### 2.5.1. Extraction of data features

The estimation of the EZN was not performed directly from the SEEG raw time series (Fig. 2C) but from their envelopes (Fig. 2E), which capture the amount of fast oscillations present in the signal simulated with a full five-dimensional Epileptor (Proix et al., 2017; V. K. Jirsa et al., 2017). This data feature extraction resulted in a set of observations matching the kind of predictions expected
<table><caption>**Table 1.** Description of the synthetic dataset.</caption><thead><tr><th>Spatial configurations</th><th>Count</th><th>Neuronal model and local coupling (Ic)</th><th>Global coupling</th></tr></thead><tbody><tr><td>EZN = {1 region}, PZN = {}</td><td>62</td><td>• NMM (Ic)</td><td>1.0</td></tr><tr><td>EZN = {1 half size region}, PZN = {}</td><td>62</td><td>• NFM with Ic = 0 (Ic_0)</td><td></td></tr><tr><td>EZN = {multiple regions}, PZN = {}</td><td>11</td><td>• NFM with Ic = 2 (Ic_2)</td><td></td></tr><tr><td>EZN = {1 region}, PZN = {1 or 2 regions}</td><td>11</td><td>• NFM with Ic = 4 (Ic_4)</td><td>adjusted</td></tr></tbody></table>

We defined 4 types of spatial configurations (1st column) to generate a set of spatial configurations across patients (2nd column). Each configuration was simulated in 4 conditions (3rd column) associated with different values for the local (3rd column) and global (4th column) coupling.

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Synthetic data for patient id003. In this spatial configuration, EZN = {Right-T2-Anterior}, PZN = {Right-T2-posterior}. Local coupling strength is equal to 2. (A) Time courses of all vertices belonging to EZN (above) and PZN (below). (B) Spatial mapping of the onset latencies of the seizure. (C) Time courses of bipolar SEEG contacts. (D) SEEG implantation of the patient. (E) Data features (envelopes of high-frequency activity) computed from the SEEG bipolar contacts. The 4 bipolar contacts (1: A14-15, 2: TB14-15, 3: B15-16, 4: GC6-7) for which data features have the highest amplitude are highlighted in (C) and (D).</figcaption></figure>

from an NMM equipped with a reduced 2D Epileptor. Indeed, with 2 state variables acting on different time scales, such a neuronal model cannot generate fast oscillations. The data features were computed separately for each SEEG bipolar contact similarly to previous studies (Hashemi et al., 2021; Vattikonda et al., 2021). The bipolar signal was high pass filtered at 50 Hz and squared. The envelope was computed with a moving average (window length = 50 ms), log-transformed, low pass filtered

at 5 Hz, and normalized between 0 and 1. For the analysis, we only considered the time period during which ictal activity was present. Even if all simulated datasets were of same duration (2000 ms), the ictal period varied across dataset and we implemented a simple automatic procedure to detect the onset and the offset of each simulated seizure separately. It detected the 2 samples corresponding to the first up-crossing and last down-crossing ($t_{up}$ and $t_{down}$) of a given threshold (here, 0.1).
From the duration of this ictal interval, $d_{ictal} = t_{down} - t_{up}$, we only considered the extracted data features within the time interval $[t_{up} - d_{ictal}/2, t_{down} + d_{ictal}/2]$. The extracted signal was visually checked and downsampled to 512 time samples in order to avoid an potential bias during the fitting due to different length of the signal across conditions.

### 2.5.2. Estimation

The probabilistic approach of the estimation procedure has already been extensively detailed in previous studies (Hashemi et al., 2020, 2021, 2023; Jha et al., 2022; Vattikonda et al., 2021; Wang et al., 2023). Given a set of observations $x$ (the data features described in previous Section 2.5.1) and using Bayes' theorem, the posterior density of parameters $\theta$ is given by:

$$p(\theta|x) = p(\theta)p(x|\theta)/p(x)$$

where $p(\theta)$ is the prior distribution of $\theta$, $p(x|\theta)$ is the likelihood function, and $p(x)$ is the marginal likelihood (independent of $\theta$). In this study, the parameters were inferred with a maximum a posteriori (MAP) approach which estimates the mode of the posterior density:

$$\theta_{MAP}(x) = \arg \max_{\theta} p(\theta|x)$$

The computation of $\theta_{MAP}(x)$ was performed in Stan's probabilistic programming language (https://mc-stan.org) based on the iterative quasi-Newton optimization algorithm L-BFGS (Nocedal & Wright, 2006). At the first step, the parameters are initialized from their prior distribution using a random seed. Then at each step, the gradient of the posterior density is evaluated for the current parameters and used to update the parameters for the next step. Convergence is controlled by a set of tolerance values (sufficiently small norm of the gradient or sufficiently small change in posterior density or in parameter values). Because the estimation of MAP depends on the initial seed, L-BFGS algorithm was run 50 times, each time with a different seed. This procedure ended up with a set of $\theta_{MAP}(x)$ which was then used to build the posterior marginal distributions of each of the parameters (described in the next Section 2.6).

In order to evaluate the likelihood function (and its gradient), the inference process requires a generative model of brain activity, here a low spatial resolution NMM brain network model equipped with the 2D reduction of the Epileptor (Proix et al., 2014, 2017). With the same notations as before, the equations read:

$$\dot{u}_{1,i} = 1 - u_{1,i}^3 - 2u_{1,i}^2 - v_i + I_1 + GC_{scaling} \sum_{j=1}^{162} C_{i,j} S(u_{1,j}, \theta_{GC})$$

$$\dot{v}_i = \frac{1}{\tau_0} (4(u_{1,i} - u_0(x_i)) - v_i)$$

Note that to be consistent with the implementation of NFM and NMM equipped with full 5D Epileptor described in previous Section 2.3, we used a global fast synaptic coupling acting on population $u_{1,i}$ rather than a slow per-mittivity coupling acting on population $v_i$.

In this study, weakly informative priors were assigned to all free parameters of the model. In particular, $u_0(x_i) \sim N(-3.6, 1)$ for each node and $GC_{scaling} \sim N(1, 1)$. To reduce the dimensionality of the parameter space and limit the complexity of the inference, initial conditions were not inferred and fixed to the values used for simulation.

## 2.6. Evaluation of performances

The estimation of the EZN was evaluated for each SEEG simulation with one measure of goodness of fit (GOF) and classification metrics. All these measures range from 0 to 1 (1 corresponds to a perfect fit).

The GOF measured the similarity between the observed and the predicted SEEG signals and was evaluated for each seed as the ratio of the explained variance (predictions) to the observed variance (equations in Supplementary Material). One GOF was computed for each of the 50 seeds, and only the optimizations for which the GOF was superior to the 3rd quartile of the 50 GOF were selected. From the selected optimizations, the GOF of the simulation was defined as the median of the GOF and the posterior distributions of each parameter were obtained from the normalized histogram of posterior modes (the area under a normalized histogram integrates to 1).

Classification aimed at attributing one class (either EZ or PZ or NIZ) to each node, given an estimation of its excitability (here, the median of the posterior marginal distribution). However, in the present context, there is no direct way to attribute a class to each node because no fixed boundary threshold based on excitability exists to differentiate PZ from NIZ nodes (Section 2.3.1.1). To circumvent this problem, we reduced the classification from 3 classes to 2 classes, either EZ_PZ or NIZ, and proposed a simple heuristic threshold to differentiate EZ_PZ from NIZ. Given such a heuristic classification threshold (explicitly defined thereafter), the general procedure to derive classification metrics followed three steps: 1) a class was assigned to each node: EZ_PZ if its estimated excitability was larger than the classification threshold, and NIZ otherwise; 2) considering NIZ as the null hypothesis, we compared for each node the estimated class with the true
class and assigned either a type of error (false positive or false negative) or the kind of correct result (true positive or true negative); 3) we computed three classification metrics: precision (the fraction of nodes truly identified EZ_PZ among all nodes identified EZ_PZ), recall (the fraction of nodes truly identified EZ_PZ among all EZ_PZ nodes), and F1 score (harmonic mean of precision and recall). This general procedure was repeated for a series of classification thresholds spanning an extended range of values (from -5 to -1 with a step of 0.01). Finally, the classification threshold for which precision was maximal was selected, and the associated precision, recall, and F1 score were reported. This technique was applied independently for each simulation.

From a clinical perspective, the maximization of precision, equivalent to the minimization of false discovery rate, corresponds to a conservative choice if we consider estimated EZ_PZ regions as candidates for surgical resection.

Other classification approaches, particularly relevant in clinical context when ground truth is not known, have been proposed based on estimated excitabilities (Hashemi et al., 2020, 2021, 2023; Jha et al., 2022) or onset times extracted from predicted data features (Vattikonda et al., 2021, 2023; Wang et al., 2023). Yet, these approaches eventually suffer from the choice of arbitrary thresholds to classify PZN. Here, because we are in a simulation study and the ground truth is known, we sidestep this problem by computing a threshold maximizing the precision. This classification bias in favor of precision but at the expense of recall is counterbalanced in the reported F1 score.

Also, for clarity and simplicity, we did not systematically distinguish EZ from PZ (NIZ is the null hypothesis) because we assumed that both EZ and PZ estimations would be globally affected the same way by the factors under consideration (forward model and local coupling). Of course, the situation would be different in a clinical context, where the distinction between EZN and PZN is essential, for instance to plan optimal surgical strategies based on seizure propagation control (Olmi et al., 2019).

# 3. RESULTS

## 3.1. Examples

In this section, we illustrate key aspects of our results with two different spatial configurations of the EZN. These two didactic examples are illustrative of common phenomena observed repeatedly across simulations. For each configuration, and across the 4 simulated conditions (Ic_none: NMM without local coupling, Ic_0: NFM with no local coupling, Ic_2 and Ic_4: NFM with local coupling scaling equal to 2 and 4 respectively), we present the data features of the SEEG simulated observations,

the predictions of the NMM (equipped with the reduced 2D Epileptor) with highest GOF, the GOF for the 50 optimizations (corresponding to 50 different random seeds of the optimization procedure, Section 2.5.2), and the posterior distributions of the excitability of the regions along with the ground truth of the configuration.

The first example (Fig. 3) corresponds to patient id003 and EZN = {Right-T2-Anterior}, PZN = {Right-T2-posterior}. The synthetic raw data for the condition Ic_2 are illustrated on Figure 2. Observed and predicted SEEG time courses are presented respectively on the first two columns (Fig. 3). In NMM simulated data (Fig. 3A, first column), we clearly see around time sample 100 the initiation of the seizure triggered in EZN (Fig. 3A, last column, node 45) and best captured by contact B15-16, followed around time sample 200 by the propagation of the seizure in PZN (Fig. 3A, last column, node 46), this time very well captured by contact FCA13-14. Considering now classification metrics (Fig. 3, last column), they can be measured visually by comparing across nodes the true excitability (represented as a dot) and the median of the posterior distribution (represented by the intermediate horizontal bar of each violin plot): the closer they are, the better the classification metrics. Once the classification threshold (represented as a pink horizontal dotted line) has been estimated (Section 2.6), all nodes for which the median of posterior distribution is higher than the threshold are positive: either true positives if they belong to EZN or PZN or false positives for those in NIZ. Precision is then the ratio of true positives among all positives, and recall is the ratio of true positives among all EZN or PZN nodes (sum of true positives and false negatives). In this example, the estimation of the EZN is very good in all 4 conditions with precision = 1 and recall = 1 except for the condition Ic_0 (Fig. 3B, last column) for which recall = 0.5 (the posterior median excitability is higher in region 51 (NIZ) compared to region 46 (PZN)). However, considering the NFM simulations, two points are noteworthy. First, the underestimation of the excitability in PZN (region 46) is related to the data features themselves, which show a much higher quantity of signal (energy) recorded from EZN (Fig. 3B, first column: contacts A14-15, TB14-15, and B15-16) compared to PZN (Fig. 3B, first column: contact GC6-7). And second, excitability has been systematically overestimated in regions 39 (Right-T1-lateral-anterior) and 51 (Right-Collateral-sulcus), close to EZN but belonging to NIZ.

The second example (Fig. 4) corresponds to patient id013 with EZN = {Left-F1-lateral-prefrontal}, PZN = %. While the estimation of the EZN is correct (region 16) in the control condition Ic_none (Fig. 4A, last column), region 10 (Left-Middle-frontal-sulcus), which is spatially very close to region 16, is identified as being the EZN in
estimated parameters

<figure><img src="image_6.png" /></figure>

GOF

<figure><img src="image_7.png" /></figure>

prediction

<figure><img src="image_8.png" /></figure>

B)

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

C)

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

D)

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 3.** Estimation of EZN for patient id003 and EZN = {Right-T2-Anterior}, PZN = {Right-T2-posterior}. Each row corresponds to the brain network model used to simulate the data: (A) NMM (Ic), (B) NFM with local coupling = 0 (Ic_0), (C) NFM with local coupling = 2 (Ic_2), and (D) NFM with local coupling = 4 (Ic_4). The first column shows the data features extracted from the SEEG simulated observations. The second column shows the predictions of the NMM for the optimization seed with highest goodness of fit (GOF). In the legend, SEEG bipolar contacts are ordered according to their highest amplitude (colors are the same across rows). The third column shows the GOF of each optimization (50 seeds). The GOF is represented as a green dot if it is higher than the 3rd quartile, and in blue otherwise. The last column shows the posterior distributions of the excitability for a selection of nodes (regions). Excitability values used to simulate the data (ground truth) are represented as dots for each node. Colors highlight the ground truth of each node (red for EZN, yellow for PZN, green for NIZ). Please see Supplementary Figure 1 for the estimation of all the nodes from the right hemisphere. The synthetic raw data of condition C are illustrated in Figure 2.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 4.** Estimation of EZN for patient id013 and EZN = {Left-F1-lateral-prefrontal}. Notations are the same as for **Figure 3**. Please see Supplementary Figure 2 for the estimation of all the nodes from the left hemisphere. (A) NNM (lc), (B) NFM with local coupling = 0 (lc_0), (C) NFM with local coupling = 2 (lc_2), and (D) NFM with local coupling = 4 (lc_4).</figcaption></figure>
data simulated with the NFM (Fig. 4B, C and D, last column). As in the previous example, this false positive can be explained by the characteristics of the data features which in this condition are particularly different between NMM and NFM, highlighting the important role of the gain matrix when computing a forward solution. Here, no significant difference is observed in the estimated parameters between conditions lc_0, lc_2, and lc_4 (Fig. 4B, C and D, last column), meaning that the effect of local coupling is limited.

We detailed these two didactic examples to present the main steps of our methodology and also to illustrate two phenomena observed recurrently across simulations. First, estimations from data simulated with NMM were good (Figs. 3A and 4A). And second, the differences observed between NMM and NFM simulated data, mostly due to the different forward models at work, could result in poor estimations of PZN (Fig. 3B) but also of EZN (Fig. 4B, C and D), depending on the spatial configuration.

## 3.2. Group results

Distributions of GOF, precision, recall, and F1 score were computed across all simulations grouped by the 4 types of spatial configuration of the EZN and by the specific brain network model used to simulate the data: lc_none, lc_0, lc_2, and lc_4 (Fig. 5). Overall, regardless of the type of spatial configuration, the medians of the 4 evaluation metrics were systematically higher (or equal) for the control condition lc_none, compared to the 3 other conditions. The first type of spatial configuration (EZN = {1 region}, PZN = {}, Fig. 5A) disclosed the most significant differences across the 4 conditions (Kruskal-Wallis H tests with p-value < 1e-06 for each metric). In particular, GOF was significantly higher for lc_none (median = 0.85) compared to the other conditions (median = 0.50, 0.48, and 0.49 for conditions lc_0, lc_2, and lc_4 respectively; Wilcoxon signed-rank tests with p-value < 1e-10 for each condition). Importantly, even if the difference was also significant for F1 score (Wilcoxon signed-rank tests with p-value < 2e-4 for each condition), the reproducible median value of 1 indicated a good estimation of the EZN for all conditions. Please note that for this configuration, because there is only one region to identify and precision is maximized, recall has to be 1 (this also applies to the next spatial configuration). For the second type of spatial configuration (EZN = {1 half size region}, PZN = {}, Fig. 5B), performances were very similar and only slightly decreased (median GOF = 0.45, 0.43, and 0.43 for conditions lc_0, lc_2, and lc_4 respectively). For the third spatial configuration (EZN = {multiple regions}, PZN = {}, Fig. 5C), the estimation of the EZN was degraded in all conditions as

illustrated by the F1 score (median = 0.55, 0.29, 0.29, and 0.40 for each condition respectively). For the last configuration (EZN = {1 region}, PZN = {1 or 2 regions}, Fig. 5D), the F1 scores indicated a better estimation of the EZN/PZN, compared to the previous configuration (median = 1.00, 0.67, 0.67, and 0.67 for each condition respectively). For the 3rd and 4th types of spatial configurations (Fig. 5C, D), the tendency of having a global decrease of performance for conditions lc_0, lc_2, and lc_4 compared to lc_none was conserved, but no longer significant (Kruskal-Wallis H tests with p-value = 0.1, 0.4, 0.6, and 0.2 for each metric of the 3rd configuration and p-value = 0.01, 0.3, 0.3, and 0.1 for each metric of the 4th configuration).

To summarize, we first observed that overall performances of the NMM inversion were good for data simulated with NMM and NFM for the 1st, 2nd, and 4th spatial configurations and worse for the more complex 3rd spatial configuration (for which an average of 5.3 regions were included in EZN/PZN). And second, the estimations based on NFM simulated data were systematically deteriorated compared to the estimations based on NMM simulated data. This tendency was similar across configurations, suggesting it was independent of the relative complexities of the spatial configurations. In addition, these results effectively disentangled the effects of using a simplified forward model from the effects of absence of local coupling during the NMM inversion. Comparing lc_none with lc_0 showed that the effects of the forward model were quite pronounced while comparing lc_0 with lc_2 and lc_4 showed that the effects of local coupling were limited (Fig. 5).

## 3.3. Effects of differences between NMM and NFM gain matrices

In this section, we studied the differences between the two very different approaches used by the NMM and NFM to compute the forward solution. For the NMM, the activity of a point-like region is projected to each SEEG contact by considering only a (weighted inversed squared) distance between the region and the SEEG contact (Section 2.3.2.2). For the NFM, the activity of each vertex is projected to an SEEG contact by also considering the orientation of the dipole at that vertex, with respect to the SEEG contact (Section 2.3.1.2).

In the previous section, we demonstrated with simulated data that the estimation of the EZN was significantly affected by the simplified gain matrix used by the NMM to compute the forward solution. In fact, the effects due to these different gain matrices are already present and visible in the synthetic SEEG data and the extracted data features. Also, as mentioned in Section 2.4, they were specifically pronounced when comparing conditions
<figure><img src="image_19.png" /><figcaption>**Fig. 5.** Evaluation of the estimation of EZN and PZN based on goodness of fit and classification metrics. Each line corresponds to one type of spatial configuration of the EZN. (A) EZN = {1 region}, PZN = {}, 62 simulations; (B) EZN = {1 half size region}, PZN = {}, 62 simulations; (C) EZN = {multiple regions}, PZN = {}, 11 simulations; and (D) EZN = {1 region}, PZN = {1 or 2 regions}, 11 simulations. Columns represent the median goodness of fit (GOF), the precision, the recall, and the F1 score respectively for each of the 4 conditions lc_none (blue), lc_0 (orange), lc_2 (green), and lc_4 (red). Each violin plot shows the distribution of a given measure across all simulations performed for one type of spatial configuration of the EZN and computed with a specific brain network model.</figcaption></figure>

Ic_none and Ic_0. This is illustrated in Figure 6A where the sensitivity of the SEEG contacts to the activity of EZN and PZN differed between NMM and NFM (different gain matrices applied to an identical source activity): EZN and PZN activities are better detected on contacts B15-16 and FCA13-14 with NMM simulated data, and on contacts A14-15 and GC6-7 with NFM simulated data. Actually, a simple inspection of the gain matrices confirmed these differences. In Figure 6B, we compared the gains of the EZ (top panel) and PZ (bottom panel) regions for NMM (blue) and NFM (orange). For the NFM, these gains were obtained by summing the gains across all vertices belonging to the region of interest. This summation is justified by

the fact that in the condition lc_0 (no local propagation and same initial conditions), all vertices of one region have the exact same activity. With this comparison, we confirmed the discrepancy that we observed in the data features: for NMM, the bipolar contacts most sensitive to the activity of the EZ and the PZ regions were B15-16 and FCA13-14 respectively, while for NFM, these bipolar contacts were A14-15 and GC6-7.

We extended this comparison to the first set of 62 spatial configurations for which EZN = {1 region} and PZN = {} (Table 1) to assess whether the similarity of gain matrices could affect the estimation of the EZN (Fig. 6C). For each simulation, we computed a Pearson correlation
<figure><img src="image_20.png" /><figcaption>**Fig. 6.** Differences between NMM and NFM gain matrices. (A) SEEG data features for patient id003 with EZN = {Right-T2-Anterior}, PZN = {Right-T2-posterior} in the condition lc (left), and lc_0 (right). (B) Normalized gains (vertical axis) of EZN = {Right-T2-Anterior} (top) and PZN = {Right-T2-posterior} (bottom) for SEEG bipolar contacts (horizontal axis) corresponding to NMM (blue) and to NFM (orange). (C) The four evaluation metrics (goodness of fit, precision, recall, and F1 score) plotted against the correlation between gains of NMM and gains of NFM. Each blue dot corresponds to one of the 62 simulations with EZN = {1 region} and PZN = {}, and for the condition lc_0. The orange curve represents the metric averaged over intervals of 0.1 of correlation.</figcaption></figure>

coefficient between the NMM and NFM gains (normalized and rectified) for the EZ region. These coefficients ranged from 0.34 to 0.98 (all p-values for the test of zero correlation were inferior to 1e-5) and were pretty good predictors of the different evaluation measures (goodness of fit and classification metrics). In other words, the more similar the NMM and NFM gain matrices (for EZN), the better the estimation of the EZN (for condition lc_0). On the contrary, an EZN estimation located in a region where gain matrices really differ should be taken more cautiously.

### 3.4. Effects of local coupling

In this section, we analyzed one additional simulation which emphasized the importance of local propagation in seizure spread and the necessity to consider it in the estimation of the EZN. In this simulation, we considered EZN = {Left-Temporal-pole} and PZN = {Left-T2-anterior}, $GC_{scaling} = 3$ and $LC_{scaling} = 4$ (Fig. 7). The key difference with the previous simulations (the 11 spatial configurations for which the PZN was not empty, Table 1) is the cause of the seizure spread from EZN to PZN. In the previous simulations, the spread was mainly due to the global coupling between the EZN and the PZN. Indeed, we adjusted $GC_{scaling}$ separately for each configuration such that even without local coupling (condition lc_0), the seizure effectively propagated from EZN to PZN and

conditions lc_2 and lc_4 corresponded to a supplementary contribution of the local coupling, in addition to the global coupling. In the present simulation, the situation is very different because there is no seizure spread in absence of local coupling (even for higher values of $GC_{scaling} = 3$), which reflects a low coupling strength between EZN and PZN coming from SC, and seizure propagates to PZN only when local coupling is introduced. The same procedure as previously described was used to infer the EZN and PZN in this simulation (Fig. 7). While the 2 regions were clearly identified among others, with a very good accuracy, it should still be noticed that posterior excitability values $u_0(x_i)$ were overestimated (and higher than the critical value -2.05 for region 45 in PZN, Fig. 7F). Indeed, in the absence of global coupling between the EZN and the PZN, and without the possibility to propagate the seizure through local coupling, one hypothesis could be that the NMM used for the estimation of the EZN had no other way to explain the seizure in the PZN than to consider the region part of the EZN (the region could seize autonomously, without the effects of the network). A very recent study (Vattikonda et al., 2023) proposed a new approach for inferring the EZN using an NFM (instead of the NMM we used in the present study) which confirmed this hypothesis. By using the exact same simulation as the one presented in this section, the method demonstrated that as soon as the local
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 7.** Simulation of seizure spread through local coupling. (A) Time courses of the vertices belonging to the EZN = {Left-Temporal-pole} (top row) and to the PZN = {Left-T2-anterior} (bottom row) simulated using an NFM with $GC_{scaling}$ = 3 and $LC_{scaling}$ = 4. (B) Spatial mapping of the onset latencies of the seizure. The seizure starts in the blue region and propagates towards the red region. (C) Time courses of bipolar SEEG contacts. (D) SEEG implantation of the patient. Onset latencies are mapped with arbitrary units (a.u.) on bipolar contacts to show the propagation of the seizure (1: TP'9-10, 167 a.u.; 2: A'13-14, 270 a.u.; 3: B'12-13, 296 a.u.; 4: GPH'14-15, 320 a.u.). (E) Observed and predicted data features of the SEEG bipolar contacts. The 4 bipolar contacts for which observed data features show the highest amplitudes are highlighted in (D). (F) Estimation of the EZN. Notations are the same as for **Figure 3**. The EZN = {Left-Temporal-pole} and the PZN = {Left-T2-anterior} correspond to regions 35 and 45 respectively. Please see Supplementary Figure 3 for the estimation of all the nodes from the left hemisphere. (G) Estimation of the excitability of all vertices belonging to PZN with an NFM (Vattikonda et al., 2023) demonstrating a clear overestimation when local coupling is ignored.</figcaption></figure>

coupling was considered in the estimation procedure, the excitability of the region part of the PZN was no longer overestimated (Fig. 7G).

# 4. DISCUSSION

## 4.1. Summary

This study evaluated the estimation of an EZN in synthetic data with a low spatial resolution NMM, with a particular focus on disentangling the effects of using a simplified forward solution from the effects of absence of

local coupling, two essential modeling features that the NMM approach is intrinsically lacking. For this purpose, we designed an ensemble of spatial configurations of the EZN, typical of what is observed in patients with epilepsy. We simulated these data with an NMM (control condition, lc_none) and with an NFM for increasing strengths of local coupling (lc_0, lc_2, and lc_4). This is to our knowledge the first time that data simulated with NFM are used for an extensive evaluation of inversion procedure based on NMM.

Our results indicated that, whatever the type of spatial configuration, the major difference was found
between condition lc_none and other conditions (lc_0, lc_2, and lc_4). No real effect was observed between the latter conditions, emphasizing the general importance of the forward problem. For the simplest spatial configurations (Fig. 5A, B) for which the EZN was composed of one single region (full and half size), the effect of the forward problem was significant but still limited, with a median F1 score of 1 in all conditions. In such clearly identified clinical situations, it means that the estimation of EZN in empirical data using an NMM inversion can be considered as an efficient option. For more complex spatial configurations (Fig. 5C, D), a comparable tendency was present, not significant this time due to the lower number of configurations (11 < 62, Table 1), in addition to a general and non-negligible decrease of general performance for all conditions, especially more pronounced when the EZN was composed of multiple (not necessarily adjacent and eventually bilateral) regions (median F1 score = 0.55, 0.29, 0.29, and 0.40 for each condition respectively). For these more impacted configurations, which, moreover, correspond to clinically very relevant and more complex cases, inversion based on NMM should be interpreted more cautiously and deeper investigation may be necessary.

## 4.2. Effects of forward gain matrix

The drop of inference performance due to the gain matrix, which was identified in all spatial configurations of the EZN, was also unsurprisingly noticed in the extracted SEEG data features (Fig. 6A). Indeed, many studies already demonstrated the importance of using a high spatial resolution brain network approach, and the ensuing accurate forward gain matrix, for the modeling of SEEG activities. For instance in Cosandier-Rimele et al. (2007a), the authors used a comparable model and found optimal parameters of two volume conductor models (infinite and 3-shell spherical head) to simulate a spike in the left middle temporal gyrus, and mimic an empirical one with high-fidelity. Interestingly, they characterized the attenuation of the electric potential with respect to the distance between the brain source and the SEEG contact and showed that for NFM, this attenuation (hyperbolic, 1/r) was less steep compared to NMM (parabolic, 1/r²). From this difference, they concluded that “neocortical sources of epileptic interictal activity have an extended nature that can be hardly represented by only one equivalent current dipole.” Another noteworthy consequence of using more realistic NFM (and associated gain matrix), and which needs to be considered for an accurate estimation of EZN, is the fact that the SEEG signal may not vary monotonously with respect to the size of the activated cortical patch while a monotonic variation seems

to be more plausible for simulated scalp-EEG signals (Cosandier-Rimele et al., 2007b). Such a very high sensitivity of SEEG to the surrounding geometry of the brain represents a counterpart to its very high spatial specificity.

## 4.3. Effects of local propagation

While the absence of local coupling modeling during the NMM inversion did not really impact the inference performance in general, we showed one specific example where it had a clear and undesired effect, namely an overestimation of the excitability. Very importantly, this may translate in the context of presurgical planning into the misidentification of a candidate region for resection. This indicates that local coupling must be considered very carefully during the inversion procedure. And indeed, preliminary results of such approaches already showed promising results (Fig. 7G).

## 4.4. Limitations

### 4.4.1. Simulated data

The comparison between conditions lc_none and lc_0 was introduced to study the specific effect of the forward gain matrix while conditions lc_2 and lc_4 carried the effects of both the gain matrix and the local coupling. To focus on the specific effects of the local coupling, one option was to modify the NMM used for the inversion by artificially assigning the activity of each region node to the set of vertex nodes associated to the region, and thus enabling the use of an identical gain matrix for both simulation and inversion. We did not follow this option because results showed that the reduction of performance was essentially imputable to the gain matrix and because the main purpose of the study was the evaluation of the NMM per se, not its refinement.

### 4.4.2. Cortical and subcortical regions

In all spatial configurations, we build EZN and PZN only from neocortical regions and did not select subcortical regions, which were systematically considered as point-like structures (NMM) in both the simulation and the inversion procedures. The reason for this is that there is no consensus in the community on how to compute the forward solution for subcortical structures. Given that the effects we observed for neocortical regions should be comparable, but also possibly even more pronounced, for subcortical structures, it remains very important to address such question in the future, and more so as subcortical structures are very often involved in EZN. Our team is currently working in this direction for the particular
case of the hippocampus which, despite its complex
organization, allows a surfacic representation approach
similar to that of neocortical structures. This should enable
its short-term inclusion as a spatially extended brain
structure and result in a more accurate NFM approach.

4.4.3. Forward model

In this study, we only considered the simple infinite
homogeneous model (IHM) for the forward model of
NFM because previous studies (Caune et al., 2014;
Cosandier-Rimele et al., 2006, 2007a) have shown that
when SEEG contacts were “close to the epileptic foci and
relatively far from the skull” (Caune et al., 2014), the IHM
and the one-sphere model provided solutions relatively
similar to the ones of a finite element model (FEM). That
being said, the use of more realistic and advanced SEEG
forward models, for instance based on anisotropic FEM
(Bangera et al., 2010; Medani et al., 2023), needs to be
addressed in the future, and especially in the context of
presurgical planning (Section 4.5). At the same time, we
also should not forget that not only the type of forward
model (boundary element/FEM, isotropic/anisotropic)
plays an important role but the precise localization of the
SEEG electrodes and contacts in anatomical images is
also essential (Zwick et al., 2022).

4.4.4. Inference techniques

In this study, no particularly informative prior was set for
regions of the EZN and PZN and we used a simple infer-
ence technique, that is, the quasi-Newton optimization
algorithm implemented in Stan's probabilistic program-
ming language. Because the algorithm is pretty sensitive
to its random seed (see the multimodal distributions of
GOF on Figs. 3 and 4), we removed the less effective opti-
mizations, probably due to local minima, with a threshold
on the GOF. This resulted in estimates of EZN/PZN with
high confidence, while estimates of NIZ were pretty simi-
lar to their weakly informative prior (Figs. 3 and 4). This in
turn justified the computation of classification metrics
based on point estimates (median of the posterior distri-
butions) rather than full posterior distributions.

It is true that more advanced techniques would defi-
nitely have provided better estimation performance,
and in particular for the configuration EZN = {multiple
regions}. Indeed, well-defined prior combined with
Markov Chain Monte Carlo (MCMC) sampling have
shown to be determinant when fitting both simulated or
empirical data (Hashemi et al., 2020, 2021, 2023; Jha
et al., 2022; Vattikonda et al., 2021; Wang et al., 2023).
In particular, for such fully Bayesian techniques (MCMC,
or Simulation Based Inference) which estimate the full

posterior distributions, there is no additional need to
consider the GOF which is taken into account through
the likelihood term (Bayes theorem, Section 2.5.2). With
that said, Posterior Predictive Check is still one com-
mon diagnostic of Bayesian techniques, similar to the
GOF, which is important to perform to check if the fitted
model is compatible with the observed data.

However, the purpose here was not to look for the best
possible estimation results, but rather to have a simple
framework in order to compare between the different con-
ditions in terms of spatial resolution and to really highlight
the limitations due to the usage of NMM. Also, given the
total number of estimations, the computation time related
to MCMC sampling could have been a limiting factor.

5. CONCLUSIONS AND FUTURE DIRECTIONS

It is fundamental to evaluate the tools and their limitations
for good research in epilepsy (Litt, 2022) and especially
when they are dedicated to a clinical usage. The estima-
tion of EZN using VEP based on NMM already gives very
acceptable results. Still, even in simplest case, the F1
score is not perfect (Fig. 5A) and according to local brain
geometry, some regions may be misclassified. We there-
fore make the hypothesis that an estimation procedure
fully based on NFM defined with a good spatial resolution
(close to the mm² and independent of any parcellation),
embedding a realistic forward model (Section 4.4) and
implementing local propagation, may be required to
obtain more accurate estimations, at least in some spe-
cific configurations. Such a tool will be available in the
near future (Vattikonda et al., 2023). It will enable the
straightforward comparison between NMM and NFM
inversions based on a same shared (simulated and empir-
ical) dataset and provide the ideal and rigorous framework
to draw new conclusions about the relative benefits of
NMM and NFM in estimating the EZN. In addition, flexible
NFM inversion itself will of course offer the opportunity to
compare different forward models as well as bring new
light on the role and importance of local propagation. In
the context of clinical translational research, it will also be
essential in a future study to evaluate how such an NFM
approach could extend the conclusions of previous ret-
rospective studies (Makhalova et al., 2022; Wang et al.,
2023), and especially, in which circumstances it could
help to predict and achieve better surgical outcomes.
Clearly, this step forward will find immediate applications,
such as the development of complementary non-invasive
diagnostic tools, for instance, based on electroencepha-
lography or magnetoencephalography for which an accu-
rate forward problem is also a required component, and
the exploration of innovative therapeutic solutions, such
as, for instance, non-invasive stimulation.
## DATA AND CODE AVAILABILITY

Simulated data and code supporting the findings of this study are available from the authors on reasonable request.

## AUTHOR CONTRIBUTIONS

J.-D.L., H.E.W., and V.J. designed the study. J.-D.L., P.T., and H.E.W. processed the data. J.-D.L. wrote the manuscript with feedback from all authors.

## DECLARATION OF COMPETING INTERESTS

The authors declare no conflict of interest.

## ACKNOWLEDGMENTS

The study was funded through the EU's Horizon Europe Programme SGA 101147319 (EBRAINS 2.0) and Amidex Recherche Blanc, AMX-22-RE-AB-135 (HR-VEP).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00153.

## REFERENCES

Bangera, N. B., Schomer, D. L., Dehghani, N., Ulbert, I., Cash, S., Papavasiliou, S., Eisenberg, S. R., Dale, A. M., & Halgren, E. (2010). Experimental validation of the influence of white matter anisotropy on the intracranial EEG forward solution. *J Comput Neurosci*, 29, 371–387. https://doi.org/10.1007/s10827-009-0205-z

Bartolomei, F., Chauvel, P., & Wendling, F. (2008). Epileptogenicity of brain structures in human temporal lobe epilepsy: A quantified study from intracerebral EEG. *Brain*, 131, 1818–1830. https://doi.org/10.1093/brain/awn111

Bartolomei, F., Lagarde, S., Wendling, F., McGonigal, A., Jursa, V., Guye, M., & Bénar, C. (2017). Defining epileptogenic networks: Contribution of SEEG and signal analysis. *Epilepsia*, 58, 1131–1147. https://doi.org/10.1111/epi.13791

Bartolomei, F., Nica, A., Valenti-Hirsch, M. P., Adam, C., & Denuelle, M. (2018). Interpretation of SEEG recordings. *Neurophysiol Clin*, 48, 53–57. https://doi.org/10.1016/j.neucli.2017.11.010

Bernabei, J. M., Li, A., Revell, A. Y., Smith, R. J., Gunnarsdottir, K. M., Ong, I. Z., Davis, K. A., Sinha, N., Sarma, S., & Litt, B. (2023). Quantitative approaches to guide epilepsy surgery from intracranial EEG. *Brain*, 146, 2248–2258. https://doi.org/10.1093/brain/awad007

Besson, P., Bandt, S. K., Proix, T., Lagarde, S., Jursa, V. K., Ranjeva, J.-P., Bartolomei, F., & Guye, M. (2017). Anatomic consistencies across epilepsies: A stereotactic-EEG informed high-resolution structural connectivity study. *Brain*, 140, 2639–2652. https://doi.org/10.1093/brain/awx181

Cardinale, F., Rizzi, M., Vignati, E., Cossu, M., Castana, L., d'Orio, P., Revay, M., Costanza, M. D., Tassi, L., Mai, R., Sartori, I., Nobili, L., Gozzo, F., Pelliccia, V., Mariani, V., Lo Russo, G., & Francione, S. (2019). Stereoelectroencephalography: Retrospective analysis of 742 procedures in a single centre. *Brain*, 142, 2688–2704. https://doi.org/10.1093/brain/awz196

Caune, V., Ranta, R., Le Cam, S., Hofmanis, J., Maillard, L., Koessler, L., & Louis-Dorr, V. (2014). Evaluating dipolar source localization feasibility from intracerebral SEEG recordings. *Neuroimage*, 98, 118–133. https://doi.org/10.1016/j.neuroimage.2014.04.058

Cosandier-Rimele, D., Badier, J. M., Chauvel, P., & Wendling, F. (2007a). A physiologically plausible spatio-temporal model for EEG signals recorded with intracerebral electrodes in human partial epilepsy. *IEEE Trans Biomed Eng*, 54, 380–388. https://doi.org/10.1109/TBME.2006.890489

Cosandier-Rimele, D., Badier, J. M., Chauvel, P., & Wendling, F. (2007b). Modeling and interpretation of scalp-EEG and depth-EEG signals during interictal activity. In *Presented at the 2007 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, Lyon, France* (pp. 4277–4280). IEEE. https://doi.org/10.1109/IEMBS.2007.4353281

Cosandier-Rimele, D., Badier, J.-M., & Wendling, F. (2006). A realistic spatiotemporal source model for EEG activity: Application to the reconstruction of epileptic depth-EEG signals. In *Presented at the 2006 International Conference of the IEEE Engineering in Medicine and Biology Society, New York, NY* (pp. 4253–4256). IEEE. https://doi.org/10.1109/IEMBS.2006.260678

da Silva, F. L., Blanes, W., Kalitzin, S. N., Parra, J., Suffczynski, P., & Velis, D. N. (2003). Epilepsies as dynamical diseases of brain systems: Basic models of the transition between normal and epileptic activity. *Epilepsia*, 44, 72–83. https://doi.org/10.1111/j.0013-9580.2003.12005.x

Deco, G., Jursa, V. K., Robinson, P. A., Breakspear, M., & Friston, K. (2008). The dynamic brain: From spiking neurons to neural masses and cortical fields. *PLoS Comput Biol*, 4, e1000092. https://doi.org/10.1371/journal.pcbi.1000092

El Houssaini, K., Bernard, C., & Jursa, V. K. (2020). The Epileptor model: A systematic mathematical analysis linked to the dynamics of seizures, refractory status epilepticus, and depolarization block. *eNeuro*, 7, ENEURO.0485-18.2019. https://doi.org/10.1523/ENEURO.0485-18.2019

Fisher, R. S., Acevedo, C., Arzimanoglou, A., Bogacz, A., Cross, J. H., Elger, C. E., Engel, J., Forsgren, L., French, J. A., Glynn, M., Hesdorffer, D. C., Lee, B. I., Mathern, G. W., Moshé, S. L., Perucca, E., Scheffer, I. E., Tomson, T., Watanabe, M., & Wiebe, S. (2014). ILAE official report: A practical clinical definition of epilepsy. *Epilepsia*, 55, 475–482. https://doi.org/10.1111/epi.12550

Goodfellow, M., Rummel, C., Abela, E., Richardson, M. P., Schindler, K., & Terry, J. R. (2016). Estimation of brain network ictogenicity predicts outcome from epilepsy surgery. *Sci Rep*, 6, 29215. https://doi.org/10.1038/srep29215

Hashemi, M., Vattikonda, A. N., Jha, J., Sip, V., Woodman, M. M., Bartolomei, F., & Jursa, V. K. (2023). Amortized Bayesian inference on generative dynamical network models of epilepsy using deep neural density estimators.
Neural Netw, 163, 178-194. https://doi.org/10.1016/j.neunet.2023.03.040

Hashemi, M., Vattikonda, A. N., Sip, V., Diaz-Pier, S., Peyser, A., Wang, H., Guye, M., Bartolomei, F., Woodman, M. M., & Jirsa, V. K. (2021). On the influence of prior information evaluated by fully Bayesian criteria in a personalized whole-brain model of epilepsy spread. *PLoS Comput Biol*, 17, e1009129. https://doi.org/10.1371/journal.pcbi.1009129

Hashemi, M., Vattikonda, A. N., Sip, V., Guye, M., Bartolomei, F., Woodman, M. M., & Jirsa, V. K. (2020). The Bayesian Virtual Epileptic Patient: A probabilistic framework designed to infer the spatial map of epileptogenicity in a personalized large-scale brain model of epilepsy spread. *Neuroimage*, 217, 116839. https://doi.org/10.1016/j.neuroimage.2020.116839

Jehi, L. (2018). The epileptogenic zone: Concept and definition. *Epilepsy Curr*, 18, 12-16. https://doi.org/10.5698/1535-7597.18.1.12

Jha, J., Hashemi, M., Vattikonda, A. N., Wang, H., & Jirsa, V. (2022). Fully Bayesian estimation of virtual brain parameters with self-tuning Hamiltonian Monte Carlo. *Mach Learn Sci Technol*, 3, 035016. https://doi.org/10.1088/2632-2153/ac9037

Jirsa, V., Wang, H., Triebkorn, P., Hashemi, M., Jha, J., Gonzalez-Martinez, J., Guye, M., Makhalova, J., & Bartolomei, F. (2023). Personalised virtual brain models in epilepsy. *Lancet Neurol*, 22, 443-454. https://doi.org/10.1016/S1474-4422(23)00008-X

Jirsa, V. K. (2009). Neural field dynamics with local and global connectivity and time delay. *Philos Trans R Soc Math Phys Eng Sci*, 367, 1131-1143. https://doi.org/10.1098/rsta.2008.0260

Jirsa, V. K., Proix, T., Perdikis, D., Woodman, M. M., Wang, H., Gonzalez-Martinez, J., Bernard, C., Bénar, C., Guye, M., Chauvel, P., & Bartolomei, F. (2017). The Virtual Epileptic Patient: Individualized whole-brain models of epilepsy spread. *Neuroimage*, 145, 377-388. https://doi.org/10.1016/j.neuroimage.2016.04.049

Jirsa, V. K., Stacey, W. C., Quilichini, P. P., Ivanov, A. I., & Bernard, C. (2014). On the nature of seizure dynamics. *Brain*, 137, 2210-2230. https://doi.org/10.1093/brain/awu133

Kalitzin, S., Petkov, G., Suffczynski, P., Grigorovsky, V., Bardakjian, B. L., Lopes da Silva, F., & Carlen, P. L. (2019). Epilepsy as a manifestation of a multistate network of oscillatory systems. *Neurobiol Dis*, 130, 104488. https://doi.org/10.1016/j.nbd.2019.104488

Kini, L. G., Bernabei, J. M., Mikhail, F., Hadar, P., Shah, P., Khambhati, A. N., Oechsel, K., Archer, R., Boccanfuso, J., Conrad, E., Shinohara, R. T., Stein, J. M., Das, S., Kheder, A., Lucas, T. H., Davis, K. A., Bassett, D. S., & Litt, B. (2019). Virtual resection predicts surgical outcome for drug-resistant epilepsy. *Brain*, 142, 3892-3905. https://doi.org/10.1093/brain/awz303

Lagarde, S., Buzori, S., Trebuchon, A., Carron, R., Scavarda, D., Milh, M., McGonigal, A., & Bartolomei, F. (2019). The repertoire of seizure onset patterns in human focal epilepsies: Determinants and prognostic values. *Epilepsia*, 60, 85-95. https://doi.org/10.1111/epi.14604

Liou, J., Smith, E. H., Bateman, L. M., Bruce, S. L., McKhann, G. M., Goodman, R. R., Emerson, R. G., Schevon, C. A., & Abbott, L. (2020). A model for focal seizure onset, propagation, evolution, and progression. *eLife*, 9, e50927. https://doi.org/10.7554/eLife.50927

Litt, B. (2022). Engineers drive new directions in translational epilepsy research. *Brain*, 145, 3725-3726. https://doi.org/10.1093/brain/awac375

Lüders, H. O., Najm, I., Nair, D., Widdess-Walsh, P., & Bingman, W. (2006). The epileptogenic zone: General principles. *Epileptic Disord*, 8(Suppl. 2), S1-S9. https://doi.org/10.1684/j.1950-6945.2006.tb00204.x

Lytton, W. W. (2008). Computer modelling of epilepsy. *Nat Rev Neurosci*, 9, 626-637. https://doi.org/10.1038/nrn2416

Makhalova, J., Medina Villalon, S., Wang, H., Giusiano, B., Woodman, M., Bénar, C., Guye, M., Jirsa, V., & Bartolomei, F. (2022). Virtual epileptic patient brain modeling: Relationships with seizure onset and surgical outcome. *Epilepsia*, 63, 1942-1955. https://doi.org/10.1111/epi.17310

Martinet, L.-E., Fiddyment, G., Madsen, J. R., Eskandar, E. N., Truccolo, W., Eden, U. T., Cash, S. S., & Kramer, M. A. (2017). Human seizures couple across spatial scales through travelling wave dynamics. *Nat Commun*, 8, 14896. https://doi.org/10.1038/ncomms14896

Medani, T., Garcia-Prieto, J., Tadel, F., Antonakakis, M., Erdbrügger, T., Höltershinken, M., Mead, W., Schrader, S., Joshi, A., Engwer, C., Wolters, C. H., Mosher, J. C., & Leahy, R. M. (2023). Brainstorm-DUNEuro: An integrated and user-friendly finite element method for modeling electromagnetic brain activity. *Neuroimage*, 267, 119851. https://doi.org/10.1016/j.neuroimage.2022.119851

Medina Villalon, S., Paz, R., Roehri, N., Lagarde, S., Pizzo, F., Colombet, B., Bartolomei, F., Carron, R., & Bénar, C.-G. (2018). EpiTools, A software suite for presurgical brain mapping in epilepsy: Intracerebral EEG. *J Neurosci Methods*, 303, 7-15. https://doi.org/10.1016/j.jneumeth.2018.03.018

Mercier, M. R., Dubarry, A.-S., Tadel, F., Avanzini, P., Axmacher, N., Cellier, D., Vecchio, M. D., Hamilton, L. S., Hermes, D., Kahana, M. J., Knight, R. T., Llorens, A., Megevand, P., Melloni, L., Miller, K. J., Piai, V., Puce, A., Ramsey, N. F., Schwiedrzik, C. M., ... Oostenveld, R. (2022). Advances in human intracranial electroencephalography research, guidelines and good practices. *Neuroimage*, 260, 119438. https://doi.org/10.1016/j.neuroimage.2022.119438

Nocedal, J., & Wright, S. J. (2006). *Numerical optimization*. Springer series in operations research (2nd ed.). Springer.

Noebels, J. L., Avoli, M., Rogawski, M. A., Olsen, R. W., & Delgado-Escueta, A. V. (Eds.). (2012). *Jasper's Basic Mechanisms of the Epilepsies* (4th ed.). National Center for Biotechnology Information (US). https://doi.org/10.1093/med/9780199746545.001.0001

Olmi, S., Petkoski, S., Guye, M., Bartolomei, F., & Jirsa, V. (2019). Controlling seizure propagation in large-scale brain networks. *PLoS Comput Biol*, 15, e1006805. https://doi.org/10.1371/journal.pcbi.1006805

Proix, T., Bartolomei, F., Chauvel, P., Bernard, C., & Jirsa, V. K. (2014). Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy. *J Neurosci*, 34, 15009-15021. https://doi.org/10.1523/JNEUROSCI.1570-14.2014

Proix, T., Bartolomei, F., Guye, M., & Jirsa, V. K. (2017). Individual brain structure and modelling predict seizure propagation. *Brain*, 140, 641-654. https://doi.org/10.1093/brain/awx004

Proix, T., Jirsa, V. K., Bartolomei, F., Guye, M., & Truccolo, W. (2018). Predicting the spatiotemporal diversity of seizure propagation and termination in human focal epilepsy. *Nat Commun*, 9, 1088. https://doi.org/10.1038/s41467-018-02973-y

Rosenow, F., & Lüders, H. (2001). Presurgical evaluation of epilepsy. *Brain*, 124, 1683-1700. https://doi.org/10.1093/brain/124.9.1683
Saggio, M. L., Crisp, D., Scott, J. M., Karoly, P., Kuhlmann, L., Nakatani, M., Murai, T., Dämpelmann, M., Schulze-Bonhage, A., Ikeda, A., Cook, M., Gliske, S. V., Lin, J., Bernard, C., Jirsa, V., & Stacey, W. C. (2020). A taxonomy of seizure dynamotypes. *eLife*, 9, e55632. https://doi.org/10.7554/eLife.55632

Sarvas, J. (1987). Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem. *Phys Med Biol*, 32, 11–22. https://doi.org/10.1088/0031-9155/32/1/004

Schlafly, E. D., Marshall, F. A., Merricks, E. M., Eden, U. T., Cash, S. S., Schevon, C. A., & Kramer, M. A. (2022). Multiple sources of fast traveling waves during human seizures: Resolving a controversy. *J Neurosci*, 42, 6966–6982. https://doi.org/10.1523/JNEUROSCI.0338-22.2022

Sinha, N., Dauwels, J., Kaiser, M., Cash, S. S., Brandon Westover, M., Wang, Y., & Taylor, P. N. (2017). Predicting neurosurgical outcomes in focal epilepsy patients using computational modelling. *Brain*, 140, 319–332. https://doi.org/10.1093/brain/aww299

Spiegler, A., & Jirsa, V. (2013). Systematic approximations of neural fields through networks of neural masses in the virtual brain. *Neuroimage*, 83, 704–725. https://doi.org/10.1016/j.neuroimage.2013.06.018

Vattikonda, A. N., Hashemi, M., Sip, V., Woodman, M. M., Bartolomei, F., & Jirsa, V. K. (2021). Identifying spatio-temporal seizure propagation patterns in epilepsy using Bayesian inference. *Commun Biol*, 4, 1244. https://doi.org/10.1038/s42003-021-02751-5

Vattikonda, A. N., Woodman, M. M., Lemarechal, J.-D., Daini, D., Hashemi, M., Bartolomei, F., & Jirsa, V. (2023).

Improving epileptogenic zone estimation using Bayesian inference on neural field models. *Health Informatics*. https://doi.org/10.1101/2023.10.01.23296377

Wang, H. E., Scholly, J., Triebkorn, P., Sip, V., Medina Villalon, S., Woodman, M. M., Le Troter, A., Guye, M., Bartolomei, F., & Jirsa, V. (2021). VEP atlas: An anatomic and functional human brain atlas dedicated to epilepsy patients. *J Neurosci Methods*, 348, 108983. https://doi.org/10.1016/j.jneumeth.2020.108983

Wang, H. E., Woodman, M., Triebkorn, P., Lemarechal, J.-D., Jha, J., Dollomaja, B., Vattikonda, A. N., Sip, V., Medina Villalon, S., Hashemi, M., Guye, M., Makhalova, J., Bartolomei, F., & Jirsa, V. (2023). Delineating epileptogenic networks using brain imaging data and personalized modeling in drug-resistant epilepsy. *Sci Transl Med*, 15, eabp8982. https://doi.org/10.1126/scitranslmed.abp8982

Wendling, F., Benquet, P., Bartolomei, F., & Jirsa, V. (2016). Computational models of epileptiform activity. *J Neurosci Methods*, 260, 233–251. https://doi.org/10.1016/j.jneumeth.2015.03.027

Wendling, F., Chauvel, P., Biraben, A., & Bartolomei, F. (2010). From intracerebral EEG signals to brain connectivity: Identification of epileptogenic networks in partial epilepsy. *Front Syst Neurosci*, 4, 154. https://doi.org/10.3389/fnsys.2010.00154

Zwick, B. F., Bourantas, G. C., Safdar, S., Joldes, G. R., Hyde, D. E., Warfield, S. K., Wittek, A., & Miller, K. (2022). Patient-specific solution of the electrocorticography forward problem in deforming brain. *Neuroimage*, 263, 119649. https://doi.org/10.1016/j.neuroimage.2022.119649