A computer-implemented method for inferring epileptogenicity of brain regions

Abstract

The present invention relates to a method for inferring the epileptogenicity of brain regions of an epileptic patient's brain that are not observed as eliciting or not eliciting seizure activity, the method comprising the steps of: providing a computerized model that models various regions of a primate brain and the connectivity between said regions; providing said computerized model with a model that is capable of reproducing the dynamics of epileptic seizures in the primate brain; providing structural data of the epileptic patient's brain and personalizing said computerized model using said structural data to obtain a virtual epileptic patient (VEP) brain model; converting a state-space representation of the virtual epileptic patient (VEP) brain model into a probabilistic programming language (PPL) using probabilistic state transitions to obtain a probabilistic virtual epileptic patient brain model (BVEP); and obtaining electroencephalogram or magnetoencephalogram data of the patient's brain and fitting the probabilistic virtual epileptic patient brain model to said data to infer the epileptogenicity of the unobserved brain regions.

Description

Technical Field

[0001] The present invention relates to a probabilistic method for inferring the epileptogenicity of brain regions that are not observed or not observed as having seizure activity in the brains of epileptic patients.

Background Art

[0002] Inverting a model, i.e., finding a set of model parameters that provides the best possible fit to the observed data, is a difficult task in statistical inference. The Bayesian framework provides a powerful and principle-based method for inferring parameters from experimental data and predicting models for a variety of applications. Within the context of neuroimaging, Bayesian methods are widely used to infer the intrinsic parameters of neuronal populations and / or the interactions between neuronal populations in a pre-specified neural network from neurophysiological data. In addition, Japanese Patent Application Publication No. 2019-527105 discloses a method for identifying an epileptogenic zone (EZ) in the brain of an epilepsy patient and modulating the epileptogenicity of the patient's brain in order to achieve an appropriate depiction of the EZ, which is essential for the success of an intervention as a resection operation.

[0003] Gradient-free sampling algorithms such as metropolis-hastings, gibbs sampling, and slice sampling are generally known to be unable to efficiently explore the parameter space when applied to large-scale inverse problems, as often faced in the use of whole-brain imaging for clinical diagnosis. In particular, traditional Markov chain Monte Carlo (MCMC) does not mix well in high-dimensional parameter spaces containing correlated variables. In contrast, gradient-based algorithms such as Hamiltonian Monte Carlo (HMC) are computationally expensive but far superior to gradient-free sampling algorithms in terms of the number of independent samples generated per unit of computational time. This class of sampling algorithms provides efficient convergence and exploration of the parameter space, even in very high-dimensional spaces that can exhibit strong correlations. Nevertheless, the efficiency of gradient-based sampling methods such as HMC is highly sensitive to algorithm parameters specified by the user. More advanced MCMC sampling algorithms, such as the No-U-Turn sampler (NUTS), a self-tuning variant of HMC, address these issues by adaptively adjusting algorithm parameters. These algorithms have been shown to efficiently sample from high-dimensional target distributions and can solve complex inverse problems conditioned on large datasets as observations.

[0004] MCMC has the advantage of being non-parametric and asymptotically exact in the limit of long / infinite runs. Among other alternatives, variational inference (VI) turns Bayesian inference into an optimization problem, which is usually a much faster computation than the MCMC method. However, traditional derivations of VI require major model-specific tasks such as defining a variational family suitable for the probabilistic model, computing the corresponding objective function, computing the gradient, and running a gradient-based optimization algorithm. Automatic Differentiation Variational Inference (ADVI) automatically solves these problems.

[0005] Probabilistic programming languages (PPLs) each provide an efficient implementation of automatic Bayesian inference in user-defined probabilistic models by featuring next-generation MCMC sampling and VI algorithms such as NUTS and ADVI. With the help of PPLs, these algorithms utilize automatic differentiation for computing derivatives in computer programs to avoid random walk behavior and sensitivity to correlation parameters. In particular, Stan and PyMC3 are high-level statistical modeling tools for Bayesian inference and probabilistic machine learning, providing advanced inference algorithms such as NUTS and ADVI enhanced with extensive and reliable diagnostic capabilities. Although PPLs enable automatic inference, the performance of these algorithms can be sensitive to the form of parameterization. The appropriate form of reparameterization in probabilistic models to improve the inference efficiency of system dynamics (governed by a series of non-linear probabilistic differential equations) remains a difficult problem.

[0006] On the other hand, personalized large-scale brain network modeling has recently received attention because it has the potential to improve medical treatment strategies. In individualized whole-brain modeling approaches, patient-specific information such as anatomical connectivity obtained from non-invasive imaging techniques is combined with mean-field models of local neuronal activity to simulate individual spatio-temporal brain activity at the macroscopic scale. The Virtual Brain (TVB) is an open-access computational framework written in Python for reconstructing and evaluating personalized architectures of the brain by using individual subject data. This neuroinformatics platform integrates computational modeling of the brain and multimodal neuroimaging data to systematically simulate individual spatio-temporal brain activity. However, currently, there is no specific workflow for automatic model inversion and data fitting verification in the preparation of TVB.

[0007] More recently, a novel approach for brain intervention based on personalized brain network models derived from non-invasive structural data of individual patients, namely, Virtual Epilepsy Patient (VEP), has been proposed. The VEP model is a large-scale computational model of an individual brain that incorporates personal data such as the onset location of seizures, subject-specific brain connectivity, and MRI lesions to inform patient-specific clinical monitoring and improve surgical outcomes. It has previously been shown that the VEP model can realistically simulate the progression of epileptic seizures in bilateral temporal lobe epilepsy patients. However, the inverse problem of such large-scale brain network models is a difficult task due to the inherent non-linear dynamics of each brain network node, as well as the large number of model parameters involved and the observations commonly encountered in brain imaging settings. SUMMARY OF THE INVENTION

[0008] Therefore, in order to systematically predict the onset location of seizures in virtual epilepsy patients, it is necessary to establish a useful link between the most popular probabilistic programming tools (e.g., Stan / PyMC3) and personalized brain network modeling (e.g., VEP model). The present invention enables the construction of a Bayesian Virtual Epilepsy Patient (BVEP), in particular, as a probabilistic framework designed to infer the hidden / non-observed dynamics of a personalized large-scale brain model of the spread of epilepsy generated by TVB.

[0009] According to a first aspect, the present invention is a method for inferring the epileptogenicity of brain regions that are not observed as having had seizure activity in the brain of an epilepsy patient or not observed as not having had seizure activity, comprising: providing a computerized model that models various regions of the primate brain and the connectivity between said regions; providing to the computerized model a model capable of reproducing the dynamics of epileptic seizures in the primate brain, said model being a function of parameters of the epileptogenicity of a region of said brain; To obtain a virtual epilepsy patient (VEP) brain model, providing brain structure data of an epilepsy patient and personalizing the computerized model using the structure data; To obtain a probabilistic virtual epilepsy patient brain model (BVEP), converting a state space representation of a virtual epilepsy patient (VEP) brain model into a probabilistic programming language (PPL) using probabilistic state transitions; To infer the epileptogenicity of a brain region that is not observed or not observed as having elicited seizure activity in a patient's brain, obtaining electroencephalogram data or magnetoencephalogram data of the patient's brain and fitting a probabilistic virtual epilepsy patient brain model to the data; Relates to a method comprising.

[0010] Preferably, - the probabilistic programming language is a Bayesian programming language, the probabilistic virtual epilepsy patient brain model is a Bayesian virtual epilepsy patient (BVEP) brain model, and the epileptogenicity of brain regions that are not observed as having occurred or not observed as not having occurred is inferred using Bayesian inference. - The structural data of the brain of an epilepsy patient includes non-invasive T1-weighted imaging data and / or diffusion MRI image data. - A model that can reproduce the dynamics of epileptic seizures in the primate brain is a model that reproduces the dynamics of the start, progression, and end of seizure events. It includes state variables that couple two oscillatory dynamical systems with three different time scales: the fastest time scale at which the state variable describes high-frequency discharges during the seizure state of the seizure period, the intermediate time scale at which the state variable represents slow spike-wave oscillations, and the slowest time scale at which the state variable is involved in the transition between the interictal state and the seizure state. The degree of epileptogenicity of a brain region is represented through the value of the excitability parameter. - To obtain a probabilistic virtual epilepsy patient brain model, a spatial map of the epileptogenicity of the patient's brain is provided. The spatial map of epileptogenicity classifies the brain regions of the patient's brain into epileptogenic zones (EZs) that may autonomously trigger epileptic seizures, propagation zones (PZs) that do not autonomously trigger seizures but may appear during the progression of seizures, and normal zones (HZs) that do not autonomously trigger seizures. - The probabilistic virtual epilepsy patient brain model is generated according to a generative model based on the state-space representation of the virtual epilepsy patient. - The state-space representation of the virtual epilepsy patient is in the form of

Number

Number

Number

Number

Number

Number

Number

Number

Number

[0011] Other features and aspects of the present invention will be apparent from the following description and the accompanying drawings.

Brief Description of the Drawings

[0012]

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Number

Number

DETAILED DESCRIPTION OF THE INVENTION

[0013] The present invention relates to a method for inferring the epileptogenicity of brain regions that are not observed or not observed as having elicited seizure activity in the brains of epileptic patients. This is a computerized probabilistic method for inferring a spatial map of epileptogenicity across various brain regions of a personalized epileptic brain patient, where seizures may start in a virtual region and propagate to candidate brain regions.

[0014] The method according to the present invention includes various steps implemented by a computer. A computer-readable medium is encoded with computer-readable instructions for performing the steps of the method according to the present invention.

[0015] This includes providing a computerized model that models various regions of the primate brain and the connectivity between said regions. This brain is The Virtual Brain. This is a neuroinformatics platform for whole-brain network simulations using biologically realistic connectivity. This simulation environment enables model-based inference of neurophysiological mechanisms across various brain scales underlying the generation of macroscopic neuroimaging signals, including functional magnetic resonance imaging (fMRI), electroencephalography (EEG), and magnetoencephalography (MEG). This allows for the reproduction and evaluation of personalized configurations of the brain by using individual subject data.

[0016] This further includes providing to the computerized model a model capable of reproducing the dynamics of epileptic seizures in the primate brain, which is a function of the parameters of epileptogenicity of brain regions.

[0017] Preferably, a model capable of reproducing the dynamics of epileptic seizures in the primate brain is a model that reproduces the dynamics of the onset, progression, and termination of seizure events, and includes state variables that couple two oscillatory dynamical systems on three different time scales: the fastest time scale on which the state variables describe the high-speed discharges during the seizure state of the seizure period, the intermediate time scale on which the state variables represent slow spike-wave oscillations, and the slowest time scale on which the state variables are involved in the transition between the interictal state and the seizure state. The degree of epileptogenicity of a brain region is represented through the value of the excitatory parameter.

[0018] In addition, the method according to the present invention includes the step of providing structural data of the brain of an epileptic patient and personalizing the computerized model using the structural data in order to obtain a virtual epileptic patient (VEP) brain model. The structural data is, for example, image data of the patient's brain obtained using magnetic resonance imaging (MRI), diffusion-weighted magnetic resonance imaging (DW-MRI), nuclear magnetic resonance imaging (NMRI), or magnetic resonance tomography (MRT). Preferably, the structural data of the brain of an epileptic patient includes non-invasive T1-weighted imaging data and / or diffusion MRI image data.

[0019] The method according to the present invention further includes the step of converting the state space representation of a virtual epileptic patient (VEP) brain model into a probabilistic programming language (PPL) using probabilistic state transitions in order to obtain a probabilistic virtual epileptic patient brain model (BVEP).

[0020] Preferably, the probabilistic programming language is a Bayesian programming language, the probabilistic virtual epileptic patient brain model is a Bayesian virtual epileptic patient (BVEP) brain model, and Bayesian inference is used to infer the epileptogenicity of brain regions that are not observed whether they have occurred or not.

[0021] Preferably, to obtain a probabilistic virtual epilepsy patient brain model, an epilepsyogenic spatial map of the patient's brain is provided, and the epilepsyogenic spatial map classifies the brain regions of the patient's brain into an epilepsyogenic zone (EZ) that may autonomously trigger an epileptic seizure, a propagation zone (PZ) that does not autonomously trigger a seizure but may appear during the progression of a seizure, and a normal zone (HZ) that does not autonomously trigger a seizure.

[0022] Preferably, the probabilistic virtual epilepsy patient brain model is generated according to a generative model based on the state space representation of the virtual epilepsy patient.

[0023] More preferably, the state space representation of the virtual epilepsy patient is

Number

Number

Number

Number

[0024] Preferably, to obtain a probabilistic virtual epilepsy patient (BVEP) model, the state space representation of the virtual epilepsy patient (VEP) model is incorporated into the probabilistic virtual epilepsy patient (BVEP) model as a state transition probability.

[0025] More preferably, the state transition probability is as follows:

Number

Number

[0026] More preferably, the generative model is defined by terms of the likelihood and prior density of the model parameters, and their product is the joint density:

Number

Number

Number

[0027] The method according to the present invention further includes the steps of acquiring electroencephalogram data or magnetoencephalogram data of a patient's brain, fitting a probabilistic virtual epilepsy patient brain model to the data, and inferring the excitability of the brain regions that are not observed whether or not the seizure activity of the patient's brain has occurred.

[0028] More preferably, a sampling algorithm is implemented to infer the excitability of the brain regions that are not observed whether or not the seizure activity of the patient's brain has occurred.

[0029] Even more preferably, the sampling algorithm is a Markov chain Monte Carlo or variational inference algorithm.

[0030] Example 1: Materials and Methods In one example, the method according to the present invention is based on personalized brain network modeling and Bayesian inference, as schematically shown in FIG. 1. This method enables the construction of BVEP according to two main steps: constructing a VEP model and then embedding the VEP into a PPL tool to infer and verify model parameters. The following steps are performed to construct the VEP model: First, non-invasive brain imaging (MRI, DTI) is performed on the patient. Based on these images, the anatomical structure of the brain network, including brain region segmentation and the patient's connectome, is provided from the reconstruction pipeline. Next, a neural cell population model is selected for each brain region to define the network model. In VEP, the Epileptor™ model is defined at each network node connected through the structural connectivity derived from diffusion tractography. Such a model is disclosed, for example, in the published document “On the nature of seizure dynamics”, Jirsa et al., Brain 2014, 137, 2210-2230, which is incorporated herein by reference. This includes five state variables that function on three different time scales. On the fastest time scale, the state variables describe the high-speed discharges during seizures. On the slowest time scale, the dielectric state variable describes slow processes such as extracellular ion concentration, energy consumption, and tissue oxygenation. The system exhibits high-frequency oscillations during the seizure state through the state variables. The autonomous switching between the interictal state and the seizure state is realized through the dielectric state variable via saddle-node and homoclinic bifurcation mechanisms for the start and end of the seizure, respectively. The switching is accompanied by direct current (DC) shifts recorded in vitro and in vivo. On the intermediate time scale, the other state variables describe the spike-and-wave electrographic pattern observed during seizures and the interictal and pre-seizure spikes when excited by the fastest system through coupling. In summary, the TVB simulation enables the simulation of empirical neuroimaging signals.Next, model fitting is performed using, for example, the NUTS / ADVI algorithm within the PPL tool (in this example, brain source activity as the observed values and the VEP model as the generative model are converted in Stan / PyMC3). Finally, cross-validation is performed, for example, by WAIC / LOO from existing samples to evaluate the predictive ability of the model for new data, and thus the network pathology can be refined. In other words, the workflow for constructing BVEP consists of two main steps: the step of constructing a VEP, that is, a personalized brain network model of epilepsy spread, and then the step of embedding the VEP model into a Bayesian framework to infer and validate model parameters. Following the formulation of VEP in state-space representation, a probabilistic reparameterization of the system dynamics is demonstrated. The proposed probabilistic reparameterization in BVEP has been shown to be able to efficiently invert non-linear state-space equations to infer the system dynamics. This approach enables the accurate estimation of the epileptogenic spatial map in a personalized brain network model of epilepsy spread using PPL. The Virtual Brain is used for brain network simulation, and Stan and PyMC3 are used to invert the simulated whole-brain model.

[0031] Below, it is shown step by step how to construct the BVEP model for a specific patient in order to fit the constructed brain model to in-silico data and validate the inference. The accuracy and reliability of the estimation are verified by several convergence diagnostics and post-hoc behavioral analyses.

[0032] Individual patient data In this study, two patients were selected: a 23-year-old woman with drug-resistant occipital lobe epilepsy (Patient 1) and a 24-year-old woman with drug-resistant temporal and frontal lobe epilepsy (Patient 2). The patients underwent standard clinical evaluations, the details of which have been described in a previous study (Proix et al., Individual brain structure and modelling predict seizure propagation. Brain 140, 641-654, 2017). The evaluations included non-invasive T1-weighted imaging (MPRAGE sequence, repetition time = 1900 ms, echo time = 2.19 ms, 1.0×1.0×1.0 mm, 208 slices), and diffusion MRI images (DTI-MR sequence, 64-direction set of angular gradients, repetition time = 10.7 s, echo time = 95 ms, 2.0×2.0×2.0 mm, 70 slices, 1000 s mm -2 of b-weighting). The images were acquired using a Siemens Magnetom Verio (trademark) 3T MR scanner.

[0033] Anatomical structure of the network A structural connectome was constructed in a reconstruction pipeline using publicly available neuroimaging software. The current version of the pipeline has evolved from the version described above (Proix et al., 2017). First, the recon-all command from the Freesurfer (trademark) package version v6.0.0 was used to reconstruct and segment the anatomical structure of the brain from the T1-weighted images. Then, the T1-weighted images were aligned with the diffusion-weighted images using the linear registration tool flirt (trademark) from the FSL package version 6.0 with a correlation ratio cost function of 12 degrees of freedom.

[0034] Next, for tractography, the MRtrix™ package version 0.3.15 was used. Spherical deconvolution with the dwi2fod tool was used to estimate the orientation distribution of nerve fibers from DWI, along with the response function estimated by the dwi2response tool using the Tournier algorithm. Next, the tckgen tool using the probabilistic tractography algorithm iFOD2 was used to generate 15 million fiber tracts. Finally, the connectome matrix was constructed with the tck2connectome tool using the Desikan-Killiany parcellation generated by FreeSurfer in the previous step. The connectome was normalized so that the maximum value was equal to 1.

[0035] Network model Typically, to construct a personalized brain network model, brain regions are defined using a parcellation scheme, and a set of equations are used to model local brain activity. Adopting such a data-driven approach to incorporate subject-specific brain anatomical information, network edges are represented by the structural connectivity of the brain obtained from non-invasive imaging data of individual patients. In the VEP model, the dynamics of brain network nodes are governed by the Epileptor equation (Jirsa et al., On the nature of seizure dynamics. Brain 137, 2210-2230, 2014) coupled through a structural connectivity matrix derived from diffusion-weighted MRI (dMRI) technology (Jirsa et al., The virtual epileptic patient: Individualized whole-brain models of epilepsy spread, 2017).

[0036] The Epileptor is a dynamic model of seizure progression and can realistically reproduce the onset, progression, and onset dynamics of events such as seizures. The Epileptor contains five state variables that couple two oscillator systems on three different timescales. On the fastest timescale, the variables x1 and y1 describe the high-frequency discharges during the seizure state of the seizure period. On the intermediate timescale, the variables x2 and y2 represent the slow spike-wave oscillations. On the slowest timescale, the dielectric state variable z is involved in the transition between the interictal state and the seizure state. Additionally, interictal and pre-seizure spikes are generated via the term g(x1).

[0037] Following Jirsa et al. (2014), the dynamics of the full Epileptor model are described by:

Equation

Equation

[0038] Following Jirsa et al. (2017), the full VEP brain model equation (N coupled Epileptors) is as follows:

Equation

Equation

[0039] Under the assumption of timescale separation, Proix et al. (2014) showed that the influence of the second neuron population of the Epileptor (i.e., the variables x2 and y2) can be ignored by averaging the coupled Epileptor equations to obtain a 2D reduced VEP model of the form:

Equation

[0040] Depending on the value of the excitability parameter η, the 2D Epileptor exhibits different stable regimes. When η < η c , the trajectories in the phase plane are attracted to a 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 normal, which means that it does not cause seizures without external input. As the value of η increases, the z-nullcline moves down, and a saddle-node bifurcation occurs at η = η c . When η > η c , the system exhibits an unstable fixed point where seizures can occur (the Epileptor is said to be epileptogenic). In this example, a 2D reduction of the VEP model is used for Bayesian inference of the epileptogenic spatial map, reducing the computational cost associated with model parameter estimation. The 2D reduction of the Epileptor allows for faster inversion while enabling prediction of the envelope of high-frequency discharges during the seizure state of the seizure period (i.e., the onset, propagation, and termination of the seizure pattern) (Proix et al., 2014; Jirsa et al., 2017).

[0041] Epileptogenic spatial map In addition to the patient's connectome that structurally constrains individual variability, the dynamics of the brain network model can be further constrained by formulating hypotheses regarding functional brain network components to generate more specific patterns of each individual's brain activity. In the case of epilepsy, clinical hypotheses regarding the location of epileptogenic regions or lesions enable refining the network pathology to better predict the onset and propagation of seizures in individual patients.

[0042] In the BVEP brain model, each network node can trigger seizures depending on its connectivity and excitability values. The parameter η controls the excitability of the tissue, and thus its spatial distribution is the subject of parameter fitting. In this study, depending on the excitability values, various brain regions are classified into three main types: - Epileptogenic zone (EZ): η > η c In this case, Epileptor may autonomously trigger seizures (brain regions involved in the initiation and initial organization of epileptic activity). - Propagation zone (PZ): η c - Δη < η < η c In this case, Epileptor does not autonomously trigger seizures, but may appear during the transition of seizures because their equilibrium states are close to the critical value. - Healthy zone (HZ): η < η c - Δη, Epileptor does not autonomously trigger seizures.

[0043] Based on the above dynamic characteristics, the epileptogenic spatial map across various brain regions includes the excitatory value of the EZ (high excitatory value), the PZ (smaller excitatory value), and all other regions classified as non-epileptogenic (HZ). However, it should be noted that since propagation is also determined by various other factors including connectivity and brain state dependence, intermediate excitatory values do not guarantee that seizures will appear in this region as part of the propagation region. In the BVEP brain model, clinical hypotheses can be formulated as prior knowledge about the spatial distribution of excitatory parameters. In this study, assuming no clinical hypothesis regarding specific brain regions, the same prior distribution was assigned to the excitatory parameters of all brain regions included in the analysis.

[0044] Probabilistic model An important component in constructing a probabilistic brain network model within the Bayesian framework is the generative model. The generative model is a probabilistic description of the mechanism by which, given a set of observations, the observed data is generated through some hidden states and unknown parameters. Thus, here the generative model has a mathematical formulation derived from a dynamic model that describes the evolution of the state variables of the model, given the parameters over time. This specification is necessary to construct the likelihood function. The full generative model is then completed by specifying the prior beliefs regarding the possible values of the unknown parameters.

[0045] The BVEP brain model presented in this study is constructed in two main steps. First, the VEP model equation, which provides the basic form of the data generation process that describes how epileptic seizures are generated. Second, the formulation of hypotheses regarding the spatial map of epileptogenicity in the brain as prior knowledge. The latter component provides a model that uses hypotheses regarding the spatial distribution of excitatory parameters across various brain regions.

[0046] The generative model in BVEP is formulated based on a system of non-linear probabilistic differential equations of the following form (so-called state-space representation):

Equation

Number

Number

Number

[0047] Considering the 2D reduction of the VEP model (see Equation (3)),

Number

Number

[0048] The state - space representation (see Equation (4)) that defines the dynamics of the hidden state x(t) has a state - transition probability:

Number

Number

[0049] A careful choice of re - parameterization has been previously shown to increase the effective sample size and reduce divergence, especially in regions of extreme curvature. To avoid the biased estimates caused by pathological samples and thus the strong correlations between parameters in the centered - form parameterization, the advantages of location - scale transformation are utilized to invert the non - linear state - space equation. This can decouple the correlations of the parameters representing the state variables at consecutive time steps.

[0050] The non - centered - form re - parameterization of the above distribution is as follows:

Number

[0051] Inference / Prediction The generative model is characterized by the joint probability distribution of the model parameters and the observed values

Number

Number

Number

Number

Number

Number

Number

Number

[0052] Posterior density

Number

[0053] The prior density of the excitatory parameter for all brain regions included in the analysis was considered to be a normal distribution with mean -2.5 and standard deviation 1.0, i.e., N(-2.5, 1.0). Furthermore, for the initial conditions of the system and the global coupling parameter K, a prior density with little information was set as a normal distribution centered on the ground truth with a standard deviation of 1.0. The prior density of the hyperparameters was considered to be a common prior density with little information N(0, 1.0).

[0054] After fitting a Bayesian model, it is often necessary to measure the predictive accuracy of the inferred model. Information criteria and leave-one-out cross-validation (LOO) are two rigorous methods for evaluating the predictive ability of a model on new data. By adopting existing simulations derived from the log-likelihood evaluated with the posterior distribution of parameter values, the widely applicable information criterion (WAIC) and Pareto smoothed importance sampling (PSIS) LOO can efficiently estimate the predictive accuracy of the fitted Bayesian model within a negligible computational time relative to the cost of model fitting.

[0055] Inferential diagnosis After running an MCMC sampling algorithm, it is necessary to perform some statistical analyses to evaluate the convergence of the MCMC samples. One simple way to evaluate the performance of the MCMC algorithm based on the posterior samples is to visualize how well the chain is mixed (i.e., whether the MCMC sampler efficiently explores all modes within the parameter space). This can be monitored in various ways, including trace plots (the evolution of parameter estimates from repeated MCMC), pair plots (to identify collinearity between variables), and autocorrelation plots (to measure the degree of correlation between draws of the MCMC samples). A more quantitative method for evaluating the convergence of MCMC to the stationary distribution is the potential scale reduction factor and the effective sample size N based on samples of the posterior model probability. eff to estimate. 's diagnosis provides an estimate of how much the variance can be reduced by running the chain longer. Each MCMC estimate is associated with its ​​​ has a statistical value that is essentially the ratio of the between-chain variance to the within-chain variance. [Number] If it is approximately less than 1.1, MCMC convergence has been achieved (approaching 1.0 in the case of infinite samples); otherwise, the chain needs to be run longer. Additionally, the eff statistical value of N indicates the number of independent samples represented by the chain. The larger the effective sample size, the higher the accuracy of the MCMC estimation. Note that these are necessary but not sufficient conditions for the convergence of MCMC samples.

[0056] In addition to the general MCMC diagnostics described above, NUTS-specific diagnostics can be used to monitor the convergence of the samples, the number of divergent leapfrog transitions (due to large changes in the curvature of the posterior distribution), the step size used by NUTS in Hamiltonian simulation (if the step size is too small, the sampler becomes inefficient; if the step size is too large, the Hamiltonian simulation diverges), and the depth of the tree used by NUTS, which is related to the number of leapfrog steps taken during Hamiltonian simulation.

[0057] Evaluation of Posterior Fit When using synthetic data for fitting, the ground truth of the inferred parameters is known, so the inference can be verified. Therefore, standard error metrics can be used to measure the similarity between the inferred parameters and the parameters used to generate the data. The metrics used to verify this inference are the confusion matrix, posterior shrinkage, and posterior z-score.

[0058] The confusion matrix is a metric for evaluating the accuracy of classification. Element q i,j is equal to the number of observations known to be in class i but predicted to be in class j, [Number] where Q is the total number of classes. In the BVEP model, three groups, namely HZ, PZ, and EZ, and thus Q = 3, were defined to classify brain regions.

[0059] Furthermore, to quantify the accuracy of the inference, the posterior z - scores (denoted by z) are plotted against the posterior shrinkage (denoted by s), which are defined as follows:

Number

Number

Number

Number

Number

Number

[0060] Synthetic dataset and model inversion To verify the inference using BVEP, the simulation function of The Virtual Brain (TVB) is utilized to generate a synthetic dataset. TVB is an open-source neuroinformatics tool written in Python for simulating large-scale brain network models based on individual subject data. This platform is widely used to simulate common neuroimaging signals including functional MRI (fMRI), EEG, SEEG, and MEG, and has a wide range of clinical applications from Alzheimer's disease, chronic stroke to human focal epilepsy (Jirsa et al., 2017).

[0061] In this study, TVB is used to reconstruct a personalized brain network model. To verify the inference of spatial epileptogenicity, the epileptic seizures of two patients are simulated. One is a simulation where the seizure spreads to all brain nodes designated as PZ (Patient 1), and the other is a simulation where the seizure spreads to some of the brain nodes designated as PZ (Patient 2). These datasets are generated using two different structural connectivity matrices and separate spatial maps of epileptogenicity.

[0062] The seizure activity of Patient 1 is simulated by setting two regions as EZ and three regions as PZ, where

Number

Number

[0063] To simulate the seizure activity of Patient 2, the nodes [Number] and [Number] Two brain regions were selected as EZ and five regions as PZ. For the regions selected as EZ, the excitability value was set to η ez = -1.5. The excitability value of PZ was set to η pz = -2.4, and all other regions were defined as HZ with η hz = -3.4.

[0064] In both synthetic datasets, to simulate the VEP model as a system of stochastic differential equations, the Euler - Maruyama integration scheme was used with an integration step of 0.04. Additive white Gaussian noise was introduced into 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 a mean of zero and a variance of (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.

[0065] Finally, to invert the BVEP for the simulated dataset, two popular open-source PPL tools for flexible probabilistic inference, namely, Stan and PyMC3, were used. The Stan language can be executed with various interfaces, while PyMC3 provides several MCMC algorithms for directly specifying models in native Python code. By specifying the model density function with these tools, the gradient of the function is calculated through automatic differentiation, a powerful technique for algorithmic calculation of derivatives, and the log posterior density is efficiently approximated by NUTS and ADVI. The calculation of independent MCMC chains can also be run in parallel on separate processors. In this example, the Stan command-line interface was used, while all the code for simulation and posterior-based analysis was implemented in Python. The simulation of the model and the estimation of the parameters were performed on a Linux machine with a 3.0 GHz Intel Xeon processor and 32 GB of memory.

[0066] Example 2: Results The results of the workflow in the BVEP model for estimating the epileptogenic spatial map across various brain regions for Patient 1 are shown in FIGS. 2A - 2E. The reconstructed brain parcellation and the patient's brain network are shown in FIGS. 2A and 2B, respectively. According to the Desikan - Killiany parcellation used in the reconstruction pipeline, the patient's brain is divided into 68 cortical regions and 16 subcortical structures. FIG. 2C shows the structural connectivity matrix derived from the patient's diffusion tractography. Following the visualization of the patient's brain, TVB was used to simulate the reconstructed VEP brain network model. The simulated time series of the fast - activity variable in the full VEP brain model is shown in FIG. 2D. Different brain node types, i.e., HZ, PZ, and EZ, are encoded in green, yellow, and red, respectively. When the Epileptor is isolated (i.e., K = 0, no network coupling), seizures are only induced in the region defined as EZ, while no seizure propagation can be observed in other regions (see FIGS. 2A and 2D). However, by coupling the Epileptor through the patient's structural connectivity matrix (see FIG. 2C), a spatial emergence pattern can be observed in the candidate brain regions defined as PZ (see FIG. 2D). In contrast to Patient 2, in whom only one of the PZs appears (see FIG. 2F), here, due to the strong coupling connections to the regions designated as PZ and the high excitability values of these nodes, seizures propagate to all other candidate brain regions designated as PZ (node numbers 6, 12, and 28). The average value of the fast - activity variable inferred by inverting the reduced VEP model is shown as a dashed line in FIG. 2D. It can be seen that there is a significant similarity between the simulated and predicted seizures with respect to the onset, propagation, and termination of seizures. The simulation exemplifies the fast - activity variable (i.e., x 1,i (t) in Equation (2)) in the full VEP brain model, while it should be noted that the inferred envelope of the time series shows the trajectory from the inversion of the reduced VEP model (see Equation (3)). The excitability parameter η for different brain node types iThe estimated density is shown in Fig. 2E. From this figure, it is observed that the true value of the excitability parameter (vertical dashed line) lies under the support of the estimated posterior density across various brain regions.

[0067] The accuracy of the epileptogenicity estimated spatial maps across various brain regions for Patient 1 by the BVEP implementation in Stan is presented in Figs. 3A - 3D. Similar results were obtained from the BVEP implementation in PyMC3. Fig. 3A is a comparison of the observed source activity and the inferred source activity for three brain node types designated as HZ, PZ, and EZ (node numbers 1, 6, and 7, respectively). The simulated data consists of 120 s of fast activity variables (i.e., x 1,i (t) in Eq. (2)) in a full VEP brain model sampled at 1000 Hz, which is downsampled by a factor of 10 to reduce the computational cost of the Bayesian inversion. The observed data is shown by the dashed line, while the shaded area indicates the range from the 5th percentile to the 95th percentile of the posterior predictive distribution. The activities of the selected brain nodes at HZ, PZ, and EZ are shown in green, yellow, and red, respectively. It is observed that the predicted time series based on samples from the posterior predictive distribution is in very good agreement with the simulation. Fig. 3B shows a violin plot of the estimated density of the excitability parameter for all 84 brain regions included in the analysis. The filled black circles represent the true parameter values used to generate the simulated data. It can be seen that the ground truth of the excitability parameter for all brain regions lies under the support of the estimated posterior density. As shown in Fig. 3C, the distribution of the posterior z - scores and posterior shrinkage for all inferred excitabilities demonstrates the reliability of the model inversion. Note that the concentration towards large shrinkage indicates that all posterior distributions in the inversion are well - discriminated, while the concentration towards small z - scores indicates that the true values are accurately included in the posterior distributions. Thus, the distribution in the lower right of the plot represents the ideal Bayesian inversion. To further confirm the accuracy of the estimated values of spatial excitability,

Number

[0068] To investigate whether BVEP is a framework independent of the platform, PyMC3 was also used to estimate the epileptogenic spatial map across various brain regions. For both patients analyzed, the same accuracy was obtained by inverting Equation 4 in Stan and PyMC3. These results indicate that the inversion of BVEP in Stan and PyMC3 leads to similar estimates of the epileptogenic spatial map across brain regions in both patients analyzed.

[0069] Furthermore, NUTS-specific diagnostics were monitored to check whether the Markov chain had converged. The diagnostic plots indicate that there are no divergent transitions in HMC, which indicates that the posterior density was efficiently explored. Also, none of the NUTS iterations reached the maximum tree depth (here the execution value of that NUTS was specified as 10.0), which indicates that the optimal number of leapfrog steps required for Hamiltonian simulation is well below the maximum value. In summary, these diagnostics verify that the samples by NUTS have converged to the target distribution.

[0070] To illustrate the mechanisms underlying the onset and propagation of seizures within the BVEP model, the phase-plane topologies of simulations (upper panel) and predictions (lower panel) characterizing the dynamics of various brain node types in the BVEP model are presented in FIG. 4. In the plotted phase plane, the x and z nullclines are colored dark gray, and the intersection of the nullclines identifies the fixed points of the system. From left to right, the columns correspond to brain nodes designated as HZ, PZ, and EZ, respectively. The black and white circles indicate stable and unstable fixed points, respectively. From FIGS. 4A and 4D, it can be observed that the trajectory of HZ (node number 1) is attracted to a stable fixed point of the system (on the left branch of the cubic x nullcline), which means that it does not trigger epileptic seizures. For PZ (node number 6), due to the coupling strength and an excitability value close to the critical value of epileptogenicity, the z nullcline moves downward and a bifurcation occurs, whereby seizures can propagate here (see FIGS. 4B and 4E). For EZ (node number 7), the system exhibits an unstable fixed point due to a high excitability value. In this regime, Epileptor has a limit cycle and autonomously triggers seizures (see FIGS. 4C and 4F). The topologies of the simulated and predicted phase-plane trajectories show a very good match except for the amplitude of the state variable z i and it should be noted that its estimation shows the results of a larger parameter recovery. Note that only the fast activity variable x 1,i is the subject of fitting as observational data.

[0071] To compare the inversion of BVEP by the NUTS scheme and the ADVI scheme, Figure 5 shows the histograms of the MCMC samples and the kernel density estimates of the posterior density generated from NUTS (left panel), and the histogram obtained by ADVI (right panel). From this figure, it is observed that NUTS and ADVI perform the estimation of the posterior density similarly, except that the mean-field ADVI slightly underestimates the variance compared to the estimation by the NUTS algorithm. However, when both methods are adopted, the true value of excitability (vertical dashed line) is below the support of the posterior density, which indicates that the parameter recovery was successful. The samples corresponding to the brain nodes designated as HZ, PZ, and EZ are shown in gray, light gray, and dark gray, respectively. Note that the prior density for all 84 brain regions included in the analysis was considered as a normal distribution centered at -2.5 with a standard deviation of 1.0 (i.e., N(-2.5, 1.0)), shown in blue. To invert the BVEP model by the NUTS algorithm, 200 sampling repetitions and 200 warm-ups were used with an expected acceptance probability of 0.95, while for running ADVI, the maximum number of repetitions and the convergence tolerance were set to 50000 and 0.001, respectively. Regarding the computation time, with these algorithm configurations, sampling by NUTS took 23993.5 seconds, while the running time of ADVI was 5392.62 seconds.

[0072] Once the model parameters are estimated, it is necessary to evaluate the convergence of the MCMC samples. To verify the reliability of the inferred estimates, the potential scale reduction factor, being the most reliable quantitative metric for MCMC convergence

Number

number

number

number

[0073] Furthermore, a scatter plot of samples extracted from the joint posterior probability distribution between hyperparameters σ and σ’ estimated by mean-field ADVI is shown in FIG. 6C. By definition, the mean-field variant of ADVI ignores the mutual correlations between parameters, so the samples extracted using mean-field ADVI do not show the correlations between hyperparameters. Finally, to check the convergence of ADVI, the variational objective function, the evidence lower bound (ELBO), is plotted against the number of iterations (see FIG. 6F). Although the algorithm appears to converge after 10,000 iterations, the algorithm is run for several thousand more iterations to ensure convergence until the change in ELBO falls below the tolerance value of 0.001.

[0074] Finally, the present invention presents a probabilistic framework for inferring a spatial map of epileptogenicity, i.e., a Bayesian virtual epilepsy patient (BVEP), to develop a personalized large-scale brain model of the spread of epilepsy (see Figure 1). The workflow for constructing the BVEP brain model consists of two main steps: In the first step, a VEP, i.e., a personalized large-scale brain network model of the spread of epilepsy, is constructed. In the VEP model, the dynamics of brain nodes are governed by a neuronal population model of epilepsy, i.e., Epileptor, which is a general model for realistically reproducing the onset, progression, and termination of seizure patterns across species and brain regions (Jirsa et al., 2014). Epileptor is coupled through the patient's connectome to combine an average-field model of abnormal neuronal activity with subject-specific anatomical information of the brain derived from non-invasive diffusion neuroimaging techniques (MRI, DTI). Along with the patient's data, a spatial map of epileptogenicity across various brain regions was provided to the VEP model. In the second step, the VEP was embedded as a generative model into a PPL tool (Stan / PyMC3) to infer and validate a spatial map of epileptogenicity across various brain regions. Using PPL with high-performance computing, by running several MCMC chains in parallel, systematic and efficient parameter inference for fitting and validating the BVEP model to the patient's data becomes possible.

[0075] To demonstrate the potential function of BVEP in predicting the onset and propagation of seizures, various spatial maps of epileptogenicity were used to simulate both simple and complex seizure spreads (see Figure S2). These synthetic data were used for fitting. This is because, given the ground truth of the model parameters, standard error metrics such as confusion matrices, posterior shrinkage, and posterior z-scores can be used to verify the accuracy of the estimation, and thus the performance of the proposed method can be evaluated. As a result, in both synthetic datasets, it was demonstrated that by inverting the large-scale brain network model with the help of PPL (Stan / PyMC3), a significant similarity between the simulated seizure activity and the predicted seizure activity can be achieved regarding onset, propagation, and termination. The simulations were generated by a full VEP model that includes five state variables at each brain node, but even the 2D reduced variant of the model was able to successfully predict important data features such as the onset, propagation, and termination of seizure patterns while considerably reducing the computational time of Bayesian inference. This 2D reduction is limited to modeling the mean value of high-speed discharges during the seizure state of the seizure period, which is a feature sufficient to properly estimate the epileptogenic spatial map as shown (see Figures 3 and S5). As a result, the BVEP model was shown to be able to accurately estimate epileptogenic spatial maps across various brain regions (see Figures 3A - 3D). The true values of excitability for all brain nodes included in the analysis were under the support of the estimated posterior density with a classification accuracy of 100% based on the confusion matrix. Furthermore, the reliability of the model inversion was confirmed by the concentration of the distribution towards small z-scores and towards large posterior shrinkage (i.e., the lower right corner of Figure 3C). It should be noted that since the accuracy of the estimation returned by that metric depends only on the mean value of the estimated posterior density, it may not be decisive to depend on the accuracy obtained by the confusion matrix. For example, consider an inference where the mean value of the posterior density is approximately the same as the ground truth but there is a large uncertainty in the estimated values.In such cases, the confusion matrix may yield high-precision performance, but plotting the post hoc z-score against post hoc shrinkage is particularly useful for identifying misbehavior in inference, such as overfitting, or an inappropriately chosen prior density that biases the estimate.

[0076] Understanding the brain dynamics of epilepsy is important for developing therapeutic approaches for brain interventions to improve surgical outcomes. Using the theory of nonlinear dynamical systems, a complete classification of epileptic seizures has been widely investigated elsewhere (Jirsa et al., 2014), along with a sufficient description of the bifurcations that give rise to the onset, termination, and transition characteristics of seizures. In the description of the parameter space of the Epileptor, the beginning and end of seizures are explained by saddle-node bifurcations and homoclinic bifurcations. The new dynamical effects in the BVEP model strongly depend on the interaction between the network node model (Epileptor), the patient-specific structural connectivity (from dMRI), and the epileptogenic spatial maps (EZ, PZ, HZ). Due to the dynamical properties of the Epileptor model, brain regions were classified into three main types: EZ (exhibiting an unstable fixed point corresponding to the brain region involved in the onset of seizures), PZ (close to the saddle-node bifurcation corresponding to the candidate brain region involved in the propagation of seizures), and HZ (exhibiting a stable fixed point corresponding to normal brain regions). This approach makes it possible to define the epileptogenic spatial map based on the excitatory parameter values that are the target of fitting.

[0077] It is important to note that excitatory values close to the seizureogenic threshold do not guarantee that seizures originating from pathological brain regions (i.e., those involved in the onset of seizures designated as the EZ) will spread to the brain regions defined as the PZ. Detailed patient evaluations have reported that individual structural connectivity is essential for predicting the spatial spread of seizures. However, it has recently been shown that purely structural information alone is not sufficient to predict seizure spread and ultimate termination. Rather, abnormal activity in the emergence region is a complex network effect that depends on the interaction between multiple factors, including the epileptogenicity (node dynamics) of brain regions, individual structural connectivity (network structure) (Jirsa et al., 2017), and brain state dependence (network dynamics). Furthermore, due to the dynamics of coupled nonlinear systems (see Fig. S2), there are nonlinearities and multiple propagation patterns that can be observed with the same set of excitatory parameters. In this work, seizure emergence is characterized by the complex spatiotemporal dynamics of large-scale brain networks, i.e., seizures originate from local networks and emerge in candidate brain regions strongly coupled to pathological regions by disrupting their stable dynamics (in the case of K = 0, seizure emergence does not occur). Among the candidate brain regions for seizure spread, due to stronger connections to the pathological region defined as the EZ, node PZ idx={28} may appear due to weak global coupling. Rather, stronger coupling is required for the appearance of seizures in all other candidate brain regions. This is consistent with the experimental observation that seizures tend to have a common spatial onset in the same patients. Based on this knowledge, a prior density with low information content centered on the ground truth was set for the global coupling parameter. Overestimation of the global coupling parameter leads to misclassification of the PZ as HZ (see Figures 2E and 3B), while underestimation of the coupling may lead to misclassification of the PZ as EZ. However, stability analysis of the network dynamics indicates that seizure propagation can be controlled by an optimal intervention on the structural connectivity matrix, which means that the connectivity of the patient-specific network predicts the seizure propagation pattern. Therefore, it has been reported that seizure propagation may not be easily controlled by simple disconnection of individual nodes, such as in surgical treatment of epilepsy, and resection does not necessarily lead to liberation from postoperative brain seizures.

[0078] In this study, to gain a better understanding of the mechanisms underlying seizure initiation and propagation within the proposed approach, an analysis of the observed and predicted phase plane trajectories across various brain regions was performed (see Figure 4). For various brain node types (e.g., EZ, PZ, and HZ), the dynamics of seizure initiation and appearance in the phase plane were well captured by the prediction. From the perspective of inference, a good correspondence was observed between the observed phase diagram of the system, including equilibrium (intersection of nullclines), stability or instability of the equilibrium, and the flow of the trajectory. These results validate the Bayesian inversion procedure for understanding the spatiotemporal evolution of seizure activity and open the way for further research on possible seizure prevention methods.

[0079] To infer the seizure origin spatial map in a personalized whole-brain model of seizure spread, both the NUTS and ADVI schemes were used. Results from both inference schemes led to similar estimates of the seizure origin spatial map across brain regions, except that ADVI slightly underestimated the variance compared to the estimates by the NUTS algorithm (see Figure 5). The similarity between the inversions using the two schemes indicates that variational approximation provides a suitable alternative to NUTS sampling in the inversion of the BVEP model. This result demonstrates that when performing inference by ADVI, the computational cost is significantly reduced compared to NUTS (4 - 5 times faster in the algorithm configuration used), which can be important when applying the BVEP approach to large datasets of patient cohorts. ADVI is generally known to be computationally more attractive than NUTS, but it can be difficult to discover algorithmic problems with this approximation. The convergence of ADVI can be evaluated by monitoring the moving average of the ELBO change, while for NUTS, several general and specific diagnostics are provided to evaluate whether the Markov chain has converged. Additionally, ADVI can stack at local minima during gradient descent optimization, and its mean-field variant cannot cover all modes of the multimodal posterior density.

[0080] Finally, the efficiency of the transformed non-central parameterization was investigated. Consistent with previous studies showing that NUTS is sensitive to parameterization, the results here show that non-central parameterization for inverting the non-linear state space equation allows for an efficient exploration of the parameter space, while central-form sampling demonstrated an inefficient exploration due to high collinearity among model parameters (see Figure 6A and Figure 6D versus Figure 6B and Figure 6E). Furthermore,

Number

[0081] A novel approach for constructing personalized in-silico brain network models based on Bayesian inference within PPL tools such as Stan and PyMC3. Although several PPL libraries for Bayesian inference have been developed, only a few of them are built around efficient sampling algorithms such as NUTS that avoid random walk behavior and sensitivity to correlation parameters. Both Stan and PyMC3 provide NUTS and ADVI with automatic differentiation for efficiently calculating gradients without the need for user intervention. Stan is a general-purpose and flexible software package, has interfaces for common data science languages, and also provides extensive diagnostics for MCMC convergence. PyMC3 provides several MCMC algorithms by directly specifying models in native Python code. Implementations in both Stan and PyMC3 yield similar estimates of the spatial map of epileptogenicity across brain regions, indicating that BVEP is a platform-independent approach. However, more warm-up repetitions were required in PyMC3 to reach the same posterior convergence achieved with the Stan implementation. This is due to differences in the NUTS implementation between Stan and PyMC3. A comparison of implementations in Stan, PyMC3, and other alternative PPL packages is beyond the scope of this note.

[0082] The present invention is the first personalized large-scale brain network modeling method for inferring epileptogenic spatial maps (node characteristics) based on patient-specific whole-brain anatomical information (i.e., network structure derived from dMRI). Dynamic causal modeling (DCM) is a well-established framework for analyzing neuroimaging modalities (such as fMRI, MEG, and EEG) by neural mass models, and can infer the coupling (effective connectivity) between brain regions to infer how changes in neuronal activity in a brain region are caused by the activity of other regions through modulation of potential couplings. Using DCM, focal seizure activity in electrocorticogram (ECoG) data has recently been studied to estimate the main synaptic parameters or coupling connections using signals observed in human subjects. In another study, the Bayesian belief update scheme of DCM has been used to estimate the synaptic drivers of cortical dynamics during seizures from EEG / ECoG recordings at a slight computational cost. Although DCM can be used to model and track changes in excitability (inhibitory balance at the onset / end of seizures), these studies are based on a single neural mass model (i.e., a small number of cortical sources are modeled), and the non-linear ordinary differential equations representing the neural mass model are approximated by its linearization, whereby only the onset or end of seizures can be modeled, but not both. In this note, the Bayesian virtual epilepsy patient (BVEP) model can characterize the spatio-temporal non-linear dynamics of seizure propagation throughout the brain. This approach enables the description of the onset and end of the seizure state, as well as the alternation between the normal and seizure states. The BVEP approach depends on patient-specific structural data rather than formulating a pure inverse problem with respect to the unknown model parameters used in DCM. It is also worth mentioning that the dynamics of the system were inferred on coupled fast and slow time scales (see Equation (3)), and thus the changes in the slow variables depend on the hidden state of the fast activity while assuming that only the fast activity variables are observed.In this study, the timescale separation in the Epileptor model enabled capturing all transitions of complex dynamics from pre-seizure, onset, in-seizure progression, to termination, without using parameters that change over time. Future extensions of the current study could also investigate the conditions of the seizure onset mechanism by explicitly analyzing the non-stationary dynamics of the network to determine whether the seizure onset is likely to occur through a deterministic parameter change like a bifurcation or is a jump phenomenon due to a noise-driven transition between bistable attractors. The Bayesian inversion in the current study is based on self-tuning algorithms such as NUTS and ADVI achieved by fast automatic differentiation for gradient calculation. This enables efficient sampling from complex high-dimensional posterior distributions using correlation parameters, compared to conventional sampling algorithms. Inference in the presented framework is also enhanced with several MCMC convergence diagnostics to evaluate the reliability of the estimates.

[0083] To improve the preoperative evaluation in the identification of EZ and thus increase the success rate of surgery, various non-invasive and invasive methods are being used. To adopt the BVEP model in clinical treatment and brain intervention, it is necessary to quantify the results of the model when fitting empirical secondary functional signals of patients such as EEG, MEG, SEEG, and fMRI signals. In this framework, it is easy to incorporate additional knowledge such as MRI lesions and clinical hypotheses regarding EZ from preoperative evaluations. Since the BVEP model can be regarded as a general approach to large-scale brain modeling, this provides a promising means for inference from invasive measurements such as clinically used non-invasive imaging signals (EEG, MEG, fMRI) and SEEG signals. The results show that the proposed method according to the present invention can fit well to the empirical SEEG data of patients (not shown). Note that in the case of empirical SEEG recordings, source localization is an ill-posed problem due to the sparsity of the lead field matrix that can affect the accuracy of the estimation. In principle, there is a possibility that the BVEP model can be used to systematically test surgical strategies, but the actual clinical application needs to be investigated and verified in future research.

[0084] In conclusion, the present invention establishes a link between probabilistic modeling and personalized brain network modeling to systematically predict the onset location of seizures in virtual epilepsy patients. Based on large-scale brain network models derived from non-invasive structural data of individual patients, it is demonstrated step by step how the proposed framework can infer epileptogenic spatial maps. The present invention is based on an advanced and efficient sampling algorithm that provides accurate and reliable estimates verified by post-behavioral analysis and convergence diagnosis. In summary, with the help of PPL, using personalized brain network models provides appropriate guidance for testing comprehensive clinical hypotheses and developing new surgical interventions.

Claims

1. A computer-implemented method for inferring the epileptogenicity of brain regions in the seizure activity of the brain of an epilepsy patient, the computer-implemented method being automatically performed by a computer commanded by computer software, the computer-implemented method comprising: providing a computerized model that models various regions of the primate brain and the connectivity between said regions; providing to the computerized model a model capable of reproducing the dynamics of epileptic seizures in the primate brain, said model being a function of parameters of the epileptogenicity of a region of said brain; providing structural data of the brain of an epilepsy patient and using said structural data to personalize the computerized model in order to obtain a virtual epilepsy patient (VEP) brain model; using probabilistic state transitions to convert the state space representation of the virtual epilepsy patient (VEP) brain model into a probabilistic programming language (PPL) in order to obtain a probabilistic virtual epilepsy patient brain model (BVEP); acquiring electroencephalogram data or magnetoencephalogram data of the brain of a patient and fitting the probabilistic virtual epilepsy patient brain model to the electroencephalogram data or magnetoencephalogram data of the brain of said patient in order to infer the epileptogenicity of brain regions in the seizure activity of the brain of said patient; A computer-implemented method comprising the above.

2. The probabilistic programming language is a Bayesian programming language, the probabilistic virtual epilepsy patient brain model is a Bayesian virtual epilepsy patient (BVEP) brain model, and the epileptogenicity of the brain region is inferred using Bayesian inference. The computer-implemented method according to claim 1.

3. The structural data of the brain of the epilepsy patient includes non-invasive T1-weighted imaging data and / or diffusion MRI image data. The computer-implemented method according to claim 1 or 2.

4. The model capable of reproducing the dynamics of epileptic seizures in the primate brain is a model that reproduces the dynamics of the start, progression, and end of seizure events, and includes state variables that couple two oscillatory dynamical systems on three different timescales: the fastest timescale on which the state variable describes high-frequency discharges during the seizure state of the seizure period, the intermediate timescale on which the state variable represents slow spike-wave oscillations, and the slowest timescale on which the state variable is involved in the transition between the seizure-free state and the seizure state. The degree of epileptogenicity of a brain region is represented through the value of an excitability parameter. The computer-implemented method according to claim 1, 2, or 3.

5. To obtain the probabilistic virtual epileptic patient brain model, a spatial map of the epileptogenicity of the patient's brain is provided, and the spatial map of the epileptogenicity classifies the brain regions of the patient's brain into epileptogenic zones (EZs) that may autonomously trigger epileptic seizures, propagation zones (PZs) that do not autonomously trigger seizures but may appear during the progression of seizures, and normal zones (HZs) that do not autonomously trigger seizures. The computer-implemented method according to one of claims 1 to 4.

6. The probabilistic virtual epileptic patient brain model is generated according to a generative model based on the state space representation of the virtual epileptic patient. The computer-implemented method according to one of claims 1 to 5.

7. The state space representation of the virtual epileptic patient is 【Number 1】 in the form of where 【Number 2】 is an n-dimensional vector of the state of a system that changes over time, and x ｔ０ is the initial state vector at time t = 0, and [Number 3] includes all unknown parameters of the virtual epileptic patient model, u(t) represents the external input, 【Number 4】 represents the measurement data targeted by the measurement error v(t), f is a vector function describing the dynamic characteristics of the system, and h represents the measurement function. The computer-implemented method according to claim 6.

8. To obtain the probabilistic virtual epileptic patient (BVEP) model, the state space representation of the virtual epileptic patient (VEP) model is incorporated into the probabilistic virtual epileptic patient (BVEP) model as the state transition probability. The computer-implemented method according to one of claims 1 to 7.

9. The state transition probability is as follows: 【Number 5】 where 【Number 6】 represents the transition probability from state x(t) to state x(t + dt). The computer-implemented method according to claim 8.

10. The generative model is defined in terms of the likelihood of the model parameters and the prior density, and their product is the joint density: 【Number 7】 resulting in where the prior density 【Number 8】 includes prior beliefs regarding hidden variables and potential parameter values, while the conditional likelihood term 【Number 9】 is the computer-implemented method according to claim 6, which represents the probability of obtaining observed values with a given set of parameter values. [

11. ] The computer-implemented method according to one of claims 1 to 10, wherein a sampling algorithm is implemented to infer the epileptogenicity of a brain region in the seizure activity of the patient's brain. [

12. ] The computer-implemented method according to claim 11, wherein the sampling algorithm is a Markov chain Monte Carlo or variational inference algorithm.

Images