Method for inferring epileptogenicity of brain regions

Abstract

The invention relates to a method for inferring epileptogenicity of a brain region not observed to be recovered or not observed to be not recovered in a seizure activity of a brain of an epilepsy patient, comprising the steps of: providing a computerized model modeling individual regions of a primate brain and connectivity between said regions; providing said computerized model with a model capable of reproducing the dynamics of a seizure in a primate brain; providing structural data of a brain of an epilepsy patient and using said structural data to individualize the computerized model in order to obtain a virtual epilepsy patient (VEP) brain model; translating a state space representation of the virtual epilepsy patient (VEP) brain model into a probabilistic programming language (PPL) using probabilistic state transitions in order to obtain a probabilistic virtual epilepsy patient brain model (BVEP); and acquiring electroencephalogram or magnetoencephalogram data of the brain of the patient and fitting the probabilistic virtual epilepsy patient brain model against said data in order to infer the epileptogenicity of said brain region not observed.

Description

Technical Field

[0001] This invention relates to a method for inferring the probability of epileptogenicity of brain regions that were not observed to be recruited or not recruited during seizures in the brains of epileptic patients. Background Technology

[0002] Model inversion (i.e., finding the set of model parameters that produces the best possible fit to observed data) is a challenging task in statistical inference. The Bayesian framework provides a powerful and principled approach to parametric inference and model prediction from experimental data with a wide range of applications. In the context of neuroimaging, the Bayesian approach has been widely used to infer intrinsic parameters of neuronal populations and / or interactions between neuronal populations in pre-specified neural networks from neurophysiological data.

[0003] It is well known that gradient-free sampling algorithms (e.g., Metropolis-Hastings, Gibbs sampling, and slice sampling) often fail to efficiently probe the parameter space when applied to large-scale inverse problems, as is frequently encountered in applications such as whole-brain imaging for clinical diagnosis. In particular, traditional Markov chain Monte Carlo (MCMC) behaves poorly in high-dimensional parameter spaces involving correlated variables. In contrast, gradient-based algorithms (e.g., Hamiltonian Monte Carlo (HMC)), while computationally expensive, are far superior to gradient-free sampling algorithms in terms of the number of individual samples produced per unit of computation time. These sampling algorithms provide efficient convergence and exploration of the parameter space, even in high-dimensional spaces that can exhibit strong correlations. However, the efficiency of gradient-based sampling methods (e.g., HMC) is highly sensitive to user-specified algorithm parameters. More advanced MCMC sampling algorithms (e.g., No-U-Turn Samplers (NUTS) – a self-tuned variant of HMC) address these problems by adaptively tuning the algorithm parameters. These algorithms have been shown to efficiently sample from high-dimensional target distributions that allow for solving complex inverse problems conditioned on large datasets of observations.

[0004] MCMC has the advantage of being nonparametric and asymptotically accurate in the limits of long / infinite runs. Among other alternatives, variational inference (VI) transforms Bayesian inference into an optimization problem, which typically yields much faster computation than MCMC methods. However, the classical derivation of VI requires model-specific work on defining a variational series suitable for the probabilistic model, computing the corresponding objective function, calculating the gradient, and running a gradient-based optimization algorithm on the master model. Automatic Differential Variational Inference (ADVI) solves these problems automatically.

[0005] Probabilistic programming languages ​​(PPLs) provide efficient implementations of automated Bayesian inference on user-defined probabilistic models, characterized by generating next-generation MCMC sampling and VI algorithms (e.g., NUTS and ADVI). With PPLs, these algorithms utilize automatic differentiation methods to compute derivatives in computer programs, avoiding random walk behavior and sensitivity to relevant parameters. In particular, Stan and PyMC3 are high-level statistical modeling tools for Bayesian inference and probabilistic machine learning, providing high-level inference algorithms rich in extensive and reliable diagnostics, such as NUTS and ADVI. While PPLs allow for automated inference, the performance of these algorithms can be sensitive to the form of parameterization. Reparameterizing probabilistic models to improve the inference efficiency of system dynamics (managed by a set of nonlinear stochastic differential equations) in appropriate forms remains a challenging problem.

[0006] On the other hand, personalized large-scale brain network modeling has gained popularity in recent years due to its potential to improve medical treatment strategies. In personalized whole-brain modeling, patient-specific information (e.g., anatomical connectivity) obtained from non-invasive imaging techniques is combined with mean-field models of local neuronal activity to simulate an individual's spatial-temporal brain activity on a macroscopic scale. Virtual Brain (TVB) is an open-access computational framework written in Python that reproduces and evaluates personalized configurations of the brain using individual subject data. This neuroinformatics platform integrates computational brain modeling and multimodal neuroimaging data to systematically simulate an individual's spatiotemporal brain activity. However, currently, there is no specific workflow for automated model inversion and data fit validation prepared for TVB.

[0007] A novel approach to brain intervention based on personalized brain network models derived from non-invasive structural data of individual patients has recently been proposed: Virtual Epilepsy Patient (VEP). VEP models are large-scale computational models of the individual brain that incorporate personal data (e.g., location of seizure onset, subject-specific brain connectivity, and MRI lesions) to inform specific clinical monitoring and improve surgical outcomes. Previous studies have shown that VEP models can realistically simulate the evolution of seizures in patients with bilateral temporal lobe epilepsy. However, the inverse problem of such large-scale brain network models is a challenging task due to the intrinsic nonlinear dynamics of each brain network node, the large number of associated model parameters, and observations typically encountered in brain imaging settings. Summary of the Invention

[0008] Accordingly, there is a need to establish a useful link between the most popular probabilistic programming tools (e.g., Stan / PyMC3) and personalized brain network modeling (e.g., VEP models) in order to systematically predict the location of seizure onset in virtual epilepsy patients. This invention specifically allows for the construction of a Bayesian Virtual Epilepsy Patient (BVEP) as a probabilistic framework designed to infer the hidden / unobserved dynamics of a personalized, large-scale brain model of epilepsy spread generated by TVB.

[0009] According to a first aspect, the present invention relates to a method for inferring the epileptogenicity of brain regions that were not observed to be restored or not observed to be unrestored during seizure activity in the brain of an epileptic patient, comprising the following steps:

[0010] Provides a computerized model for modeling the various regions of the primate brain and the connectivity between said regions;

[0011] The computerized model is provided to reproduce the dynamics of epileptic seizures in the primate brain, the model being a function of parameters of epileptogenicity of brain regions;

[0012] Structural data of the brains of epilepsy patients are provided, and the structural data is used to personalize computerized models in order to obtain virtual epilepsy patient (VEP) brain models;

[0013] The state-space representation of the Virtual Epilepsy Patient (VEP) brain model is transformed into a probabilistic programming language (PPL) using probabilistic state transitions to obtain a probabilistic Virtual Epilepsy Patient (BVEP) brain model; and

[0014] Obtain electroencephalogram (EEG) or magnetoencephalogram (MEG) data of the patient's brain, and fit the data to a probabilistic virtual brain model of an epileptic patient in order to infer the epileptogenicity of the brain regions that were not observed to be restored or not observed to be restored during the patient's seizure activity.

[0015] Prior to this, - 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 Bayesian inference is used to infer the epileptogenicity of brain regions not observed as relapsed or not observed as unrelapsed; - the structural data of the epilepsy patient brain includes non-invasive T1-weighted imaging data and / or diffuse MRI image data; - the model capable of reproducing the seizure dynamics of the primate brain is a model that reproduces the onset, progression, and cancellation seizure events, including state variables coupled to two oscillatory dynamic systems at three different time scales: the fastest time scale, where the state variables record the rapid discharge during the burst state; the intermediate time scale, where the state variables represent slow spikes and wave oscillations; and the slowest time scale, where the state variables are responsible for the transition between the inter-burst period and the burst state, and where the degree of epileptogenicity of brain regions is represented by the value of an excitability parameter; - To obtain a probabilistic virtual epilepsy patient brain model, a spatial map of the patient's brain's epileptogenicity is provided, which classifies the brain regions into epileptogenic zones (EZs) (which can spontaneously trigger seizures), propagation zones (PZs) (which do not spontaneously trigger seizures but can recover during seizure evolution), and healthy zones (HZs) (which do not spontaneously 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 has the following form

[0016]

[0017] in, It is an n-dimensional vector of the system state that evolves over time, x t0 It is the initial state vector at time t = 0. Includes all unknown parameters of the virtual epilepsy patient model, where u(t) represents the external input. Let f represent the measurement data subject to measurement error v(t), f be a vector function describing the dynamic properties of the system, and h represent the measurement function; - To obtain the 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 state transition probabilities; - State transition probabilities are, for example:

[0018]

[0019] in, This represents the transition probability from state x(t) to x(t + dt); - The generative model is defined based on the likelihood and prior model parameters (their product produces the following joint density):

[0020]

[0021] Among them, the prior distribution This includes prior beliefs about the values ​​of hidden variables and latent parameters, while the conditional likelihood term... - To obtain the probability of observation using a given set of parameter values; - To implement a sampling algorithm to infer the epileptogenicity of brain regions that were not observed to be relapsed or not observed to be relapsed during seizure activity in a patient's brain; - The sampling algorithm is a Markov chain Monte Carlo or variational inference algorithm; and - The method is implemented by a computer. Attached Figure Description

[0022] Other features and aspects of the invention will become apparent from the following description and accompanying drawings, wherein:

[0023] Figure 1 This is a schematic diagram of the method according to the present invention;

[0024] Figure 2A , Figure 2B , Figure 2C , Figure 2D and Figure 2E The results of spatial maps obtained according to the method of the present invention for estimating epileptogenicity across different brain regions of a patient are shown. More specifically, Figure 2A This shows the segmentation of the patient's reconstructed brain. Figure 2B The patient's brain network, consisting of 84 regions, is shown (gray: HZ, light gray: PZ, dark gray: EZ). The thickness of the lines indicates the strength of the connections. For illustrative purposes, only the 10% of connections with weights higher than the maximum weight are shown. Figure 2C The structural connectivity matrix is ​​shown. Figure 2D This is a demonstration simulation of a complete VEP model of source-level brain activity and prediction envelope (dashed line). Figure 2E Excitability parameters for different brain node types are shown. The estimated density, where the vertical dashed line indicates the true value;

[0025] Figure 3A , Figure 3B , Figure 3C ,and Figure 3D The accuracy of the results obtained by the method according to the invention is shown, such as the spatial map of epileptogenicity estimation across different brain regions using the NUTS algorithm for patient 1. More specifically, Figure 3A The diagram shows examples of observed data (dotted lines) and predictions for three brain node types defined as HZ (grey), PZ (light gray), and EZ (dark gray). The shaded areas depict the range between the 5th and 95th percentiles of the posterior prediction distribution. Figure 3B It shows the indication of 84 brain regions. A plot of the estimated density. The actual values ​​are shown by solid black circles. Figure 3C The distribution of posterior z-scores and posterior contraction is shown, which implies an ideal Bayesian inversion. Figure 3D The confusion matrix of the estimated spatial map of epileptogenicity is shown. Predicts accurately the predefined classes of all brain nodes labeled HZ, PZ, and EZ (accuracy = 1.0, misclassification = 0.0).

[0026] Figure 4 This diagram compares simulated (top row) and predicted (bottom row) phase planes for different brain node types in the BVEP model. From left to right, the columns correspond to brain nodes designated as HZ, PZ, and EZ, respectively. The trajectories of these brain regions are shown in green, yellow, and red, respectively. In each phase plane, the intersection of the x and z zero tidal lines (colored in dark gray), depending on the excitability parameter, determines the fixed points of the system. Full circles and empty circles indicate stable and unstable fixed points, respectively.

[0027] Figure 5 The diagram shows the estimated spatial map of epileptogenicity obtained by the NUTS algorithm compared to ADVI. The exemplary histograms and kernel density estimates of the samples obtained by NUTS are shown in plane A, compared to the approximation of the mean-field variant obtained by ADVI shown in plane B. For all brain nodes included in the analysis, the prior (shown in dark gray) is assumed to be N(-2.5, 1.0). Vertical dashed lines indicate the true values; and...

[0028] Figure 6 The diagram illustrates the convergence diagnostics of NUTS and ADVI. In plane A, samples generated by NUTS from the joint posterior probability distribution between hyperparameter pairs (σ, σ') are shown. In this case, the parameterized non-centered form generates individual samples from the posterior distribution. In plane B, the centered form of sampling results in high correlation between hyperparameters, indicating that the sampler is not effectively probing the posterior distribution. In plane C, samples from an approximate joint posterior probability distribution are shown using the mean-field variant of ADVI. Plane D shows the non-centered form of sampling. A value below 1.05 for all estimated hidden states and parameters indicates that the MCMC has converged. In plane E, the value returned by the centered sampling method... The high-value indicator chain has not yet converged. In plane F, the variational objective function (ELBO) and the number of ADVI iterations are shown. Detailed Implementation

[0029] This invention relates to a method for inferring the epileptogenicity of brain regions that were not observed to be remission or not observed to be non-remission during seizure activity in the brains of epileptic patients. This is a computerized probabilistic method for inferring the spatial map of epileptogenicity of different brain regions in individualized epileptic patients whose seizures begin across hypothetical regions and can propagate to candidate brain regions.

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

[0031] It includes steps for providing computerized models that model the various regions of the primate brain and the connectivity between said regions. This brain is a virtual brain. It is a neuroinformatics platform used for whole-brain network simulation using biorealistic connectivity. This simulation environment enables model-based reasoning across neurophysiological mechanisms at different brain scales, which are the basis for generating macroscopic neuroimaging signals, including functional magnetic resonance imaging (fMRI), EEG, and magnetoencephalography (MEG). It allows for the reproduction and evaluation of personalized configurations of the brain using individual subject data.

[0032] It further includes the step of providing the computerized model with a model capable of reproducing the dynamics of epileptic seizures in the primate brain, the model being a function of parameters of the epileptogenicity of brain regions.

[0033] Preferably, a model capable of reproducing the epileptic seizure elevator in the primate brain is a model that reproduces the dynamics of the onset, progression, and cancellation of seizure events, comprising state variables of two oscillating dynamic systems coupled at three distinct timescales: the fastest timescale, where the state variables describe the rapid discharge during the burst state; the intermediate timescale, where the state variables represent slow spikes and wave oscillations; and the slowest timescale, where the state variables are responsible for the transition between the inter-burst period and the burst state, and where the epileptogenicity of brain regions is represented by the value of an excitability parameter.

[0034] Furthermore, the present invention also includes the steps of: providing structural data of the brain of an epileptic patient, and using said structural data to personalize a computerized model to obtain a virtual epileptic patient (VEP) brain model. The structural data is, for example, image data of the patient's brain acquired using magnetic resonance imaging (MRI), diffusion-weighted magnetic resonance imaging (DW-MRI), or nuclear magnetic resonance imaging (NMRI) computed tomography (MRT). Preferably, the structural data of the epileptic patient's brain includes non-invasive T1-weighted imaging data and / or diffusion MRI image data.

[0035] The method according to the invention further includes the step of: 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).

[0036] Prior to this, 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 Bayesian inference is used to infer the epileptogenicity of unobserved (neither unobserved as restored nor unobserved as unrestored) brain regions.

[0037] Prior to obtaining a probabilistic virtual brain model of an epilepsy patient, an epileptogenic spatial map of the patient's brain is provided, which classifies brain regions of the patient's brain into epileptogenic zones (EZ) (which can autonomously trigger seizures), propagation zones (PZ) (which do not autonomously trigger seizures but are able to recover during seizure evolution) and healthy zones (HZ) (which do not autonomously trigger seizures).

[0038] Prioritize generating a probabilistic virtual epilepsy patient brain model based on a generative model of the state space representation of a virtual epilepsy patient.

[0039] More preferably, the state-space representation of a virtual epilepsy patient has the following form

[0040]

[0041] in, It is an n-dimensional vector of the system state that evolves over time, x t0 It is the initial state vector at time t = 0. Includes all unknown parameters of the virtual epilepsy patient model, where u(t) represents the external input. Let f represent the measurement data subject to measurement error v(t), f be a vector function describing the dynamic properties of the system, and h be the measurement function.

[0042] Prior to obtaining the 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 the state transition probability.

[0043] More preferably, the state transition probabilities are, for example:

[0044]

[0045] in, This represents the probability of transition from state x(t) to x(t + dt).

[0046] More preferably, the generative model is defined based on the likelihood and prior model parameters (their product yielding the following joint density):

[0047]

[0048] Among them, the prior distribution This includes prior beliefs about the values ​​of hidden variables and latent parameters, while the conditional likelihood term... This indicates the probability of obtaining an observation using a given set of parameter values.

[0049] The method according to the invention further includes the steps of: acquiring electroencephalogram (EEG) or magnetoencephalogram (MEG) data of the patient's brain, and fitting a probabilistic virtual epilepsy patient brain model against the data in order to infer the excitability of the brain regions not observed (neither observed as recovery nor observed as non-recovery) during the patient's seizure activity.

[0050] Prior to this, a sampling algorithm is implemented to infer the excitability of brain regions that were not observed to be restored or not observed to be restored during the patient's seizure activity.

[0051] More preferably, the sampling algorithm is Markov chain Monte Carlo or variational inference algorithm.

[0052] Example 1: Materials and Methods

[0053] In the example, the method according to the invention is based on personalized brain network modeling and Bayesian inference, such as... Figure 1As shown in the diagram. This approach allows for the construction of a BVEP based on two main steps that constitute a VEP model, and then the VEP is embedded in a PPL tool to infer and validate model parameters. To construct the VEP model, the following steps are performed: First, the patient undergoes non-invasive brain imaging (MRI, DTI). Based on these images, a brain network anatomy is provided from the reconstruction pipeline, including brain segmentation and the patient's connectome. Then, a neural population model is selected for each brain region to define the network model. In the VEP, an Epiliptor™ model is defined at each network node connected by structural connectivity derived from diffuse beam imaging. Such a model is disclosed, for example, in a publication entitled "On the nature of seizure dynamics" (Jirsa et al., Brain 2014, 137, 2210-2230), which is incorporated herein by reference. It comprises five state variables that operate at three different timescales. At the fastest timescale, the state variables record the rapid discharges during an attack. At the slowest timescale, capacitive state variables describe slow processes, such as changes in extracellular ion concentration, energy expenditure, and tissue oxygenation. The system exhibits rapid oscillations during the burst state via these state variables. Autonomous switching between the burst interval and the burst state is achieved via capacitive state variables through the onset of the attack and the canceling saddle node and isoclinic bifurcation mechanisms, respectively. Switching is accompanied by a direct current (DC) offset, which has been recorded in vitro and in vivo. At intermediate timescales, other state variables describe the spikes and electrocochleography patterns observed during the attack, as well as the burst interval and pre-burst spike when excited by the fastest system via coupling. In summary, TVB simulation allows for the simulation of empirical neuroimaging signals. Model fitting is then performed using, for example, the NUTS / ADVI algorithm within the PPL tool (in this example, as observed brain-derived activity and as a VEP model transformed into a generative model in Stan / PyMC3). Finally, cross-validation can be performed from existing samples, for example via WAIC / LOO, to evaluate the model's ability to predict on new data, thus refining the network pathology. In other words, the workflow for constructing the BVEP consists of two main steps: constructing the VEP, i.e., a personalized brain network model of epileptiform spread; and then embedding the VEP model into a Bayesian framework to infer and validate model parameters. After the VEP is formulated in the state-space representation, the probabilistic reparameterization of the system dynamics is demonstrated. It is shown that the proposed probabilistic reparameterization in the BVEP can effectively invert the nonlinear state-space equations to infer system dynamics. This approach allows for accurate estimation of the epileptogenic spatial graph in the personalized brain network model of epileptiform spread by leveraging PPL. A virtual brain is used for brain network simulation, and Stan and PyMC3 are used to invert the simulated whole-brain model.

[0054] The following sections step-by-step demonstrate how to construct a BVEP model for a specific patient, fit the constructed brain model against in-silico data, and validate our reasoning. The accuracy and reliability of the estimates are validated through several convergent diagnostic and posterior behavioral analyses.

[0055] Individual patient data

[0056] For 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 prefrontal lobe epilepsy (Patient 2). Patients underwent standard clinical evaluation, details of which are described in a previous study (Proix et al., Individual brain structure and modelling predict seizure propagation. Brain 140, 641-654, 2017). 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 diffuse MRI (DTI-MR sequence, angular gradient ensemble in 64 directions, repetition time = 10.7 s, echo time = 95 ms, 2.0 × 2.0 × 2.0 mm, 70 slices, 1000 s mm). -2 (b-weighted). Images acquired on a Siemens Magnetom Verio™ 3T MR scanner.

[0057] Network Anatomy

[0058] The reconstructed connectome is built using a reconstruction pipeline employed with commonly available neuroimaging software. The current version of this pipeline is an evolution of a previously described version (Proix et al., 2017). First, the command `recon-all` from the Freesurfer™ package in version v6.0.0 is used to reconstruct and segment brain anatomy from T1-weighted images. Then, the T1-weighted images are co-registered with the diffusion-weighted images using the linear registration tool `flirt™` from version 6.0, employing a correlation ratio cost function with 12 degrees of freedom.

[0059] The MRtrix™ package in version 0.3.15 is used for beam imaging. Fiber orientation is estimated from DWI using spherical deconvolution via the dwi2fod tool, where the response function is estimated using the Tournier algorithm via the dwi2response tool. Subsequently, 15 million fiber bundles are generated using the iFOD2 probabilistic beam imaging algorithm via the tckgen tool. Finally, a connectivity group matrix is ​​constructed using the Desikan-Killiany segmentation generated by FreeSurfer in the previous step via the tck2connectome tool. The connectivity groups are normalized such that the maximum value equals one.

[0060] Network Model

[0061] Typically, to construct personalized brain network models, segmentation schemes are used to define brain regions, and a set of mathematical equations is used to model regional brain activity. This data-driven approach incorporates subject-specific brain anatomical information, with network edges represented by structural connectivity derived from non-invasive imaging data of individual patients. In the VEP model, the dynamics of brain network nodes are governed by the Epiliptor equation (Jirsa et al., On the nature of seizure dynamics. Brain 137, 2210-2230, 2014), and these nodes are coupled through a structural connectivity matrix derived from diffusion-weighted MRI (dMRI) (Jirsa et al., The virtual epileptic patient: Individualized whole-brain models of epilepsy spread, 2017).

[0062] Epileptor is a dynamic model of seizure evolution, realistically reproducing the dynamics of onset, progression, and the initial onset of seizure-like events. Epileptor comprises five state variables coupled to two oscillatory dynamic systems across three distinct timescales: at the fastest timescale, variables x1 and y1 describe the rapid discharge during the burst state. At the intermediate timescale, variables x2 and y2 represent slow spikes and wave oscillations. At the slowest timescale, the capacitive state variable z is responsible for the transition between the burst interval and the burst state. Additionally, the burst interval and pre-burst spike are generated via the term g(x1).

[0063] According to Jirsa et al. (2014), the dynamics of the complete Epipuleptor model are described by the following equation:

[0064]

[0065] in

[0066]

[0067] Where τ0 = 2857, τ1 = 1, τ2 = 10, Ι1 = 3.1, Ι2 = 0.45, and γ = 0.01. The degree of epileptogenicity is represented by the value of the excitability parameter η. If , where η C If the threshold for epileptogenicity is reached, then the epileptor will spontaneously exhibit seizure activity and is called epileptogenic; otherwise, the epileptor is in its (healthy) equilibrium state and does not spontaneously trigger seizures.

[0068] According to Jirsa et al. (2017), the complete VEP brain model equation (N-coupled Epiliptor) is as follows:

[0069]

[0070] The network nodes are approximated by a linear method using permeability coupling. (It includes the global scaling factor K and the patient's connectome C) ij They are coupled.

[0071] Under the assumption of timescale separation, Proix et al. (2014) have shown that by averaging the coupled Epileptor equations (which yields a 2D reduction of the VEP model as follows), the second set of Epileptor neurons (i.e., variables x2 and y2) can be neglected:

[0072]

[0073] Depending on the value of the excitability parameter η, the 2D Epileptor exhibits different levels of stability. For The trajectory in the phase plane is attracted to a single stable fixed point of the system on the left branch of the cubic x-zero slope line (nullcline). In this system, the epipuleptor is said to be healthy, meaning that it does not trigger seizures in the absence of external input. With As the value increases, the zero slope of z shifts downward, and the saddle node bifurcates at the point corresponding to the onset of the attack. It happened. For The system exhibits unstable fixed behavior.

[0074] This allows for the occurrence of seizures (which Epileptor calls epileptogenic). In this example, we use a 2D reduction of the VEP model for Bayesian inference of the epileptogenic spatial map to reduce the computational cost associated with model parameter estimation. The 2D reduction of Epileptor allows for faster inversion while enabling us to predict the envelope of rapid emissions during burst seizure states (i.e., the onset, progression, and cancellation of seizure patterns) (Proix et al., 2014; Jirsa et al., 2017).

[0075] Epileptigenic Spatial Map

[0076] In addition to the patient's connectome (which structurally restricts individual variability), the dynamics of brain network models can be further constrained by formulating hypotheses related to functional brain network components to produce more specific patterns of brain activity across individuals. In the case of epilepsy, clinical hypotheses about the location of epileptogenic zones or lesions allow for the subdivision of network pathology to better predict the onset and progression of seizures in individual patients.

[0077] In the BVEP brain model, each network node can trigger a seizure based on its connectivity and excitability values. The parameter η controls tissue excitability, and its spatial distribution is therefore the target for parameter fitting. In this study, depending on the excitability value, different brain regions are classified into three main types:

[0078] - Epileptogenic Zone (EZ): If Then the epileptor can autonomously trigger seizures (the brain region responsible for the origin and early organization of epileptic activity).

[0079] - Propagation Zone (PZ): If If the epipileptor does not trigger a seizure spontaneously, it can recover during the seizure's evolution because its equilibrium state is close to the critical value.

[0080] - Health Zone (HZ): If Then Epipuleptor does not spontaneously trigger seizures.

[0081] Based on the above dynamic properties, the spatial map of epileptogenicity across different brain regions includes the excitability values ​​of the EZ (high excitability value), PZ (lower excitability value), and all other regions classified as HZ (non-epileptic). However, it should be noted that intermediate excitability values ​​do not guarantee that a seizure will revert to a part of the propagation zone, as propagation is also determined by various other factors, including connectivity and brain state dependence. In the BVEP brain model, clinical hypotheses can be formulated as prior knowledge about the spatial distribution of excitability parameters. In this study, it is assumed that there are no clinical hypotheses about specific brain regions, and the same prior distribution is assigned to excitability parameters across all brain regions included in the analysis.

[0082] Probabilistic Model

[0083] A key component in probabilistic brain network models within the Bayesian framework is the generative model. Given a set of observations, the generative model is a probabilistic description of the mechanism by which it generates the observed data through some hidden states and unknown parameters. Here, the generative model will therefore have a mathematical formula guided by a dynamic model that describes the evolution of the model's state variables and the given parameters over time. This specification is required to construct the likelihood function. The complete generative model is then completed by specifying prior beliefs about the possible values ​​of the unknown parameters.

[0084] The BVEP brain model presented in this study was established in two main steps. First, the VEP model equations provide the basic form of the data generation process describing how epileptic seizures occur. Second, hypotheses about the spatial map of epileptogenicity in the brain were formulated as prior knowledge. The latter component informs the model using hypotheses about the spatial distribution of excitability parameters across different brain regions.

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

[0086]

[0087] in, It is an n-dimensional vector of the system state that evolves over time, x t0 It is the initial state vector at time t = 0. Includes all unknown parameters of the virtual epilepsy patient model, where u(t) represents the external input. This represents the measurement data subjected to measurement error v(t). (The data is obtained through...) and The process (dynamic) noise and measurement noise are assumed to follow a mean of zero and a variance of zero, respectively. and The Gaussian distribution of the system is used. Colored and non-Gaussian dynamic noise can be captured in the term w(t), while additional terms appear in the presence of multiplicative noise (i.e., noise whose intensity depends on the system state) or multiplicative feedback (where the system state further influences the intensity of the driving noise), which can lead to solutions of different quality. Furthermore, f(.) is a vector function describing the dynamic properties of the system, and h(.) denotes the measurement function. In the source localization problem, h(.) is called the guiding field matrix. Note that current work focuses on observing the potential brain sources of activity to avoid the inevitable inconsistencies associated with the mapping from the source dipole to the electrode contact (i.e., h(.) is a linear function here).

[0088] Considering the 2D reduction of the VEP model (see equation (3)), then , where n = 2N, and N equals the number of brain regions. Accordingly, Where p = 3N + 3. The patient is virtualized using a reconfiguration pipeline, as described in Section 2.2, where N = 84.

[0089] The dynamic state-space representation of the hidden state x(t) is defined (see equation (4) and incorporated into the BVEP model as the state transition probability):

[0090]

[0091] in, This represents the probability of transition from state x(t) to x(t + dt). However, the above parameterization, called central parameterization, can exhibit a pathological geometry that produces biased estimates.

[0092] It has been previously shown that careful selection of reparameterization increases the effective sample size and reduces divergence, especially for regions with extreme curvature. To avoid pathological samples and thus biased estimations caused by strong correlations between parameters in the centralized form of the parameterization, the advantages of position-scale transformation are used to invert nonlinear state-space equations, which allows for uncorrelation of the parameters representing state variables at continuous time steps.

[0093] The non-centered reparameterization of the above distribution is described below:

[0094]

[0095] Example 3 shows that by avoiding biased estimations caused by strong correlations between parameters, using a parameterized non-centered form to infer system dynamics greatly improves sampling performance.

[0096] Reasoning / Prediction

[0097] Generative models are characterized by the joint probability distribution of model parameters and observations. Where Y represents the observed variable, and This includes the system's hidden variables and model parameters. Given only the observed responses and prior beliefs related to the underlying generation process, Bayesian techniques infer the distribution of unknown parameters of the underlying data generation process. Through the product rule, the generative model can be defined based on the likelihood and prior model parameters (their product yielding a joint density as follows):

[0098]

[0099] Among them, the prior distribution This includes prior beliefs about the values ​​of hidden variables and latent parameters, while the conditional likelihood term... This represents the probability of obtaining an observation using a given set of parameter values. In Bayesian inference, this involves seeking the posterior density. It is the conditional distribution of the model parameters given the observations. Bayes' theorem expresses this posterior density based on likelihood and prior information as follows:

[0100]

[0101] Among them, the denominator It represents the probability of the data, and it is called the evidence or marginal likelihood (which is actually just a normalization term).

[0102] In order to obtain the posterior density Sampling reveals that the performance of HMC is highly sensitive to the step size and number of steps in the leapfrog integrator used to update the position and momentum variables in the Hamiltonian dynamics simulation. If the number of steps in the leapfrog integrator is chosen too small, HMC exhibits undesirable random walk behavior similar to the Metropolis-Hastings algorithm, and thus the algorithm poorly explores the parameter space. If the number of leapfrog steps is chosen too large, the associated Hamiltonian trajectory may loop back to the neighborhood of the initial state, and the algorithm wastes computational effort. NUTS extends HMC by adaptively tuning the step size and number of steps in the leapfrog integrator to efficiently sample from the posterior distribution. In the alternative, ADVI assumes a set of densities, automatically calculates the gradient, and then finds the closest member (measured by the Kullback-Leibler divergence). In this study, we use NUTS, a self-tuned variant of HMC, and ADVI to approximate the posterior distribution of the model parameters (see Equation (3)).

[0103] Prior excitability parameters for all brain regions included in the analysis were assumed to follow a normal distribution with a mean of -2.5 and a standard deviation of 1.0, i.e., N(-2.5, 1.0). Furthermore, weakly informative priors were imposed on the system initial conditions and the global coupling parameter K, assuming a ground-centric normal distribution with a standard deviation of 1.0. Priors on hyperparameters were considered to be the general weakly informative prior N(0, 1.0).

[0104] After fitting a Bayesian model, it is typically necessary to measure the predictive accuracy of the inferred model. The information criterion and leave-one-out cross-validation (LOO) are two rigorous methods for evaluating a model's ability to predict new data. Broadly applicable information criterion (WAIC) and Pareto smooth importance sampling (PSIS) LOO allow for efficient estimation of the predictive accuracy of the fitted Bayesian model within a negligible computational time relative to the cost of model fitting, leveraging existing simulation extractions from the log-likelihood evaluated in the posterior of the parameter values.

[0105] Reasoning and Diagnosis

[0106] After running the MCMC sampling algorithm, a statistical analysis is required to evaluate the convergence of the MCMC samples. A simple way to evaluate the performance of the MCMC algorithm based on posterior samples is to visualize the degree of chain mixing (i.e., how effectively the MCMC sampler probes all patterns in the parameter space). This can be monitored in different ways, including tracking plots (from the evolution of parameter estimates from iterative MCMC extractions), pairwise plots (identifying collinearity between variables), and autocorrelation plots (measuring the degree of correlation between extractions of MCMC samples). A more quantitative way to evaluate convergence to a stationary MCMC distribution is to estimate the potential scaling factor based on samples with posterior model probabilities. and effective sample size N eff . The diagnostics provide an estimate of how much variance might be reduced by running the chain longer. Each MCMC estimate has its associated... Statistically, this is basically the ratio of inter-chain variance to intra-chain variance. If... If the value is approximately less than 1.1, then MCMC convergence has been achieved (approaching 1.0 in the case of infinite samples); otherwise, the chain will require a longer run time. Furthermore, N... eff The statistics give the number of individual samples represented in the chain. The larger the effective sample size, the higher the accuracy of the MCMC estimate. Note that these are necessary but not sufficient conditions for the convergence of the MCMC sample.

[0107] In addition to the general MCMC diagnostics described above, NUTS-specific diagnostics can also be used to monitor sample convergence; the number of divergent frog jump transitions (due to changes in posterior curvature at height); the step size used by NUTS in its Hamiltonian simulation (if the step size is too small, the sampler becomes inefficient, and 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 frog jump steps taken during the Hamiltonian simulation.

[0108] Evaluation of posterior fit

[0109] Using synthetic data for fitting allows us to validate our inference because the ground reality for the inferred parameters is known. Therefore, the standard error metric can be used to measure the similarity between the inferred parameters and those used to generate the data. The metrics used to validate our inference are the confusion matrix, posterior shrinkage, and posterior z-score.

[0110] The confusion matrix is ​​a metric for evaluating classification accuracy. The element q... i,j This equals the number of observations known to be in class i but predicted to be in class j, where Where Q is the total number of classes. In the BVEP model, we define three groups (i.e., HZ, PZ, and EZ) to classify brain regions, so Q = 3.

[0111] Furthermore, to quantify the precision of the inference, the posterior z-score (represented by z) relative to the posterior contraction (represented by s) is plotted, and they are defined as:

[0112]

[0113] in, and These are the estimated average and the actual ground conditions, respectively. and These indicate the variance (uncertainty) of the prior and posterior, respectively. The posterior z-score quantifies the degree to which the posterior distribution reflects the actual ground reality, while the posterior contraction quantifies the degree to which the posterior distribution shrinks from the initial prior distribution.

[0114] Synthetic datasets and model inversion

[0115] To validate inferences using BVEP, the simulation capabilities of the Virtual Brain (TVB) were utilized to generate synthetic datasets. TVB is an open-source neuroinformatics tool written in Python that simulates 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, which have broad clinical applications ranging from Alzheimer's disease and chronic stroke to focal epilepsy in humans (Jirsa et al., 2017).

[0116] In this study, TVB was used to reconstruct personalized brain network models. To validate inferences about spatial epileptogenicity, epileptic seizures were simulated in two patients: one simulation where the seizure spread to all brain nodes designated as PZ (patient 1), and another simulation where the seizure extended to some brain nodes designated as PZ (patient 2). These datasets were generated using two different structural connectivity matrices and different spatial graphs of epileptogenicity.

[0117] Patient 1's seizure activity was simulated by setting two regions to EZ and three regions to PZ, where, and , where are excitability values ​​respectively. and All other brain nodes were fixed to be non-epileptic, i.e. HZ.

[0118] To simulate the seizure activity of Patient 2, two brain regions were selected as the EZ, and five regions were located at the nodes. and The region is selected as the PZ. For the region selected as the EZ, the excitability value is set to... The excitability of PZ is set to... And all other areas are defined as HZ.

[0119] In both synthetic datasets, the Euler-Maruyama integration scheme is used to model the VEP model as a system of stochastic differential equations, with an integration step size of 0.04. Additive white Gaussian noise is 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), where the mean is zero and the variance is (0.01, 0.01, 0.0, 0.0015, 0.0015). Initial conditions are chosen for each state variable in the interval (-2.0, 5.0).

[0120] Finally, to invert the BVEP of the simulated dataset, two popular open-source PPL tools were used for flexible probabilistic inference: Stan and PyMC3. The Stan language can run in different interfaces, while PyMC3 provides several MCMC algorithms directly in native Python code for model specification. By specifying the model density function in these tools, the gradient of the function is computed via automatic differentiation (i.e., powerful techniques for arithmetic computation of derivatives to efficiently approximate the log-posterior density using NUTS and ADVI). Computation of individual MCMC chains can also be performed in parallel on independent processors. In this example, the Stan command-line interface is used, while all code for simulation and posterior-based analysis is implemented using PyTon. Model simulation and parameter estimation were performed on a Linux machine with a 3.0 GHz Intel Xeon processor and 32GB of memory.

[0121] Example 2: Result

[0122] The results of the workflow for estimating the spatial map of epileptogenicity across different brain regions in Patient 1 using the BVEP model are as follows: Figures 2A-2E The diagram shows the segmentation of the reconstructed brain and the patient's brain network. Figure 2A and Figure 2B As shown in the diagram, after the Desikan-Killiany segmentation used in the reconstruction pipeline, the patient's brain was divided into 68 cortical areas and 16 subcortical structures. Figure 2C The structural connectivity matrix derived from diffusion-beam imaging of the patient is shown. Following the virtualization of the patient's brain, TVB was used to simulate a reconstructed VEP brain network model. Simulation time series of fast-activity variables in the complete VEP brain model are shown in... Figure 2D As shown in the diagram. Different brain node types (i.e., HZ, PZ, and EZ) are coded in green, yellow, and red, respectively. When the Epipuleptor is isolated (i.e., K = 0; no network coupling), seizures are triggered only in the regions defined as EZ, and no seizure propagation is observed in other regions (see [reference]). Figure 2A and Figure 2D However, through the structural connectivity matrix of the patient, Epipuleptor is coupled (see...). Figure 2C Spatial restoration patterns can be observed in candidate brain regions defined as PZ (see...). Figure 2D ). And 2 of the patients who recovered only one of PZ (see Figure 2E In contrast, here, due to the strong coupling connections to the region designated as PZ and the high excitability values ​​of these nodes, the seizure propagates to all other candidate brain regions designated as PZ (node ​​numbers 6, 12, and 28). The mean of the fast activity variables inferred from the VEP model by inversion reduction is shown in the figure. Figure 2DThe data is shown by dashed lines. It can be seen that there are significant similarities between the simulated and predicted seizures related to the onset, spread, and termination of the seizure. Note that the simulation illustrates the activity of the rapid variables in the full VEP brain model (i.e., x in equation (2)). 1,i (t)), while the inferred envelope model of the time series comes from the inverted trajectory of the reduced VEP model (see Equation (3)). The estimated density of excitability parameters for different brain node types is in Figure 2E As shown in the figure, the true values ​​of the excitability parameters (vertical dashed lines) are supported by estimated posterior densities across different brain regions.

[0123] The accuracy of the spatial map estimating epileptogenicity across different brain regions in patient 1, achieved through BVEP in Stan, is [not specified]. Figures 3A to 3D The results are presented in [the PyMC3 implementation]. A similar result can be obtained from the BVEP implementation in PyMC3. Figure 3A The observed and inferred source activities were compared for three brain node types designated as HZ, PZ, and EZ (node ​​numbers 1, 6, and 7, respectively). The simulation data consisted of fast variables from a full VEp brain model sampled at 1000 Hz (i.e., x in equation (2)). 1,i The activity of (t) is composed of 120 s, which is downsampled by a factor of 10 to reduce the computational cost of Bayesian inversion. Observed data are shown by dashed lines, while shaded areas indicate the range between the 5th and 95th percentiles of the posterior prediction distribution. The activity of selected brain nodes in HZ, PZ, and EZ is shown in green, yellow, and red, respectively. It was observed that the predicted time series based on samples from the posterior prediction distribution is in very good agreement with the simulation. Figure 3B A violin plot showing the estimated density of excitability parameters for all 84 brain regions included in the analysis. Solid black circles represent the actual parameter values ​​used to generate the simulation data. It can be seen that the ground truth of excitability parameters for all brain regions is supported by the estimated posterior distribution. Figure 3C As shown, the distributions of the posterior z-scores and posterior contractions for all inferred excitabilities confirm the reliability of the model inversion. Note that the concentration toward large contractions indicates that all posteriors in the inversion are fully identified, while the concentration toward small z-scores indicates that the true values ​​are precisely included in the posteriors. Therefore, the distribution in the lower right of the plot suggests an ideal Bayesian inversion. To further confirm the accuracy of the estimation in spatial excitability, based on the inferred z-scores of i ∈ {1, 2, …, 84}... The calculated confusion matrix is ​​in Figure 3D As shown in the diagram, the diagonal values ​​in the confusion matrix labeled HZ, PZ, and EZ indicate that the predefined class of all brain nodes was accurately predicted (accuracy = 1.0, misclassification = 0.0).

[0124] To investigate whether BVEP is a platform-independent framework, we also used PyMC3 to estimate spatial maps of epileptogenicity across different brain regions. For the two patients analyzed, the same precision was obtained by inverting Equation 4 in both Stan and PyMC3. These results indicate that BVEP inversion in Stan and PyMC3 yields similar estimates of spatial maps of epileptogenicity across brain regions in both analyzed patients.

[0125] In addition, NUTS-specific diagnostics were monitored to check whether the Markov chain had converged. The diagnostic plots showed that there was no divergent transition in the HMC indicating an effective probe of the posterior density. Furthermore, no NUTS iterations reached the maximum tree depth (the value for running NUTS is specified here as 10.0), indicating that the optimal number of frog-jump steps required for the Hamiltonian simulation is sufficiently lower than the maximum. These diagnostics collectively validate that the NUTS samples have converged to the target distribution.

[0126] To illustrate the mechanisms underlying seizure initiation and progression within the BVEP model, the phase plane topology characterizing the dynamics of different brain node types in the BVEP model is analyzed using simulated (apical) and predicted (basal) phase plane topologies. Figure 4 Presented in the diagram. In the plotted phase plane, the x and z zero-slope lines are colored dark gray, with the intersection of these lines marking the system's fixed points. From left to right, the columns correspond to brain nodes designated as HZ, PZ, and EZ, respectively. Full circles and empty circles indicate stable and unstable fixed points, respectively. Figure 4 A and Figure 4 D observes that the trajectory of HZ (node ​​number 1) is attracted to the system's stable fixed point (on the left branch of the cubic x-zero slope line), meaning that seizures are not triggered. For PZ (node ​​number 6), due to the coupling strength and the excitability value approaching the epileptogenic threshold, the z-zero slope line shifts downward, causing a bifurcation, thus allowing seizure propagation here (see...). Figure 4 B and Figure 4 E). For EZ (node ​​number 7), the system exhibits an unstable fixed point due to its high excitability. In this system, Epipuleptor has a limiting cycle and spontaneously triggered seizures (see...). Figure 4 C and Figure 4 F). Note that the topology of the simulated and predicted phase plane trajectories shows very good agreement, except for the state variable z. i Beyond the magnitude (from which the estimate indicates the result of a larger parameter recovery). Note that only the fast variable x 1,i The activity is used as the target for fitting the observed data.

[0127] To compare the BVEP inversion using the NUTS and ADVI schemes Figure 5The histogram of the MCMC samples is shown, along with the posterior kernel density estimate generated from NUTS (left plane) and the kernel density estimate obtained via ADVI (right plane). This plot shows that NUTS and ADVI perform similarly in their posterior estimates, except that the mean-field ADVI slightly underestimates the variance compared to the estimate obtained via the NUTS algorithm. However, the true values ​​of excitability (vertical dashed lines) in both methods, supported by the posterior density, indicate successful parameter recovery. Samples corresponding to brain nodes designated as HZ, PZ, and EZ are shown in gray, light gray, and dark gray, respectively. Note that the priors for all 84 brain regions included in the analysis are assumed to be a normal distribution centered at -2.5 with a standard deviation of 1.0, as shown in blue. To invert the BVEP model via the NUTS algorithm, 200 sampling iterations and 200 warm-ups were used with an expected acceptance probability of 0.95, while for running ADVI, the maximum number of iterations and convergence tolerance were set to 50,000 and 0.001, respectively. In terms of computation time, for these algorithm configurations, sampling via NUTS took 23993.5 seconds, while ADVI took 5392.62 seconds.

[0128] Once the model parameters have been estimated, the convergence of the MCMC samples needs to be evaluated. To test the reliability of the inferred estimates, we monitor the potential size reduction factor. This is because it is the most reliable quantitative measure of MCMC convergence. Additionally, posterior samples from the joint posterior probability distribution are plotted to demonstrate the efficiency of the transformed non-centered parameterization compared to the centered form of the parameterization. Figure 6 Top row indicates information from hyperparameters and The posterior samples of the joint posterior probability distribution are given, where the hyperparameters are the standard deviations of the process (dynamic) noise and the measurement noise, respectively (see Equation (4)). In this figure, the left and middle columns show the results of sampling by NUTS using parameterized non-centered and centered forms, respectively. For ease of comparison with NUTS, the last column shows the results from the mean-field variant of ADVI. The points in each scatter plot represent 200 samples extracted from the joint posterior probability distribution. Figure 6 A and Figure 6 As can be clearly seen in B, there is no correlation between the posterior samples extracted from the non-centered parameterization, while samples from the centered form exhibit high collinearity among the hyperparameters. This high collinearity leads to inefficient posterior probing, which results in a reduction in the number of effective samples and an increase in... The values ​​are observed in terms of quantity. All hidden states and parameters estimated through the non-centered form. The value is below 1.05 (see Figure 6D), while over 82% of estimates using the median form have a value higher than 1.1. Value (see) Figure 6 E). This indicates that the Markov chain converges to the parameterized non-centered form but not to the centered form. The effective number of samples returned by the centered form through NUTS is related to the number of iterations (N). eff =N iter The ratio of ) is less than 0.001 for all estimated parameters. This indicates poor sampling from the parameterized center form, as it generates a small number of individual samples per Markov chain.

[0129] Furthermore, the scatter plot extracted from the joint posterior probability distribution between the hyperparameters σ and σ' estimated by mean-field ADVI is shown in... Figure 6 As shown in C. Since, by definition, the mean-field variant of ADVI ignores cross-correlation between parameters, the samples extracted using mean-field ADVI do not exhibit correlation between hyperparameters. Finally, to check the convergence of ADVI, the lower bound of evidence (ELBO), i.e., the variational objective function versus the number of iterations, is plotted (see [reference]). Figure 6 F). Although the algorithm appears to converge in 10,000 iterations, it runs thousands of additional iterations to ensure convergence until the change in ELBO drops below a tolerance of 0.001.

[0130] Finally, this invention provides a probabilistic framework (i.e., a Bayesian virtual epilepsy patient (BVEP)) to infer a spatial map of epileptogenicity for the development of personalized, large-scale brain models of epilepsy spread (see [link]). Figure 1 The workflow for constructing the BVEP brain model consists of two main steps: In the first step, the VEP, i.e., a personalized large-scale brain network model of epileptiform spread, is constructed. In the VEP model, the dynamics of brain nodes are managed by a neural population model of epileptiform spread (i.e., Epileptor), a general model that realistically reproduces the onset, progression, and cancellation of seizure patterns across species and brain regions (Jirsa et al., 2014). Epileptor is coupled through the patient's connectome to combine a mean-field model of aberrant neuronal activity with subject-specific brain anatomy information obtained from non-invasive diffusion neuroimaging techniques (MRI, DTI). Along with the patient data, the VEP model is then equipped with a spatial map of epileptiform spread across different brain regions. In the second step, the VEP is embedded as a generative model in the PPL tool (Stan / PyMC3) to infer and validate the spatial map of epileptiform spread across different brain regions. Using PPL along with high-performance computing to run several MCMC chains in parallel enables system and effective parameter inference to fit and validate the BVEP model against patient data.

[0131] To demonstrate the potential functionality of BVEP in predicting seizure onset and spread, different spatial maps of epileptogenicity were used to simulate simple and complex seizure spread (see [link]). Figures 2A-2E These synthetic data were used for fitting because, given ground-real conditions with the model parameters, standard error metrics such as posterior shrinkage of the confusion matrix and posterior z-scores could be used to validate the accuracy of the estimation, thus evaluating the performance of the proposed approach. Results demonstrate that significant similarities between simulated and predicted seizure activity related to onset, propagation, and termination can be achieved on both synthetic datasets by inverting large-scale brain network models using PPL (Stan / PyMC3). While the simulations were generated from a full VEP model including five state variables for each brain node, a 2D-reduced variant of the model still successfully predicted key data features such as onset, propagation, and cancellation of seizure patterns, while significantly reducing the computational time of Bayesian inference. This 2D reduction was limited to modeling the average of rapid emissions during burst seizure states, as illustrated (see [link to documentation]). Figures 3A-3D and Figure 5 ( ) is a sufficient feature for correctly estimating the spatial map of epileptogenicity. The results indicate that the BVEP model can accurately estimate the spatial map of epileptogenicity across different brain regions (see Figures 3A-3D The true values ​​of excitability for all brain nodes included in the analysis are supported by estimated posterior density with 100% classification accuracy based on the confusion matrix. Furthermore, the distribution converges towards the set of small z-scores along with the set that contracts towards the large posterior (i.e., Figure 3C (The bottom right corner of the image) together confirms the reliability of the model inversion. Note that the accuracy obtained relying on the confusion matrix may be uncertain, as the accuracy of the estimate returned by this metric depends only on the mean of the estimated posterior density. For example, consider inference where the mean of the posterior is almost identical to the ground truth, but there is significant uncertainty in the estimate. In such cases, the confusion matrix can produce high accuracy performance, and plotting the posterior z-score and posterior shrinkage is particularly useful for identifying failures in inference (such as overfitting or improperly chosen priors that bias the estimate).

[0132] Understanding brain dynamics in epilepsy is crucial for developing therapeutic approaches to brain interventions to improve surgical outcomes. The complete classification of epileptic seizures has been extensively investigated elsewhere using theories of nonlinear dynamic systems, with thorough descriptions of the bifurcations that cause onset, cancellation, and seizure evolution characteristics (Jirsa et al., 2014). In the parametric space description of Epipuleptor, seizure onset and cancellation are described via saddle nodes and co-dipping bifurcations. The acute dynamic effects in the BVEP model depend primarily on the interactions between the network node model (Epileptor), patient-specific structural connectivity (from dMRI), and the spatial maps of epileptogenicity (EZ, PZ, HZ). Based on the dynamic properties of the Epipuleptor model, brain regions are classified into three main types: EZ (presenting unstable fixed points corresponding to brain regions responsible for seizure onset), PZ (approaching saddle node bifurcations corresponding to candidate brain regions responsible for seizure propagation), and HZ (presenting stable fixed points corresponding to healthy brain regions). This approach allows us to define the spatial map of epileptogenicity based on excitability parameter values ​​(which are the target values ​​for fitting).

[0133] It is important to note that excitability values ​​approaching the epileptogenic threshold do not guarantee the propagation of seizures originating from pathological brain regions (i.e., those responsible for seizure initiation designated as EZ) to such brain regions defined as PZ. Through detailed patient assessments, individual structural connectivity has been reported to be crucial for predicting seizure spatial propagation. However, recent studies have shown that purely structural information is insufficient to predict seizure propagation and eventual cessation. Abnormal activity in the recovery region is instead a complex network effect, dependent on the interaction of multiple factors, including epileptogenicity of brain regions (node ​​dynamics), individual structural connectivity (network structure) (Jirsa et al., 2017), and brain state dependence (network dynamics). Furthermore, nonlinear and multiple propagation patterns exist, which can be observed for the same set of excitability parameters due to coupled nonlinear system dynamics (see [link to relevant documentation]). Figures 2A-2E In this work, seizure remission is characterized by the complex spatiotemporal dynamics of large-scale brain networks, i.e., seizures originate from local networks and remit strongly coupled candidate brain regions to pathological regions by disrupting their stable dynamics (if K = 0, seizure remission does not exist). Among the candidate brain regions of seizure propagation, nodes PZ are more strongly coupled to pathological regions defined as EZ due to stronger connections. idx ={28} can be restored through weak global coupling. Stronger coupling is required for seizure restoration for all other candidate brain regions. This is consistent with experimental observations that seizures tend to have a common spatial origin in the same patient. Based on this knowledge, it imposes a weak informational prior on the ground-based global coupling parameters. Overestimation of the global coupling parameters leads to misclassification of the PZ as the HZ (see [link to relevant documentation]). Figure 2E and Figure 3BHowever, underestimating coupling can lead to misclassification of the PZ as the EZ. But stability analysis of network dynamics indicates that seizure propagation is controlled through optimal intervention of the structural connectivity matrix, meaning that patient-specific network connectivity can predict seizure propagation patterns. Therefore, seizure propagation may not be easily controlled by simple cutting of individual nodes, as it has been reported that in surgical treatment of epilepsy, resection does not necessarily induce postoperative seizures in the brain.

[0134] In this study, the analysis of the observed systems and predicted phase plane trajectories was performed across different brain regions to gain a better understanding of the mechanisms underlying the proposed pattern of seizure initiation and propagation (see [link to study]). Figure 4 For different brain node types (e.g., EZ, PZ, and HZ), the dynamics of seizure onset and recovery in the phase plane were fully captured through prediction. From a reasoning perspective, a good correspondence was observed with the phase diagrams of the observed systems, including equilibrium (intersection of zero-slope lines), the stability or instability of equilibrium, and the flow of trajectories. These results validate our Bayesian inversion process to understand the spatiotemporal evolution of seizure activity, paving the way for further research into possible seizure prevention mechanisms.

[0135] Both the NUTS and ADVI schemes were used to infer the spatial map of epileptogenicity in a personalized whole-brain model of epileptiform spread. Results from both inference schemes yielded similar estimates of the spatial map of epileptogenicity across brain regions, except that ADVI slightly underestimated the variance compared to the estimate made using the NUTS algorithm (see [link to relevant documentation]). Figure 5 Using the similarity indicator variational approximation between the inversions of the two schemes provides a suitable alternative for NUTS sampling in BVEP model inversion. Our results demonstrate a significant reduction in computational cost (4-5x faster for the algorithm configuration used) in inference performed by ADVI compared to NUTS, which may be important when applying the BVEP approach to large datasets of patient populations. While ADVI is generally known to be more computationally attractive than NUTS, using this approximation to discover algorithmic problems can be challenging. The convergence of ADVI can be assessed by monitoring the running average of ELBO changes, while NUTS is equipped with several general and specific diagnostics to assess whether the Markov chain has converged. Furthermore, ADVI may get stuck in local minima during gradient descent optimization, and its mean-field variant cannot cover all modes of multimodal posterior density.

[0136] Finally, we investigate the efficiency of the non-centered parameterization. Consistent with previous studies demonstrating NUTS's sensitivity to parameterization, our results indicate that the non-centered form of the parameterized inversion of the nonlinear state-space equations yields an efficient parameter space exploration, while the sampled centered form proves to be an inefficient exploration due to high collinearity among the model parameters (see [link to relevant documentation]). Figure 6A and Figure 6 D and Figure 6 B and Figure 6 E). Additionally, based on convergence diagnostics (e.g., (This proves that, compared to the parameterized centering form, samples generated by NUTS converge faster in the non-centering parameterization.)

[0137] A new approach to building personalized in-silico brain network models based on Bayesian inference within PPL tools such as Stan and PyMC3. While several PPL libraries have been developed for Bayesian inference, only a few have been built around efficient sampling algorithms, such as NUTS, which avoids random walk behavior and is sensitive to cross-correlation parameters. Both Stan and PyMC3 provide automatic differentiation for NUTS and ADVI to efficiently compute gradients without user intervention. Stan is a general-purpose and flexible package with interfaces to common data science languages ​​and also provides extensive diagnostics for MCMC convergence. PyMC3 provides several MCMC algorithms through direct model specification in native Python code. Our implementations in Stan and PyMC3 produce similar estimates of epileptogenic spatial maps across brain regions, indicating that BVEP is platform-independent. However, PyMC3 requires a larger number of warm-up iterations to achieve the same post-acceptance convergence as our implementation in Stan. This is attributed to the differences in the implementation of NUTS in Stan and PyMC3. A comparison of implementations in Stan, PyMC3, and other alternative PPL packages is beyond the scope of this paper.

[0138] This invention is the first personalized large-scale brain network modeling approach for inferring epileptogenic spatial maps (the nature of nodes) based on patient-specific whole-brain anatomical information (i.e., network structures derived from dMRI). Dynamic causal modeling (DCM) is a well-developed framework for analyzing neuroimaging modalities (such as fMRI, MEG, and EEG) through neuroquality models, enabling inferences related to couplings (effective connectivity) between brain regions to infer how changes in neuronal activity in one brain region are caused by activity in other regions through modulation in potential couplings. Using DCM, focal seizure activity in electrocorticography (ECoG) data has recently been studied to estimate key synaptic parameters or coupling connections using observed signals in human subjects. In another study, a Bayesian belief update scheme for DCM was used to estimate synaptic drivers of cortical dynamics during seizures from EEG / ECoG recordings at minimal computational cost. While DCM can be used to model and track changes in excitability (inhibitory balance during seizure onset / offset), these studies are based on a single neuromass model (i.e., modeling a small number of cortical sources) and represent nonlinear ordinary differential equations of the neuromass model approximated by their linearization, which can only model seizure onset or offset, not both. In this paper, the Bayesian Virtual Epilepsy Patient (BVEP) model is able to characterize the whole-brain spatiotemporal nonlinear dynamics of seizure propagation. This approach allows for the description of the onset and offset of burst states and the alternation between normal and burst cycles. The BVEP approach relies on patient-specific structural data rather than formulating the inverse problem entirely based on unknown model parameters used in DCM. It is also worth mentioning that the dynamics of the system are inferred by coupling fast and slow timescales (see Equation (3)), so that the changes in slow variables depend on the hidden states of fast activity, while assuming that only the activity of fast variables is observed. In this study, timescale separation in the Epiliptor model enables reliable capture of the complete evolution of complex dynamics, ranging from pre-burst to onset, burst evolution, and offset, rather than using time-varying parameters. Future extensions to the current work could explicitly examine the non-static dynamics of the network to investigate the conditions underlying the mechanisms of outbreak initiation: whether outbreak onset is more likely to occur through deterministic parameter changes such as bifurcation or due to jumps caused by noise-driven transitions between bistable attractors. The Bayesian inversion in the current work is based on autotuning algorithms (such as NUTS and ADVI) implemented using fast automatic differentiation with computational gradients. This allows for efficient sampling from complex and high-dimensional posterior distributions with interconnected parameters compared to conventional sampling algorithms. Several MCMC convergence diagnostics are also employed to enrich the inference within the provided framework to assess the reliability of the estimates.

[0139] Various non-invasive and invasive methods have been used to improve preoperative assessment for identifying EZ (extracorporeal lesions) and thus increase surgical success rates. The adoption of BVEP models in clinical treatment and brain interventions will require quantifying model outcomes in fitting patients' empirical assistive functional signals (such as EEG, MEG, SEEG, and fMRI signals). In this framework, it is straightforward to combine further knowledge from the preoperative assessment, such as MRI lesions and clinical assumptions about EZ. Since BVEP models can be considered a general approach to large-scale brain modeling, they offer a promising avenue for inference from clinically used non-invasive imaging signals (EEG, MEG, fMRI) and invasive measurements (such as SEEG signals). Results indicate that fitting to patients' empirical SEEG data (not shown) can be successfully achieved according to the proposed method of this invention. It should be noted that in the case of empirical SEEG recordings, source localization is an ill-conditioned problem due to the sparsity of the guiding field matrix, which can affect the accuracy of the estimation. In principle, it is possible to systematically test surgical strategies using BVEP models, but practical clinical application still awaits investigation and validation in future work.

[0140] In summary, this invention establishes a link between probabilistic modeling and personalized brain network modeling to systematically predict the location of seizure onset in virtual epilepsy patients. It demonstrates step-by-step how the proposed framework allows for the inference of epileptogenic spatial maps based on large-scale brain network models derived from non-invasive structural data of individual patients. The invention utilizes advanced efficient sampling algorithms that provide accurate and reliable estimates validated through posterior behavioral analysis and convergent diagnostics. In conclusion, the use of personalized brain network models, aided by PPL, provides sound guidance for the development of comprehensive clinical hypothesis testing and novel surgical interventions.

Claims

1. A computer-based method for inferring the epileptogenicity of brain regions that were not observed to be restored or not observed to be restored during seizure activity in an epileptic patient, comprising the following steps: Provides a computerized model for modeling the various regions of the primate brain and the connectivity between said regions; The computerized model is provided with a model capable of reproducing the dynamics of epileptic seizures in the primate brain, the model being a function of parameters of the epileptogenicity of regions of the brain; Structural data of the brain of the epilepsy patient is provided, and the structural data is used to personalize the computerized model in order to obtain a virtual epilepsy patient VEP brain model; The state-space representation of the virtual epilepsy patient VEP brain model is converted into the probabilistic programming language PPL using probabilistic state transitions to obtain a probabilistic virtual epilepsy patient brain model BVEP. The state-space representation of the virtual epilepsy patient VEP model is incorporated into the probabilistic virtual epilepsy patient brain model BVEP as state transition probabilities. The method provides acquired electroencephalogram (EEG) or magnetoencephalogram (MEG) data of the brain of the epilepsy patient, and uses this data to fit a probabilistic virtual brain model of the epilepsy patient in order to infer the epileptogenicity of brain regions that were not observed to be rejuvenated or not observed to be rejuvenated during seizure activity in the epilepsy patient's brain. Specifically, the probabilistic virtual epilepsy patient brain model is generated based on a generative model derived from the state space representation of the virtual epilepsy patient. The state space representation of the virtual epilepsy patient is determined systematically based on a nonlinear stochastic differential equation of the following form: in, It is an n-dimensional vector of the system state that evolves over time. It is the initial state vector at time t = 0. Includes all unknown parameters of the virtual epilepsy patient model, where u(t) represents the external input. Let v(t) represent the measurement data subjected to measurement error, f be a vector function describing the dynamic properties of the system, and h be the measurement function of the guiding field matrix, where w(t) represents the process noise.

2. The method according to claim 1, wherein, The probabilistic programming language is a Bayesian programming language, the probabilistic virtual epilepsy patient brain model is a Bayesian virtual epilepsy patient brain model, and Bayesian inference is used to infer the epileptogenicity of the brain regions that were not observed as restored or not observed as unrestored.

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

4. The method according to claim 1, wherein, The model capable of reproducing the dynamics of the epileptic seizures in the primate brain is a model that reproduces the dynamics of the onset, progression, and cancellation of seizure events. It includes state variables coupled to two oscillating dynamic systems at three different timescales: the fastest timescale, where the state variables record the rapid discharges during the burst state; the intermediate timescale, where the state variables represent slow spikes and wave oscillations; and the slowest timescale, where the state variables are responsible for the transition between the inter-burst period and the burst state, and where the degree of epileptogenicity of brain regions is represented by the value of an excitability parameter.

5. The method according to claim 1, wherein, To obtain the probabilistic virtual brain model of the epilepsy patient, a spatial map of the epileptogenicity of the epilepsy patient's brain is provided. The epileptogenic spatial map classifies the brain regions of the epilepsy patient's brain into epileptogenic zones (EZ) that can spontaneously trigger seizures, propagation zones (PZ) that do not spontaneously trigger seizures but can recover during seizure evolution, and healthy zones (HZ) that do not spontaneously trigger seizures.

6. The method according to claim 1, wherein, The state transition probability is: in, This represents the probability of transition from state x(t) to x(t + dt).

7. The method of claim 1, wherein the generative model is defined based on likelihood and prior model parameters, the product of the likelihood and prior model parameters producing the following joint density: in, Prior distribution This includes prior beliefs about the values ​​of hidden variables and latent parameters, while the conditional likelihood term... This represents the probability of obtaining an observation using a given set of parameter values.

8. The method of claim 1, further comprising implementing at least one sampling algorithm to infer the epileptogenicity of brain regions in the patient’s brain that were not observed to be restored or not observed to be unrestored during the seizure activity.

9. The method according to claim 8, wherein, The sampling algorithm is either Markov chain Monte Carlo or variational inference algorithm.

Images