domain: scco.neur info scco sdv

text:
M Hashemi
email: meysam.hashemi@univ-amu.fr

A N N Vattikonda

V Sip

Maxime Guye

Fabrice Bartolomei

M M M Woodman

Viktor K Jirsa
email: viktor.jirsa@univ-amu.fr

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

Keywords: Bayesian Inference, Personalized Brain Network Model, Epileptic Seizures, Epileptogenicity BVEP, Bayesian virtual epileptic patient, VEP, virtual epileptic patient, TVB, The Virtual Brain, EZ, epileptogenic zone, PZ, propagation zone, HZ, healthy zone, PPL, probabilistic programming language, MCMC, Monte Carlo Markov Chain, VI, variational inference, NUTS, No-U-Turn Sampler, ADVI, automatic differentiation variational inference, HMC, Hamiltonian Monte Carlo, CV, cross validation, ELBO, evidence lower bound

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Model inversion i.e., finding a set of model parameters that yields the best possible fit to the observed data is a challenging task in statistical inference [START_REF] Gelman | Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science[END_REF]Aster et al., 2013;[START_REF] Groetsch | Inverse Problems: Activities for Undergraduates[END_REF]. Bayesian frameworks offer powerful and principled methods for parameter inference and model prediction from experimental data with a broad range of applications [START_REF] Gelman | Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science[END_REF][START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. Within neuroimaging context, the Bayesian approaches have been widely used for inference of neuronal population's intrinsic parameters and/or interactions between neuronal populations (effective connectivity) in a pre-specified neuronal network from neurophysiological data [START_REF] Friston | Dynamic causal modelling[END_REF][START_REF] David | Dynamic causal modeling of evoked responses in EEG and MEG[END_REF][START_REF] Daunizeau | Stochastic dynamic causal modelling of fmri data: Should we care about neural noise?[END_REF][START_REF] Boly | Connectivity changes underlying spectral EEG changes during propofol-induced loss of consciousness[END_REF]Friston et al., 2014a;[START_REF] Frassle | A generative model of whole-brain effective connectivity[END_REF]. The sampling from the posterior distribution provides the neuroscientist with an estimation of parameter distributions, which make the link to the underlying physiological mechanisms employed in the network model and, thus, enable its scientific interrogation, for instance for resting state and task paradigms [START_REF] David | Dynamic causal modeling of evoked responses in EEG and MEG[END_REF]Friston et al., 2014a), and clinical translation, for instance in epilepsy [START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF].

Bayesian approaches require to calculate conditional probabilities over multiple parameters given observations that is typically interactable analytically. In the field of machine learning, several approaches exist to perform Bayesian inference. Markov Chain Monte Carlo (MCMC; [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF][START_REF] Gelman | Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science[END_REF]) methods allow us to approximate the expectation values of relevant functions by drawing samples from target probability distributions [START_REF] Brooks | Handbook of Markov Chain Monte Carlo[END_REF]Betancourt, 2014b). Based on simulating artificial dataset conditioned on the sampled parameters, several likelihood-free inference methods have been proposed to perform inference when likelihood computations are prohibitively expensive [START_REF] Papamakarios | Fast -free inference of simulation models with Bayesian conditional density estimation[END_REF]Meeds andWelling, 2014, 2015;[START_REF] Gutmann | Bayesian optimization for likelihoodfree inference of simulator-based statistical models[END_REF].

It is well-known that gradient-free sampling algorithms such as Metropolis-Hastings (Metropolis et al., 1953;[START_REF] Hastings | Monte carlo sampling methods using markov chains and their applications[END_REF], Gibbs sampling [START_REF] Geman | Stochastic relaxation, gibbs distributions, and the Bayesian restoration of images[END_REF][START_REF] Gilks | Adaptive rejection metropolis sampling within gibbs sampling[END_REF] and slice-sampling [START_REF] Neal | Slice sampling[END_REF][START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF] generally fail to explore the parameter space efficiently when applied to large-scale inverse problems [START_REF] Girolami | Riemann manifold langevin and hamiltonian monte carlo methods[END_REF][START_REF] Hoffman | The No-U-turn Sampler: Adaptively setting path lengths in hamiltonian monte carlo[END_REF], as often encountered in the application of whole-brain imaging for clinical diagnoses. In particular, the traditional MCMC mix poorly in high-dimensional parameter spaces involving correlated variables [START_REF] Girolami | Riemann manifold langevin and hamiltonian monte carlo methods[END_REF][START_REF] Hoffman | The No-U-turn Sampler: Adaptively setting path lengths in hamiltonian monte carlo[END_REF]Betancourt, 2014b). In contrast, gradient-based algorithms such as Hamiltonian Monte Carlo (HMC; [START_REF] Duane | Hybrid monte carlo[END_REF][START_REF] Neal | MCMC using Hamiltonian dynamics[END_REF]), although computationally expensive, they are far superior to gradient-free sampling algorithms in terms of the number of independent samples produced per unit computational time [START_REF] Sengupta | Gradient-free MCMC methods for dynamic causal modelling[END_REF][START_REF] Sengupta | Gradient-based MCMC samplers for dynamic causal modelling[END_REF]. This class of sampling algorithms provides efficient convergence and exploration of parameter space even in very high-dimensional spaces that may exhibit strong correlations [START_REF] Girolami | Riemann manifold langevin and hamiltonian monte carlo methods[END_REF][START_REF] Hoffman | The No-U-turn Sampler: Adaptively setting path lengths in hamiltonian monte carlo[END_REF]Betancourt, 2014b). Nevertheless, the efficiency of gradient-based sampling methods such as HMC is highly sensitive to the user-specified algorithm parameters [START_REF] Girolami | Riemann manifold langevin and hamiltonian monte carlo methods[END_REF][START_REF] Hoffman | The No-U-turn Sampler: Adaptively setting path lengths in hamiltonian monte carlo[END_REF]. More advanced MCMC sampling algorithms such as No-U-Turn Sampler (NUTS; [START_REF] Hoffman | The No-U-turn Sampler: Adaptively setting path lengths in hamiltonian monte carlo[END_REF]), a self-tuning variant of HMC [START_REF] Carpenter | Stan: A probabilistic programming language[END_REF] solve these issues by adaptively tuning the algorithm parameters. It has been shown that these algorithms efficiently sample from high-dimensional target distributions that allow us to solve complex inverse problems conditioned on massive data set as the observation [START_REF] Gopalan | Efficient discovery of overlapping communities in massive networks[END_REF][START_REF] Gopalan | Scaling probabilistic models of genetic variation to millions of humans[END_REF][START_REF] Kucukelbir | Automatic differentiation variational inference[END_REF].

MCMC has the advantage of being non-parametric and asymptotically exact in the limit of long/infinite runs [START_REF] Gelman | Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science[END_REF]. Of the other alternatives, Variational Inference (VI; [START_REF] Jordan | An introduction to variational methods for graphical models[END_REF][START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF]) turns the Bayesian inference into an optimization problem, which typically results in much faster computation than MCMC methods [START_REF] Gelman | Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science[END_REF][START_REF] Kucukelbir | Automatic differentiation variational inference[END_REF]. However, the classical derivation of VI requires a major model-specific work on defining a variational family appropriate to the probabilistic model, computing the corresponding objective function, computing gradients, and running a gradient-based optimization algorithm [START_REF] Kucukelbir | Automatic variational inference in stan[END_REF][START_REF] Blei | Variational inference: A review for statisticians[END_REF]. Automatic Differentiation Variational Inference (ADVI; [START_REF] Kucukelbir | Automatic differentiation variational inference[END_REF][START_REF] Blei | Variational inference: A review for statisticians[END_REF]) solves these problems automatically.

Probabilistic programming languages (PPLs; [START_REF] Carpenter | Stan: A probabilistic programming language[END_REF][START_REF] Salvatier | Probabilistic programming in Python using PyMC3[END_REF][START_REF] Tran | Edward: A library for probabilistic modeling, inference, and criticism[END_REF][START_REF] Bingham | Pyro: Deep Universal Probabilistic Programming[END_REF]) provide efficient implementation for automatic Bayesian inference on user-defined probabilistic models by featuring the next generation of MCMC sampling and VI algorithms such as NUTS and ADVI, respectively. With the help of PPLs, these algorithms take the advantage of automatic differentiation methods for the computation of derivatives in computer programs to avoid the random walk behavior and sensitivity to correlated parameters [START_REF] Carpenter | Stan: A probabilistic programming language[END_REF][START_REF] Salvatier | Probabilistic programming in Python using PyMC3[END_REF]. In particular, Stan (Stan Development Team, 2018) and PyMC3 [START_REF] Salvatier | Probabilistic programming in Python using PyMC3[END_REF] are high-level statistical modeling tools for Bayesian inference and probabilistic machine learning, which provide the advanced inference algorithms such as NUTS and ADVI, enriched with exten-sive and reliable diagnostics. Although PPLs allow for automatic inference, the performance of these algorithms can be sensitive to the form of parameterization [START_REF] Betancourt | Hamiltonian Monte Carlo for hierarchical models[END_REF]Betancourt, 2014b;[START_REF] Monnahan | Faster estimation of bayesian models in ecology using hamiltonian monte carlo[END_REF]). An appropriate form of reparameterization in the probabilistic models to improve the inference efficiency of system dynamics (governed by a set of nonlinear stochastic differential equations) remains a challenging problem.

On the other hand, due to the potential to improve medical treatment strategies, the personalized large-scale brain network modeling has gained popularity over the recent years [START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF][START_REF] Bansal | Personalized brain network models for assessing structure-function relationships[END_REF]. In the individualized whole-brain modeling approach, the patient-specific information such as anatomical connectivity obtained from non-invasive imaging techniques is combined with the mean-field models of local neuronal activity to simulate the individual's spatio-temporal brain activity at the macroscopic scale [START_REF] Bernard | Virtual brain for neurological disease modeling[END_REF][START_REF] Proix | Individual brain structure and modelling predict seizure propagation[END_REF]. The Virtual Brain (TVB; [START_REF] Sanz Leon | The Virtual Brain: a simulator of primate brain network dynamics[END_REF]) is an open-access computational framework written in Python to reproduce and evaluate the personalized configurations of the brain by using individual subject data. This neuroinformatics platform integrates brain computational modeling and multimodal neuroimaging data to systematically simulate the individual's spatiotemporal brain activity. However, there is currently no specific workflow for automatic model inversion and data fitting validation in preparation for TVB.

More recently, [START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF] [START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF]. However, the inverse problem of such large-scale brain network models is a challenging task due to the intrinsic non-linear dynamics of each brain network node as well as the related large number of model parameters and the observation as commonly encountered in brain-imaging setting.

The aforementioned findings have motivated us to establish a useful link between the most popular probabilistic programming tools (e.g., Stan/PyMC3) and the personalized brain network modeling (e.g., the VEP model), in or-der to systematically predict the location of seizure initiation in a virtual epileptic patient. In the present note, we show how to build the Bayesian Virtual Epileptic Patient (BVEP) as a probabilistic framework designed to infer the hidden/unobserved dynamics of personalized large-scale brain model of epilepsy spread generated by TVB. We first describe the steps necessary to build the BVEP, and then we use NUTS and ADVI algorithms to infer the degree of epileptogenicity of the brain network regions.

Materials and methods

The body of the work is based on the personalized brain network modeling and Bayesian inference as schematically illustrated in Fig. 1. The workflow to build the BVEP consists of two main steps: constructing the VEP, a personalized brain network model of epilepsy spread [START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF], and then embedding the VEP model in a Bayesian framework to infer and validate the model parameters. Following the VEP formulation in state-space representation, the probabilistic reparameterization of the system dynamics is demonstrated. We show that the proposed probabilistic reparameterization in BVEP is able to efficiently invert the nonlinear state-space equations to infer the system dynamics. This novel approach allows us to accurately estimate the spatial map of epileptogenicity in a personalized brain network model of epilepsy spread by taking advantage of PPLs. In this note, we use TVB (see https://www.thevirtualbrain.org) for brain network simulations, and Stan (see https://mc-stan.org) as well as PyMC3 (see https://docs.pymc.io) for inverting the simulated whole-brain model.

In what follows, we show step by step how to build the BVEP model for a particular patient in order to fit the constructed brain model against insilico data and validate our inference. The accuracy and the reliability of the estimations are validated by several convergence diagnostics and posterior behavior analysis.

Individual patient data

For this study, we selected two patients: a 23 year-old female with drugresistant occipital lobe epilepsy (patient 1), and a 24 year-old female with drug-resistant temporo-frontal lobe epilepsy (patient 2). The patients underwent standard clinical evaluation, details of which were described in a previous study [START_REF] Proix | Individual brain structure and modelling predict seizure propagation[END_REF]. The evaluation included non-invasive T1-weighted imaging (MPRAGE sequence, repetition time = 1900 ms, echo time = 2.19 ms, 1.0 x 1.0 x 1.0 mm, 208 slices) and diffusion MRI images (DTI-MR sequence, angular gradient set of 64 directions, repetition time = 10.7 s, echo time = 95 ms, 2.0 x 2.0 x 2.0 mm, 70 slices, b-weighting of 1000 s mm -2 ). The images were acquired on a Siemens Magnetom Verio 3T MR-scanner.

Network anatomy

The structural connectome was built with a reconstruction pipeline using generally available neuroimaging software (see Fig. S1). The current version of the pipeline evolved from a previously described version [START_REF] Proix | Individual brain structure and modelling predict seizure propagation[END_REF].

First, the command recon-all from Freesurfer package [START_REF] Fischl | Freesurfer[END_REF] in version v6.0.0 was used to reconstruct and parcellate the brain anatomy from T1-weighted images. Then, the T1-weighted images were coregistered with the diffusion weighted images by the linear registration tool flirt [START_REF] Jenkinson | Improved optimization for the robust and accurate linear registration and motion correction of brain images[END_REF] from FSL package in version 6.0 using the correlation ratio cost function with 12 degrees of freedom.

The MRtrix package in version 0.3.15 was then used for the tractography.

The fibre orientation distributions were estimated from DWI using spherical deconvolution [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution[END_REF] by the dwi2fod tool with the response function estimated by the dwi2response tool using the tournier algorithm [START_REF] Tournier | Determination of the appropriate b value and number of gradient directions for high-angularresolution diffusion-weighted imaging[END_REF]. Next, we used the tckgen tool, employing the probabilistic tractography algorithm iFOD2 [START_REF] Tournier | Improved probabilistic streamlines tractography by 2nd order integration over fibre orientation distributions[END_REF] to generate 15 millions fiber tracts. Finally, the connectome matrix was built by the tck2connectome tool using the Desikan-Killiany parcellation [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF] generated by FreeSurfer in the previous step (see Table S1 for label names and indices of sub-divided brain regions).

The connectome was normalized so that the maximum value is equal to one.

Network model

Typically, to build a personalized brain network model, the brain regions are defined using a parcellation scheme and a set of mathematical equations is used to model the regional brain activity [START_REF] Sanz-Leon | Mathematical framework for large-scale brain network modeling in the virtual brain[END_REF][START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF]. Taking such a data-driven approach to incorporate the subject-specific brain's anatomical information, the network edges are then represented by structural connectivity of the brain, which are obtained from non-invasive imaging data of individual patients [START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF][START_REF] Bansal | Personalized brain network models for assessing structure-function relationships[END_REF]. In VEP model, the dynamics of brain network nodes are governed by Epileptor equations [START_REF] Jirsa | On the nature of seizure dynamics[END_REF] that are coupled through the structural connectivity matrix derived from diffusion-weighted MRI (dMRI)

techniques [START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF].

The Epileptor is a dynamical model of seizure evolution and is able to realistically reproduce the dynamics of onset, progression and offset of seizurelike events [START_REF] Jirsa | On the nature of seizure dynamics[END_REF][START_REF] Proix | Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy[END_REF]. The Epileptor comprises five state variables coupling two oscillatory dynamical systems on three different timescales: on the fastest timescale, variables x 1 and y 1 account for fast discharges during the ictal seizure states. On the intermediate timescale, variables x 2 and y 2 represent the slow spike-and-wave oscillations. On the slowest timescale, the permittivity state variable z is responsible for the transition between interictal and ictal states. In addition, the interictal and preictal spikes are generated via the term g(x 1 ).

Following [START_REF] Jirsa | On the nature of seizure dynamics[END_REF], the dynamics of full Epileptor model is described by:

ẋ1 = y 1 -f 1 (x 1 , x 2 ) -z + I 1 ẏ1 = 1 τ 1 (1 -5x 2 1 -y 1 ) ż = 1 τ 0 (4(x 1 -η) -z) (1) ẋ2 = -y 2 + x 2 -x 3 2 + I 2 + 0.002g(x 1 ) -0.3(z -3.5) ẏ2 = 1 τ 2 (-y 2 + f 2 (x 2 ))
where

f 1 (x 1 , x 2 ) = x 3 1 -3x 2 1 if x 1 < 0 (x 2 -0.6(z -4) 2 )x 1 if x 1 ≥ 0 f 2 (x 2 ) = 0 if x 2 < -0.25 6(x 2 + 0.25) if x 2 ≥ -0.25 g(x 1 ) = t -t 0 exp -γ(t-τ ) x 1 (τ )dt,
with τ 0 = 2857, τ 1 = 1, τ 2 = 10, I 1 = 3.1, I 2 = 0.45, and γ = 0.01. The degree of epileptogenicity is represented through the value of excitability parameter η. If η > η c , where η c is the critical value of epileptogenicity, Epileptor shows seizure activity autonomously and is referred to as epileptogenic; otherwise Epileptor is in its (healthy) equilibrium state and does not trigger seizures autonomously.

Following [START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF], the full VEP brain model equations (Ncoupled Epileptors) read as follows:

ẋ 1,i = y 1,i -f 1 (x 1,i , x 2,i ) -z i + I 1 ẏ 1,i = 1 τ 1 (1 -5x 2 1,i -y 1,i ) żi = 1 τ 0 (4(x 1,i -η i ) -z i -K N j=1 C ij (x 1,j -x 1,i )) ẋ 2,i = -y 2,i + x 2,i -x 3 2,i + I 2 + 0.002g(x 1,i ) -0.3(z i -3.5) ẏ 2,i = 1 τ 2 (-y 2,i + f 2 (x 2,i )) (2)
where the network nodes are coupled by a linear approximation of permittivity coupling through K N j=1 C ij (x 1,j -x 1,i ), which includes a global scaling factor K, and the patient's connectome C ij .

By applying averaging methods, [START_REF] Proix | Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy[END_REF] have shown that the effect of second neuronal ensemble of Epileptor (i.e., the variables x 2 and y 2 ) is negligible by averaging on the coupled Epileptor equations. Then, under time scale separation (τ 0 1), the fast variables (x 1 and y 1 ) rapidly collapse on the slow manifold [START_REF] Mcintosh | The hidden repertoire of brain dynamics and dysfunction[END_REF], whose dynamics is governed by the slow variable z. This approach yields the 2D reduction of VEP model as follows:

ẋ 1,i = 1 -x 3 1,i -2x 2 1,i -z i + I 1,i żi = 1 τ 0 (4(x 1,i -η i ) -z i -K N j=1 C ij (x 1,j -x 1,i )).
(3)

Depending on the value of excitability parameter η, the 2D Epileptor exhibits different stability regimes. The details regarding linear stability analysis and parameter space exploration of 2D Epileptor are provided in [START_REF] Proix | Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy[END_REF]; [START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF]. For η < η c , a trajectory in the phase plane is attracted to the single stable fixed point of the system on the left branch of the cubic x-nullcline. In this regime, the Epileptor is said to be healthy, meaning not triggering epileptic seizure without external input.

As the value of η increases, the z-nullcline moves down and a saddle-node bifurcation occurs at η = η c corresponding to a seizure onset. For η > η c , the system exhibits an unstable fixed point allowing a seizure to happen (the Epileptor is said to be epileptogenic). Isolated nodes display a bifurcation at the crtical value η c = -2.05 [START_REF] Proix | Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy[END_REF][START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF].

In this study, we use the 2D reduction of VEP model for Bayesian inference of spatial map of epileptogenicity to reduce the computational cost associated with the model parameter estimation. The 2D reduction of Epileptor allows for faster inversion while enabling us to predict the envelope of fast discharges during the ictal seizure states (i.e., onset, propagation and offset of seizure patterns) [START_REF] Proix | Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy[END_REF][START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF].

Spatial map of epileptogenicity

The individual structural connectivity imposes strong constraints on the emergent spatiotemporal dynamics of the virtual brain models [START_REF] Deco | Key role of coupling, delay, and noise in resting brain fluctuations[END_REF][START_REF] Deco | Emerging concepts for the dynamical organization of resting-state activity in the brain[END_REF][START_REF] Melozzi | Individual structural features constrain the mouse functional connectome[END_REF]. Further constraints can be established by incorporating explicit hypotheses into the network model. In the case of epilepsy, clinical hypothesis on the location of epileptogenic zone or lesion allows refining the network pathology to better predict seizure initialization and propagation in individual patients.

In the BVEP brain model, each network node can trigger seizures depending on its connectivity and the excitability value. The parameter η controls the tissue excitability, and its spatial distribution is thus the target of parameter fitting. In this study, depending on the excitability value, the different brain regions are classified into three main types:

• Epileptogenic Zone (EZ): if η > η c , the Epileptor can trigger seizures autonomously (brain region responsible for the origin and early organization of the epileptic activity).

• Propagation Zone (PZ): if η c -∆η < η < η c , the Epileptor does not trigger seizures autonomously but they may be recruited during the seizure evolution since their equilibrium state is close to the critical value.

• Healthy Zone (HZ): if η < η c -∆η, the Epileptor does not trigger seizures autonomously.

Based on the above dynamical properties, the spatial map of epileptogenicity across different brain regions comprises the excitability values of EZ (high value of excitability), PZ (smaller excitability values) and all other regions categorized as HZ (not epileptogenic). Note however, that an intermediate excitability value does not guarantee that the seizure recruits this area as part of the propagation zone, because the propagation is also determined by various other factors including connectivity and brain state dependence.

In the BVEP brain model, the clinical hypotheses can be formulated as the prior knowledge on the spatial distribution of excitability parameters. In this study, assuming no clinical hypothesis on a particular brain area, we assign the same prior distribution on the excitability parameter across all brain regions included in the analysis.

Probabilistic model

The key component in constructing a probabilistic brain network model within a Bayesian framework is the generative model. Given a set of observations, the generative model is a probabilistic description of the mechanisms by which observed data are generated through some hidden states and unknown parameters [START_REF] Daunizeau | Variational Bayesian identification and prediction of stochastic nonlinear dynamic causal models[END_REF][START_REF] Daunizeau | VBA: A probabilistic treatment of nonlinear models for neurobiological and behavioural data[END_REF]. Here, the generative model will therefore have a mathematical formulation guided by the dynamical model that describes the evolution of model's state variables, given parameters, over time. This specification is necessary to construct the likelihood function [START_REF] Cooray | Characterising seizures in anti-nmda-receptor encephalitis with dynamic causal modelling[END_REF][START_REF] Hashemi | Optimal model parameter estimation from eeg power spectrum features observed during general anesthesia[END_REF]. The full generative model is then completed by specifying prior beliefs about the possible values of the unknown parameters (Friston et al., 2014b).

The BVEP brain model presented in this study is built upon two main steps. First, the VEP model equation that provides the basic form of the data generative process describing how the epileptic seizures are generated.

Second, the hypothesis formulation on the spatial map of epileptogenicity in the brain as our prior knowledge. The later component informs the model using hypotheses about the spatial distribution of excitability parameter across different brain regions.

The generative model in the BVEP is formulated based on a system of nonlinear stochastic differential equations of the form (so-called state-space representation):

ẋ(t) = f (x(t), u(t), θ) + w(t), x(0) = x t 0 y(t) = h(x(t)) + v(t) (4)
where x(t) ∈ R n is a n-dimensional vector of system's states evolving over time, x t 0 is the initial state vector at time t = 0, θ ∈ R p contains all the unknown parameters of the VEP model, and u(t) stands for the external input. In addition, y(t) ∈ R m denotes the measured data subject to the measurement error v(t). The process (dynamical) noise and the measurement noise denoted by w(t) ∼ N (0, σ 2 ) and v(t) ∼ N (0, σ 2 ), respectively, are assumed to follow a Gaussian distribution with mean zero and variance σ 2 and σ 2 , respectively. The colored and non-Gaussian dynamical noise can be captured in the term w(t), whereas in the presence of multiplicative noise (i.e., the noise whose intensity depends upon the system's state) or multiplicative feedback (the system's state further influences the driving noise intensity), an additional term appears which can lead to qualitatively different solutions [START_REF] Pesce | Stratonovichto-itô transition in noisy systems with multiplicative feedback[END_REF][START_REF] Volpe | Effective drifts in dynamical systems with multiplicative noise: a review of recent progress[END_REF]. Moreover, f (.) is a vector function that describes the dynamical properties of the system and h(.) represents a measurement function. In source localization problem, h(.)

is known as the lead-field matrix [START_REF] Friston | Multiple sparse priors for the M/EEG inverse problem[END_REF][START_REF] Hu | Unified bayesian estimator of eeg reference at infinity: rrest (regularized reference electrode standardization technique)[END_REF][START_REF] Hu | The statistics of EEG unipolar references: Derivations and properties[END_REF].

We note that the current work focuses on the potential brain sources of observed activity to avoid the inevitable inconsistency associated with mapping from source dipoles to the measurements at electrode contacts (i.e., h(.) is a linear function here).

Considering the 2D reduction of VEP model (cf., Eq. ( 3)), then

x(t) = (x 1,1 , z 1 , x 1,2 , z 2 , ..., x 1,N , z N ) ∈ R n , with n = 2N
, where N is equal to the number of brain regions. Accordingly, θ = (x t 0 ,1 , x t 0 ,2 , ..., x t 0 ,N , η 1 , η 2 , ..., η N , K, σ, σ ) ∈ R p , where p = 3N + 3. Using the reconstruction pipeline to virtualize a patient as described in section 2.2, here N = 84.

The state-space representation (cf. Eq. ( 4)) defining the dynamics of hidden states x(t) is incorporated in the BVEP model as state transition probabilities:

T (x(t), x(t + dt)) ∼ N (x(t) + dtf (x(t), u(t), θ), σ 2 ), (5)
where T denotes the transition probability from state x(t) to x(t + dt).

However, the above parameterization referred to as centered parameterization may exhibit a pathological geometry yielding biased estimations (Betancourt, 2014b; [START_REF] Betancourt | Hamiltonian Monte Carlo for hierarchical models[END_REF].

It has been previously shown that a careful choice of reparameterization increases the effective sample size and decreases the divergences, in particular for the regions of extreme curvature [START_REF] Betancourt | Hamiltonian Monte Carlo for hierarchical models[END_REF][START_REF] Monnahan | Faster estimation of bayesian models in ecology using hamiltonian monte carlo[END_REF]. To avoid pathological samples, and therefore, the biased estimations due to strong correlation between parameters in the centered form of parameterization, we take advantage of location-scale transformation [START_REF] Betancourt | Hamiltonian Monte Carlo for hierarchical models[END_REF] to invert the nonlinear state-space equations, which allows us to decorrelate the parameters representing state variables at successive time steps.

A non-centered reparameterization of the above distribution reads as follows:

T (x * (t), x * (t + dt)) ∼ N (0, 1), x(t + dt) = f (x(t), u(t), θ) + dtdx + σx * (t + dt). (6)
In section 3, we show that using the non-centered form of parameterization to infer the system dynamics dramatically improves the performance of sampling by avoiding biased estimations due to the strong correlation between parameter.

Inference/Prediction

A generative model is characterized by the joint probability distribution of the model parameters and the observation P (Y, ϑ) where Y denotes the observed variables, and ϑ includes the system's hidden variables and the model parameters [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF][START_REF] Daunizeau | Variational Bayesian identification and prediction of stochastic nonlinear dynamic causal models[END_REF]. Bayesian techniques infer the distribution of unknown parameters of the underlying data generating process, given only observed responses and prior beliefs about the underlying generative process [START_REF] Gelman | Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science[END_REF][START_REF] Cooray | Dynamic causal modelling of electrographic seizure activity using Bayesian belief updating[END_REF]. By product rule, the generative model can be defined in terms of likelihood and prior on the model parameters, whose product yields the joint density:

P (Y, ϑ) = P (Y | ϑ)P (ϑ), (7)
where prior distribution P (ϑ) includes our prior beliefs about the hidden variables and potential parameter values, while the conditional likelihood term P (Y | ϑ) represents the probability of obtaining the observation, with a given set of parameter values. In Bayesian inference, we seek the posterior density P (ϑ | Y), which is the conditional distribution of model parameters given the observation [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF][START_REF] David | Dynamic causal modeling of evoked responses in EEG and MEG[END_REF]). Bayes's Theorem expresses this posterior density in terms of likelihood and prior as follows:

P (ϑ | Y) = P (Y | ϑ)P (ϑ) P (Y) , (8)
where the denominator P(Y) represents the probability of the data and it is known as evidence or marginal likelihood (in practice amounts to simply a normalization term [START_REF] Gelman | Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science[END_REF]).

To sample from posterior density P (ϑ | Y), the performance of HMC is highly sensitive to the step size and the number of steps in leapfrog integrator for updating the position and momentum variables in Hamiltonian dynamic simulation [START_REF] Hoffman | The No-U-turn Sampler: Adaptively setting path lengths in hamiltonian monte carlo[END_REF]. If the number of steps in the leapfrog integrator is chosen too small, then HMC exhibits an undesirable random walk behaviour similar to Metropolis-Hastings algorithm, and thus algorithm poorly explores the parameter space. If the number of leapfrog steps is chosen too large, the associated Hamiltonian trajectories may loop back to a neighbourhood of the initial state, and the algorithm wastes computation efforts [START_REF] Hoffman | The No-U-turn Sampler: Adaptively setting path lengths in hamiltonian monte carlo[END_REF][START_REF] Betancourt | The geometric foundations of Hamiltonian Monte Carlo[END_REF]. NUTS extends HMC with adaptive tuning of both the step size and the number of steps in leapfrog integration to sample efficiently from posterior distributions [START_REF] Hoffman | The No-U-turn Sampler: Adaptively setting path lengths in hamiltonian monte carlo[END_REF][START_REF] Betancourt | Generalizing the no-u-turn sampler to riemannian manifolds[END_REF][START_REF] Betancourt | The geometric foundations of Hamiltonian Monte Carlo[END_REF]. In an alternative approach, ADVI posits a family of densities, automatically computes the gradients, and then finds the closest member (measured by Kullback-Leibler divergence) to the target distribution [START_REF] Kucukelbir | Automatic differentiation variational inference[END_REF][START_REF] Blei | Variational inference: A review for statisticians[END_REF]. In this study, we use NUTS, a self-tuning variant of HMC, as well as ADVI to approximate the posterior distribution of the model parameters (cf., Eq. ( 3)).

The prior on excitability parameter for all brain regions included in the analysis was assumed as a normal distribution with a mean of -2.5 and a standard deviation of 1.0, i.e., N (-2.5, 1.0). Moreover, we placed a weakly informative prior on the system initial conditions and the global coupling parameter K, as a normal distribution centered at the ground-truth with standard deviation of 1.0 (see Tables S2 andS3). The prior on the hyperparameters was considered as a generic weakly informative prior N (0, 1.0).

After fitting a Bayesian model, it is often necessary to measure the predictive accuracy of the inferred model [START_REF] Gelman | Understanding predictive information criteria for Bayesian models[END_REF]. The information criteria and leave-one-out cross-validation (LOO; (Vehtari et al., 2017b)) are two rigorous approaches to assess the model's ability in prediction of new data. Taking the existing simulation draws from log-likelihood evaluated at the posterior of the parameter values, widely applicable information criterion (WAIC; [START_REF] Watanabe | Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory[END_REF]) and Pareto-smoothed importance sampling (PSIS; (Vehtari et al., 2017a)) LOO allow for efficiently estimating predictive accuracy of a fitted Bayesian model within a negligible computational time relative to the cost of model fitting [START_REF] Vehtari | Bayesian leave-one-out cross-validation approximations for gaussian latent variable models[END_REF].

Inference diagnostics

After running a MCMC sampling algorithm, it is necessary to carry out some statistical analysis in order to evaluate the convergence of MCMC samples [START_REF] Cowles | Markov Chain Monte Carlo convergence diagnostics: A comparative review[END_REF]. One simple way to assess the performance of MCMC algorithms based on posterior samples is to visualize how well the chain is mixing (i.e., MCMC sampler explores all the modes in the parameter space efficiently). This can be monitored in different ways including traceplot (evolution of parameter estimates from MCMC draws over the iterations), pair plots (to identify collinearity between variables), and autocorrelation plot (to measure the degree of correlation between draws of MCMC samples). A more quantitative way to assess the MCMC convergence to the stationary distribution is to estimate the potential scale reduction factor R [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF][START_REF] Brooks | General methods for monitoring convergence of iterative simulations[END_REF], and effective sample size N ef f [START_REF] Gelman | Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science[END_REF] based on the samples of posterior model probabilities (see Appendix for definitions). The R diagnostic provides estimate of how much variance could be reduced by running chains longer. Each MCMC estimation has R statistic associated with it, which is essentially the ratio of between-chain variance to within-chain variance [START_REF] Gelman | Inference from iterative simulation using multiple sequences[END_REF][START_REF] Brooks | General methods for monitoring convergence of iterative simulations[END_REF]. If R is approximately less than 1.1, the MCMC convergence has been achieved (approaching to 1.0 in the case of infinite samples); otherwise, the chains need to be run longer [START_REF] Gelman | Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science[END_REF]. Moreover, the N ef f statistic gives the number of independent samples represented in the chain. The larger the effective sample size, higher the precision of MCMC estimates. Note that these are necessary but not sufficient conditions for convergence of MCMC samples.

In addition to the general MCMC diagnostics mentioned above, the NUTS- Confusion matrix is a metric to evaluate the accuracy of a classification [START_REF] Fawcett | An introduction to ROC analysis[END_REF][START_REF] Powers | Evaluation: From precision, recall and f-measure to roc., informedness, markedness & correlation[END_REF]. The element q i,j is equal to the number of observations known to be in class i but predicted to be in class j, with i, j ∈ {1, 2, ...Q}, where Q is the total number of classes [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. In the BVEP model, we defined three groups namely HZ, PZ, and EZ to classify brain regions, thus Q = 3.

Moreover, in order to quantify the accuracy of the inference, we plot the posterior z-scores (denoted by z) against the posterior shrinkage (denoted by s), which are defined as (Betancourt, 2014a):

z = | θ -θ * σ post |, (9)
s = 1 - σ 2 post σ 2 prior , (10)
where θ and θ * are the estimated-mean and the ground-truth, respectively, whereas σ 2 prior , and σ 2 post indicate the variance (uncertainty) of the prior and the posterior, respectively. The posterior z-score quantifies how much the posterior distribution encompasses the ground-truth, while the posterior shrinkage quantifies how much the posterior distribution contracts from the initial prior distribution (Betancourt, 2014a).

Synthetic data sets and model inversion

In order to validate the inference using BVEP, we take advantage of and MEG with a wide range of clinical applications from Alzheimer disease [START_REF] Zimmermann | Differentiation of alzheimer's disease based on local and global parameters in personalized virtual brain models[END_REF], chronic stroke [START_REF] Falcon | Functional mechanisms of recovery after chronic stroke: Modeling with the virtual brain[END_REF] to human focal epilepsy [START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF].

In this study, we used TVB to reconstruct the personalized brain network model. In order to validate the inference on spatial epileptogenicity, we simulated epileptic seizures for two patients: one simulation with the seizure spread to all brain nodes specified as PZ (patient 1), and another with the seizure spread to some of the brain nodes specified as PZ (patient 2). These data sets were generated using two different structural connectivity matrices and distinct spatial map of epileptogenicity.

The seizure activity of patient 1 was simulated by setting two regions as EZ, and three regions as PZ, where EZ idx ∈ {7, 35}, and PZ idx ∈ {6, 12, 28},

with the excitability values η ez = -1.6, and η pz = -2.4, respectively. All the other brain nodes were fixed as not epileptogenic i.e., HZ with η hz = -3.6

(see Table S2).

To simulate the seizure activity of patient 2, we selected two brain regions as EZ, and five regions as PZ, at the nodes EZ idx ∈ {7, 24}, and PZ idx ∈ {10, 23, 27, 28, 35}, respectively. For the regions selected as EZ, the excitability value was set to η ez = -1.5. The excitability of PZ was set as η pz = -2.6, and all the other regions were defined as HZ with η hz = -3.4

(see Table S3).

In both synthetic data sets, to simulate the VEP model as a system of stochastic differential equations, we used an Euler-Maruyama integration scheme with an integration step of 0.04. The additive white Gaussian noise was introduced in the state variable

x(t) = (x 1,i (t), y 1,i (t), z i (t), x 2,i (t), y 2,i (t))
with zero mean and variance (0.01, 0.01, 0.0, 0.0015, 0.0015). The initial conditions were selected in the interval (-2.0, 5.0) for each state variable.

Results

The result of workflow in the BVEP model to estimate the spatial map of epileptogenicity across different brain regions for patient 1 is illustrated in Furthermore, the NUTS-specific diagnostics were monitored to check whether the Markov chain has converged (see Fig. S7). The diagnostics plot shows that there are no divergent transitions in HMC indicating that the posterior density was explored efficiently. Also, none of the NUTS iterations reached maximum tree-depth (its value to run NUTS was specified 10.0 here) indicating that the optimal number of leapfrog steps needed for the Hamiltonian simulation was sufficiently lower than the maximum. Together, these diagnostics validate that the samples by NUTS has converged to the target distribution. unstable fixed points, respectively. From Fig. 4A and D, it can be observed that a trajectory of an HZ (node number 1) is attracted to the stable fixed point of the system (on the left branch of cubic x-nullcline) meaning not triggering epileptic seizure. For a PZ (node number 6), due to the coupling strength and the value of excitability which is close to the critical value of epileptogenicity, the z-nullcline moves down, causing a bifurcation thereby allowing the seizure to propagate here (see Fig. 4B andE). For the EZ (node number 7), the system exhibits an unstable fixed point due to the high value of excitability. In this regime, Epileptor possesses a limit cycle and the seizure triggers autonomously (see Fig. 4C andF). Notice that the topology of simulated and predicted phase-plane trajectories show very good agreement, except the amplitude of state variable z i , from which the estimation indicates the result of a larger parameter recovery. Note that only the activity of fast variable x 1,i is the target of fitting as the observed data. Fig. 4: Comparison between the simulated (top row) and the predicted (bottom row) phase-plane for different brain node types in the BVEP model. From left to right, the columns correspond to the brain nodes specified as HZ, PZ, and EZ, respectively. A trajectory for these brain regions is shown in green, yellow and red, respectively. In each phase-plane, depending on the excitability parameter, the intersection of x-and znullclines (colored in blue) determines the fixed point of the system. Full circle and empty circle indicate the stable and unstable fixed points, respectively.

To compare the BVEP inversion by NUTS and ADVI schemes, Fig. 5 displays the histogram of MCMC samples and the kernel density estimates of the posteriors generated from NUTS (left panel) versus those obtained by ADVI (right panel). From this figure, we observe that NUTS and ADVI perform similarly in their estimates of the posterior, except that the meanfield ADVI slightly underestimates the variances compared to the estimations by NUTS algorithm. However, taking both approaches, the true values of excitabilities (dashed vertical lines) are under the support of the posterior densities indicating that the parameter recovery was successful. The samples corresponding to the brain nodes specified as HZ, PZ, and EZ are illustrated in green, yellow and red, respectively. Note that the prior for all 84 brain regions included in the analysis was assumed as a normal distribution centered at -2.5 with standard deviation of 1.0 (i.e., N (-2.5, 1.0) as shown in blue). To invert the BVEP model by NUTS algorithm, we used 200 number of sampling iterations and 200 warmup with the expected acceptance probability of 0.95, whereas to run ADVI, the maximum number of iterations and the convergence tolerance were set to 50000 and 0.001, respectively. In terms of computational time, for these algorithm configurations, sampling by NUTS took 23993.5 seconds whereas the running time of ADVI was 5392.62 seconds. Once the model parameters have been estimated, it is necessary to assess the convergence of MCMC samples. To verify the reliability of the inferred estimates, we monitored the potential scale reduction factor R as it is the most reliable quantitative metric for MCMC convergence. In addition,

we have plotted the posterior samples from the joint posterior probability distribution to show the efficiency of the transformed non-centered parameterization in comparison to the centered form of parameterization. Fig. 6 top row represents the posterior samples from the joint posterior probability distribution between the hyper-parameters σ and σ , which are the standard deviation of the process (dynamical) noise and the measured noise, respectively (cf., Eq. ( 4)). In this figure, the left and middle columns show the result of sampling by NUTS with non-centered and centered form of parameterization, respectively. For the sake of comparison with NUTS, the last column illustrates the result from mean-field variant of ADVI. The dots in each scatterplot represent 200 samples drawn from the joint posterior probability distribution. In Fig. 6A and B, it can be clearly seen that there is no correlation between the posterior samples drawn from non-centered parameterization, whereas the samples from the centered form show a high collinearity between hyper-parameters. Such a high collinearity leads to inefficient exploration of posterior which can be quantifiably observed in decreased numbers of effective samples and increased R values. The values of R for all of hidden states and parameters estimated by non-centered form are below 1.05 (see Fig. 6D), whereas more than 82 percent estimations by centered form has R value above 1.1 (see Fig. 6E). This indicates that the Markov chains converged for non-centered but not for centered form of parameterization. We report that the ratio of effective number of samples to the number of iteration (N ef f /N iter ) returned by centered form of NUTS is less than 0.001 for all the estimated parameters. This indicates the poor sampling from centered form of parameterization as it generates a very small number of independent samples per Markov chain.

Moreover, scatterplot of samples drawn from joint posterior probability distribution between the hyper-parameters σ and σ estimated by the mean-field ADVI is illustrated in Fig. 6C. Since by definition, the meanfield variant of ADVI ignores the cross correlation between parameters, the samples drawn using mean-field ADVI show no correlation between hyperparameters. Lastly, in order to check the convergence of ADVI, evidence lower bound (ELBO), the variational objective function, is plotted versus the number of iterations (see Fig. 6F). While the algorithm appears to have converged in 10000 iterations, the algorithm runs for another few thousand iterations to guarantee the convergence until the change in ELBO drops below the tolerance of 0.001.

Discussion

This technical note presents a probabilistic framework namely the Bayesian Virtual Epileptic Patient (BVEP) to infer the spatial map of epileptogenicity for developing a personalized large-scale brain model of epilepsy spread (cf. To demonstrate the potential functionality of the BVEP in prediction of seizure initiation and propagation, we simulated simple and complex seizure spread using different spatial maps of epileptogenicity (cf. Fig. S2). We have used these synthetic data for fitting, since given the ground-truth of model parameters, we can use the standard error metrics such as confusion matrix, posterior shrinkages and posterior z-scores to validate the accuracy of the estimations, thus to evaluate the performance of the proposed approach.

Our results demonstrated that in both synthetic data sets, by inverting the large-scale brain network model with the help of PPLs (Stan/PyMC3), we can achieve a remarkable similarity between the simulated and the predicted seizure activity regarding the initiation, propagation and termination. Although, the simulation was generated by the full VEP model comprising five state variables at each brain node, the 2D reduced variant of the model was still able to successfully predict the key data features such as onset, propagation and offset of seizure patterns, while considerably alleviating the computational time of the Bayesian inference. This 2D reduction is limited to modeling the average of fast discharges during the ictal seizure states [START_REF] Proix | Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy[END_REF], which as we showed (see Figs. In such cases, confusion matrix may result in a high accuracy performance, whereas plotting the posterior z-scores versus the posterior shrinkages is particularly useful for identifying the malfunctioning in the inference such as overfitting, or a poorly-chosen prior that biases the estimations (Betancourt, 2014a).

Understanding brain dynamics in epilepsy is critical for developing therapeutic approaches towards brain interventions to improve the surgical outcome. Using theory of nonlinear dynamic systems, the complete taxonomy of epileptic seizures with a thorough description of bifurcations that give rise to onset, offset and seizure evolution characteristics has been extensively investigated elsewhere [START_REF] Jirsa | On the nature of seizure dynamics[END_REF]. In parameter space description of Epileptor, the seizure onset and offset are described by saddle-node and homoclinic bifurcations, respectively [START_REF] Proix | Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy[END_REF][START_REF] El Houssaini | Seizures, refractory status epilepticus, and depolarization block as endogenous brain activities[END_REF]. The emergent dynamic effects in the BVEP model crucially depend on the interplay between network node model (Epileptor), patient specific structural connectivity (from dMRI), and spatial maps of epileptogenicity (EZ, PZ, HZ). According to the dynamical properties of Epileptor model [START_REF] Jirsa | On the nature of seizure dynamics[END_REF][START_REF] Proix | Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy[END_REF], we classified the brain regions into three main types: EZ (exhibiting unstable fixed point corresponding to the brain area responsible for the seizure initiation), PZ (close to saddle-node bifurcation corresponding to the candidate brain area responsible for the seizure propagation), and HZ (exhibiting stable fixed point corresponding to healthy brain area). This approach allows us to define the spatial map of epileptogenicity based on the excitability parameter value, which is the target of fitting.

It is important to note that an excitability value close to the critical value of epileptogenicity does not guarantee that the seizure originating from pathological brain areas (i.e., responsible for the seizure onset specified as EZ) propagates to such brain regions defined as PZ. By a detailed patient evaluation, it has been reported that the individual structural connectivity is essential for predicting seizure spatial propagation [START_REF] Petkov | A critical role for network structure in seizure onset: A computational modeling approach[END_REF][START_REF] Taylor | The impact of epilepsy surgery on the structural connectome and its relation to outcome[END_REF][START_REF] Proix | Individual brain structure and modelling predict seizure propagation[END_REF][START_REF] Proix | Predicting the spatiotemporal diversity of seizure propagation and termination in human focal epilepsy[END_REF]. However, it has been recently

shown that purely structural information is not sufficient to predict the propagation and eventual stopping of the seizures [START_REF] Olmi | Controlling seizure propagation in large-scale brain networks[END_REF]. Rather, the abnormal activity in the recruited regions is a complex network effect which depends on the interplay between multiple factors including the brain region's epileptogenicity (node dynamics) [START_REF] Bartolomei | Epileptogenicity of brain structures in human temporal lobe epilepsy: a quantified study from intracerebral EEG[END_REF][START_REF] Goodfellow | Estimation of brain network ictogenicity predicts outcome from epilepsy surgery[END_REF][START_REF] Lambert | Brain regions and epileptogenicity influence epileptic interictal spike production and propagation during nrem sleep in comparison with wakefulness[END_REF], the individual structural connectivity (network structure) [START_REF] Jirsa | The virtual epileptic patient: Individualized whole-brain models of epilepsy spread[END_REF][START_REF] Proix | Individual brain structure and modelling predict seizure propagation[END_REF][START_REF] Proix | Predicting the spatiotemporal diversity of seizure propagation and termination in human focal epilepsy[END_REF], and brain state dependence (network dynamics) [START_REF] Spencer | Neural networks in human epilepsy: Evidence of and implications for treatment[END_REF][START_REF] Terry | Seizure generation: The role of nodes and networks[END_REF]Kramer and Cash, 2012;[START_REF] Van Diessen | Functional and structural brain networks in epilepsy: What have we learned?[END_REF][START_REF] Burns | Network dynamics of the brain and influence of the epileptic seizure onset zone[END_REF]. Furthermore, there are nonlinearities and multiple propagation patterns that can be observed for the same excitability parameter sets due to the coupled nonlinear system dynamics (cf., Fig. S2). In this work, the seizure recruitment is characterized by complex spatio-temporal dynamics of large-scale brain network i.e., seizure originates from a local network and recruits candidate brain regions strongly coupled to the pathological areas by perturbing their stable dynamics (if K=0, then there is no seizure recruitment). Among the candidate brain regions for seizure propagation, due to stronger connection to the pathological areas defined as EZ, the node PZ idx = {28} can be recruited by a weak global coupling (see Fig. S2). Rather, a stronger coupling is required for seizure recruitment to all other candidate brain areas. This is in agreement with experimental observations that seizures tend to have a common spatial origin in the same patient [START_REF] Proix | Individual brain structure and modelling predict seizure propagation[END_REF][START_REF] Proix | Predicting the spatiotemporal diversity of seizure propagation and termination in human focal epilepsy[END_REF]. According to this knowledge, we placed a weakly informative prior on the global coupling parameter centered at the ground-truth. Overestimation of global coupling parameter leads to miss-classification of PZ as HZ (see Figs. 2E and3B), whereas an underestimation of coupling may yield to miss-classification of PZ as EZ (see Figs. S4E andS5B). However, the stability analysis of the network dynamics indicates that the seizure propagation can be controlled by an optimal intervention on the structural connectivity matrix [START_REF] Olmi | Controlling seizure propagation in large-scale brain networks[END_REF][START_REF] An | Optimization of surgical intervention outside the epileptogenic zone in the virtual epileptic patient (VEP)[END_REF] implying that the patient-specific network connectivity is predictive for seizure propagation pattern. Therefore, the seizure propagation may not be easily controlled by a simple dissection of the individual nodes [START_REF] Stam | Modern network science of neurological disorders[END_REF][START_REF] Goodfellow | Estimation of brain network ictogenicity predicts outcome from epilepsy surgery[END_REF][START_REF] Taylor | The impact of epilepsy surgery on the structural connectome and its relation to outcome[END_REF], as in the surgical treatment of epilepsy, it has been reported that the resection does not always lead to post-surgery seizure freedom in the brain [START_REF] Tonini | Predictors of epilepsy 48 surgery outcome: a meta-analysis[END_REF][START_REF] De Tisi | The long-term outcome of adult epilepsy surgery, patterns of seizure remission, and relapse: a cohort study[END_REF][START_REF] Mohammed | Impact of epilepsy surgery on seizure control and quality of life: A 26-year follow-up study[END_REF][START_REF] Taylor | The impact of epilepsy surgery on the structural connectome and its relation to outcome[END_REF].

In this study, the analysis of phase-plane trajectories of the observed system versus the prediction was carried out across different brain regions in order to gain a better understanding of the mechanisms underlying seizure initiation and propagation within the proposed approach (cf. Fig. 4). For different brain node types (e.g., EZ, PZ, and HZ), the dynamics of seizure initiation and recruitment in the phase-plane was captured well by the prediction. From inference perspective, we observed good correspondence to the phase portraits of the observed system including equilibria (the intersection of the nullclines), the stability or instability of the equilibria, and the flow of trajectories. These results validate our Bayesian inversion procedure in order to understand the spatio-temporal evolution of seizure activity, paving the way for further studies on possible seizure prevention approaches.

In this technical note, we have used both NUTS and ADVI schemes to infer the spatial map of epileptogenicity in a personalized whole-brain model of epilepsy spread. The results from both inference schemes led to similar estimation of spatial map of epileptogenicity across brain regions, except that the ADVI slightly underestimate the variances compared to the estimations by NUTS algorithm (cf. Fig. 5). The similarity between the inversions using the two schemes indicates that the variational approximation offers an appropriate alternative to NUTS sampling in BVEP model inversion. Our results demonstrated that there is a clear reduction in computational cost in performing the inference by ADVI compared to NUTS (4-5 times faster for the used algorithm configurations), which may be important when applying BVEP approach to large data sets of patient cohort. While it is generally known that ADVI is more computationally appealing than NUTS, however, it can be challenging to discover algorithmic problems with this approximation [START_REF] Yao | Yes, but did it work?: Evaluating variational inference[END_REF]. The convergence of ADVI can be assessed by monitoring the running average of ELBO changes, whereas NUTS is furnished with several general and specific diagnostics to assess whether the Markov chain has converged. In addition, ADVI may get stuck in local minima during gradient descent optimization and its mean-field variant is unable to cover all the modes of the multi-modal posterior densities [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF][START_REF] Blei | Variational inference: A review for statisticians[END_REF][START_REF] Yao | Yes, but did it work?: Evaluating variational inference[END_REF].

Lastly, we investigated the efficiency of transformed non-centered parameterization. In agreement with previous studies showing that NUTS is sensitive to the parameterization [START_REF] Betancourt | Hamiltonian Monte Carlo for hierarchical models[END_REF][START_REF] Monnahan | Faster estimation of bayesian models in ecology using hamiltonian monte carlo[END_REF], our results indicated that the non-centered form of parameterization to invert the nonlinear state-space equations yields an efficient parameter-space exploration, whereas the centered form of sampling demonstrates an inefficient exploration due to the high collinearity between model parameters (see Fig. 6A and D versus Fig. 6B andE). In addition, based on the convergence diagnostics such as R, we demonstrated that the samples generated by NUTS converged faster in the non-centered parameterization compared to the centered form of parameterization.

In this technical note, we proposed a novel approach to build personalized Comparison of the implementations in Stan, PyMC3 and other alternative PPL packages is beyond the scope of this note.

To our best knowledge, this study is the first personalized large-scale brain network modeling approach for inferring the spatial map of epileptogenicity (properties of nodes) based on patient-specific whole-brain anatomical information (i.e., network structure derived from dMRI). Dynamic Causal Modelling (DCM; [START_REF] Friston | Dynamic causal modelling[END_REF]) is a well-established framework for analyzing neuroimaging modalities (such as fMRI, MEG, and EEG) by neural mass models where inferences can be made about the coupling among brain regions (effective connectivity) to infer how the changes in neuronal activity of brain regions are caused by activity in the other regions through the modulation in the latent coupling [START_REF] Friston | Dynamic causal modelling[END_REF][START_REF] David | Dynamic causal modeling of evoked responses in EEG and MEG[END_REF][START_REF] Moran | A neural mass model of spectral responses in electrophysiology[END_REF][START_REF] Boly | Connectivity changes underlying spectral EEG changes during propofol-induced loss of consciousness[END_REF][START_REF] Kiebel | Dynamic causal modeling for EEG and MEG[END_REF]Friston et al., 2014b;[START_REF] Frassle | A generative model of whole-brain effective connectivity[END_REF]. Using DCM, focal seizure activity in electrocorticography (ECoG) data was recently studied to estimate the key synaptic parameters or coupling connections using observed signals in a human subject [START_REF] Papadopoulou | Tracking slow modulations in synaptic gain using dynamic causal modelling: Validation in epilepsy[END_REF]. In another study, Bayesian belief updating scheme for DCM has been used to estimate the synaptic drivers of cortical dynamics during a seizure from EEG/ECoG recordings with a little computational expense [START_REF] Cooray | Dynamic causal modelling of electrographic seizure activity using Bayesian belief updating[END_REF]. Although DCM can be used to model and track the changes in excitatory-inhibitory balance at seizure onset/offset, these studies are based on single neural mass model (i.e., small number of cortical sources are modelled), and the non-linear ordinary differential equation representing the neural mass model is approximated by its linearization, with which only the seizure onset or offset can be modelled but not both. In this note, the Bayesian Virtual Epileptic Patient (BVEP) model can characterize whole-brain spatio-temporal nonlinear dynamics of seizure propagation. This approach allows describing the onset and offset of ictal states as well as the alternation between normal and ictal periods. The BVEP approach relies on the patient-specific structural data rather formulating the inverse problem purely in terms of unknown model parameters used in DCM [START_REF] Papadopoulou | Tracking slow modulations in synaptic gain using dynamic causal modelling: Validation in epilepsy[END_REF][START_REF] Cooray | Dynamic causal modelling of electrographic seizure activity using Bayesian belief updating[END_REF]. It is also worth mentioning that we infer the dynamics of system with coupled fast and slow time-scales (cf., Eq. ( 3)), therefore, the variations in slow variable depend on the hidden states of fast activity while it is assumed that only the activity of fast variable is observed. In this study, the time-scale separation in Epileptor model enabled us to capture reliably full evolutions of complex dynamics, ranging from pre-ictal to onset, ictal evolution and offset [START_REF] Jirsa | On the nature of seizure dynamics[END_REF], rather using time-varying parameters [START_REF] Lopez-Cuevas | State and parameter estimation of a neural mass model from electrophysiological signals during the status epilepticus[END_REF]. Future extensions to the current work could examine explicitly the non-stationary dynamics of networks in order to investigate conditions for mechanism of seizure initiation whether the seizure onset is more likely to occur through a deterministic parameter changes as in a bifurcation [START_REF] Breakspear | A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis[END_REF][START_REF] Jirsa | On the nature of seizure dynamics[END_REF] or it is a jump phenomenon due to the noise-driven transition between bistable attractors (Kramer et al., 2012;[START_REF] Jedynak | Temporally correlated fluctuations drive epileptiform dynamics[END_REF][START_REF] Karoly | Seizure pathways: A model-based investigation[END_REF] In conclusion, we established a link between the probabilistic modeling and personalized brain network modeling in order to systematically predict the location of seizure initiation in a virtual epileptic patient. We demonstrated step by step, how the proposed framework allows one to infer the spatial map of epileptogenicity based on large-scale brain network models that are derived from noninvasive structural data of individual patients. The approach rests on advanced efficient sampling algorithms that provide accurate and reliable estimates validated by the posterior behavior analysis and convergence diagnostics. In summary, with the help of PPLs, the use of personalized brain network models offers a proper guidance for development of comprehensive clinical hypothesis testing and novel surgical intervention.

  have proposed a novel approach namely Virtual Epileptic Patient (VEP) to brain interventions based on personalized brain network models derived from non-invasive structural data of individual patients. The VEP model is a large-scale computational model of an individual brain that incorporates personal data such as the locations of seizure initiation, subject-specific brain connectivity, and MRI lesions to inform patient-specific clinical monitoring and improve surgical outcomes. It has been previously shown that the VEP model is able to realistically mimic the evolution of epileptic seizures in a patient with bitemporal epilepsy

  specific diagnostics can be used to monitor the convergence of samples; the number of divergent leapfrog transitions (due to highly varying posterior curvature), the step size used by NUTS in its Hamiltonian simulation (if the step size is too small, the sampler becomes inefficient, whereas if the step size is too large, the Hamiltonian simulation diverges), and the depth of tree used by NUTS, which is related to the number of leapfrog steps taken during the Hamiltonian simulation[START_REF] Gelman | Bayesian Data Analysis, Third Edition. Chapman & Hall/CRC Texts in Statistical Science[END_REF]; Stan Development Team (2018).2.8. Evaluate posterior fitUsing synthetic data for fitting allows us to validate the inference as we know the ground-truth of the parameters being inferred. Therefore, we can use standard error metrics to measure the similarity between the inferred parameters and those used for data generation. The metrics we used to validate our inference are confusion matrix, posterior shrinkage, and posterior z-score.

Finally, to invert

  the BVEP for the simulated data sets, we used two popular open-source PPL tools for flexible probabilistic inference: Stan (Stan Development Team, 2018), and PyMC3[START_REF] Salvatier | Probabilistic programming in Python using PyMC3[END_REF]. Stan language can be run in different interfaces, whereas PyMC3 provides several MCMC algorithms for model specification directly in native Python code. By specifying the model density functions in these tools, the gradients of functions are computed through automatic differentiation[START_REF] Baydin | Automatic differentiation in machine learning: A survey[END_REF][START_REF] Margossian | A review of automatic differentiation and its efficient implementation[END_REF], a powerful technique for algorithmic computation of derivatives, to efficiently approximate the log-posterior density by NUTS and ADVI. The computation of independent MCMC chains can also be performed in parallel on separate processors. In this study, we used Stan command line interface, whereas all the codes for simulations and posterior-based analysis were implemented in Python. The model simulation and parameter estimation were performed on a Linux machine with 3.0 GHz Intel Xeon processor and 32 GB of memory.

Fig. 1 :

 1 Fig.1: Schematic illustration of workflow in the BVEP brain model. The approach to build the BVEP comprises two main steps: constructing the VEP model, and then embedding VEP in a PPL tool to infer and validate the model parameters. To build the VEP model, we take the following steps: First, the patient undergoes non-invasive brain imaging (MRI, DTI). Based on these images, the brain network anatomy including brain parcellation and the patient's connectome are provided from the reconstruction pipeline. Then, a neural population model is selected for each brain region to define the network model. In VEP, the Epileptor model is defined on each network node that are connected through structural connectivity derived from diffusion tractography. Put together, TVB simulations allow to mimic the empirical neuroimaging signals. Then, model fitting is performed using NUTS/ADVI algorithms within a PPL tool (in this note, the brain source activity as the observation, and the VEP model as the generative model translated in Stan/PyMC3). Finally, cross validation can be performed by WAIC/LOO from the existing samples to assess the model's ability in new data prediction, thus, in order to refine the network pathology.

Fig. 2 .Fig. 2 :

 22 Fig. 2. Parcellation of the reconstructed brain and the patient's brain network are shown in Fig. 2A and B, respectively. Following Desikan-Killiany parcellation used in the reconstruction pipeline, the patient's brain is divided into 68 cortical regions and 16 subcortical structures. Fig. 2C illustrates the structural connectivity matrix derived from diffusion tractography of the patient. Following the virtualization of the patient's brain, we used TVB to simulate the reconstructed VEP brain network model. The simulated time series of fast activity variable in full VEP brain model are illustrated in Fig. 2D. The different brain node types namely HZ, PZ, and EZ are encoded in green, yellow and red, respectively. When the Epileptors are in isolation(i.e., K = 0; no network coupling), the seizures are triggered only in the regions defined as EZ , whereas no seizure propagation can be observed in other regions (see Fig.S2A and D). However, by coupling the Epileptors through the structural connectivity matrix of the patient (see Fig.2C), the spatial recruitment pattern can be observed in the candidate brain regions defined as PZ (see Fig.2D). In contrast to patient 2 (see Fig.S2F) where only one of the PZ is recruited, here, due to the strong coupling connections to regions specified as PZ as well as the high excitability value of these nodes, the seizure propagates to all other candidate brain regions specified as PZ (nodes number 6, 12, and 28). Average of fast activity variable inferred by inverting reduced VEP model is illustrated by the dashed line in Fig.2D. It can be seen that there is a remarkable similarity between the simulated and the predicted seizures regarding the seizure initiation, propagation and termination. Note that the simulation illustrates the activity of fast variable in full VEP brain model (i.e., x 1,i (t) in Eq. (2)), whereas the inferred envelope of time series demonstrates the trajectories from inversion of reduced VEP model (cf. Eq. (3)). The estimated densities of the excitability parameters η i for different brain node types are shown in Fig.2E. From this figure, we observe that the true value of excitability parameter (dashed vertical line) is under the support of the estimated posterior density across different brain regions.

Fig. 3 :

 3 Fig.3: Accuracy of the estimated spatial map of epileptogenicity across different brain regions using NUTS algorithm for patient 1. (A) Exemplary of observed data (dash-dotted lines) versus the prediction for three brain node types defined as HZ (green), PZ (yellow), and EZ (red). The shaded area depicts the ranges between the 5th and 95th percentiles of the posterior predictive distribution. (B) Violin plots of the estimated densities of the η i for 84 brain regions. The true values are displayed by the filled black circles. (C) The distribution of posterior z-scores versus posterior shrinkages implies an ideal Bayesian inversion. (D) The confusion matrix of the estimated spatial map of epileptogenicity. The pre-defined class for all the brain nodes labeled as HZ, PZ, and EZ are accurately predicted (accuracy=1.0, misclass=0.0).

Fig. 5 :

 5 Fig. 5: The estimated spatial map of epileptogenicity by NUTS algorithm in comparison to ADVI. Exemplary histogram and the kernel density estimates of the samples obtained by NUTS are illustrated in panel (A) versus the approximation by mean-field variant of ADVI shown in panel (B). For all the brain nodes included in the analysis, the prior (shown in blue) was assumed as N (-2.5, 1.0). The dashed vertical lines indicate the true values.

Fig. 6 :

 6 Fig. 6: NUTS and ADVI convergence diagnostics. (A) The samples generated by NUTS from joint posterior probability distribution between the hyperparameter pair (σ, σ ). Here, the non-centered form of parameterization yields independent samples from the posterior distribution. (B) The centered form of sampling leads to high correlation between hyper-parameters indicating that the sampler was not efficiently exploring the posterior distribution. (C) The samples from approximate joint posterior probability distribution using the mean-field variant of ADVI. (D) Shows the R values for the non-centered form of sampling, which are lower than 1.05 for all estimated hidden states and parameters implying that the MCMC has converged. (E) The high values of R returned by the centered form of sampling indicate that the chain has not converged. (F) The variational objective function (ELBO) versus the number of ADVI iterations.

  Fig. 1). The workflow to build the BVEP brain model consists of two main steps: in the first step, we constructed the VEP i.e., the personalized large-scale brain network model of epilepsy spread. In the VEP model, the dynamics of brain nodes are governed by the neural population model of epilepsy namely Epileptor, which is a generic model to realistically reproduce the onset, progression and offset of seizure patterns across species and brain regions (Jirsa et al., 2014). The Epileptors are coupled through the patient's connectome to combine the mean-field model of abnormal neuronal activity with the subject-specific brain's anatomical information derived from non-invasive diffusion neuroimaging techniques (MRI, DTI). Together with patient's data, the VEP model was then furnished with the spatial map of epileptogenicity across different brain regions. In the second step, we embedded the VEP as the generative model in PPL tools (Stan/PyMC3) to infer and validate the spatial map of epileptogenicity across different brain regions. Using the PPLs along with high-performance computing to run several MCMC chains in parallel enables systematic and efficient parameter inference to fit and validate the BVEP model against the patient's data.

  3, and S5) is a sufficient feature for correctly estimating the spatial map of epileptogenicity. Our results indicated that the BVEP model is able to accurately estimate the spatial map of epileptogenicity across different brain regions(cf. Figs. 3, and S5). The true value of excitability for all brain nodes included in analysis was under the support of estimated posterior densities with 100% classification accuracy based on the confusion matrix. In addition, the concentration of distribution toward small z-scores along with concentration towards large posterior shrinkages (i.e., the bottom right corner in Fig.3C) confirmed the reliability of model inversion. Note that relying on the accuracy obtained by confusion matrix may be inconclusive, since the accuracy of estimation return by this metric depends only on the mean value of estimated posterior densities. For instance, consider inferences where the posterior modes are nearly identical to the ground-truth, however there is a large uncertainty over the estimations.

  in-silico brain network models based on Bayesian inference within PPL tools such as Stan and PyMC3. Although several PPL libraries have been developed for Bayesian inference, only a few of them are built around efficient sampling algorithms such as NUTS that avoids the random walk behavior and and sensitivity to correlated parameters. Both Stan and PyMC3 provide NUTS and ADVI with automatic differentiation to efficiently compute gradients without requiring user intervention. Stan is a generic and flexible software package that has interfaces for common data science languages, also providing extensive diagnostics for MCMC convergence. PyMC3 provides several MCMC algorithms by model specification directly in native Python code. Our implementation in both Stan and PyMC3 result in similar estimation of spatial map of epileptogenicity across brain regions (cf. Figs. S5 and S6) indicating that BVEP is a platform-independent approach. However, we required a larger number of warm-up iterations in PyMC3 to arrive at the same posterior convergence achieved by our implementation in Stan. This is due to the differences in implementation of NUTS in Stan and PyMC3.

  . The Bayesian inversion in the current work is based on the auto-tuning algorithms such as NUTS and ADVI accomplished with fast automatic differentiation for the calculation of gradients. This allows us to efficiently sample from complex and high-dimensional posterior distributions with correlated parameters compared to the traditional sampling algorithms. The inferences in the presented framework is also enriched with several MCMC convergence diagnostics to assess the reliability of the estimations. Various noninvasive and invasive methods have been used to improve presurgical evaluation in identification of the EZ, and consequently to increases surgery success rates. Employing the BVEP model in clinical therapies and brain interventions will require quantification of the model outcomes in fitting empirical secondary functional signals of patients such as EEG, MEG, SEEG, and fMRI signals. In this framework, it is straightforward to incorporate further knowledge such as MRI lesions and clinical hypothesis on EZ from pre-surgical evaluation. Since the BVEP model can be considered as a generic approach towards large-scale brain modeling, it offers promising avenue for inference from clinically used non-invasive imaging signals (EEG, MEG, fMRI), and invasive measurements such as SEEG signals. Our pre-liminary results indicate that the proposed approach in this study is able to successfully fit against the patient's empirical SEEG data (not shown). Note that in the case of empirical SEEG recording, the source localization is an ill-posed problem due to the sparsity of lead-field matrix, which can affect the accuracy of the estimates. In principle, it is possible that the surgical strategies can be systematically tested using the BVEP model, however, the real clinical application remains to be investigated and validated in future work.

Information Sharing Statement

The main source codes needed to reproduce the presented results are available on GitHub (https://github.com/ins-amu/BVEP).

Author contributions

M.H., M.M.W., and V.K.J. designed the study. M.G. acquired the data.

M.H. performed the study. M.H, A.N.K., V.S., M.G., F.B, M.M.W., and V.K.J. wrote the manuscript.

Acknowledgements

The authors wish to acknowledge the financial support of the following

Competing interests

The authors declare no competing interests.