Method for determining onset time and cerebral regional excitability

Abstract

The present invention relates to a method for determining the onset time and excitability of brain regions that are not observed as recruited or non-recruited during seizure activity in the brain of an epilepsy patient. The method according to the present invention comprises the steps of: providing a dynamic model of the propagation of an epileptic seizure in a brain network; providing a statistical model that defines the probability of the dynamic model generating a set of observations of the state of the brain network; using the statistical model and a data set of observations of a training cohort to train the dynamic model of the propagation of the epileptic seizure; and inverting the trained dynamic model and using the statistical model to infer the onset time and excitability of a third region from the onset times observed for a first and a second region.

Description

Technical Field

[0001] The present invention relates to a method for determining the excitability and onset time of brain regions that are not observed as recruited or not recruited during seizure activity in the brain of an epilepsy patient. Background Art

[0002] A possible treatment for medically refractory epilepsy patients is surgical intervention aimed at removing one or more suspected epileptogenic zones, i.e., brain regions that are the cause of seizures. However, the success rate of these surgical interventions is only 60 - 70%.

[0003] However, the success rate can be improved using the method disclosed in document WO2018 / 015778A1. This document discloses a method for modulating epileptogenicity in the brain of an epilepsy patient, comprising the steps of: providing a virtual brain; providing a model of an epileptogenic zone and a propagation zone and loading the model into the virtual brain to create a virtual epileptic brain; acquiring data of the brain of the epilepsy patient; identifying the location of at least one possible epileptogenic zone in the data; fitting the virtual epileptic brain to the data collected from the epilepsy patient and parameterizing the at least one possible epileptogenic zone in the virtual epileptic brain as an epileptogenic zone; and simulating the effect of network modulation of a clinical intervention on the patient's brain within the virtual epileptic brain. Simulating a clinical intervention on the patient's brain can then allow for the development of improved surgical strategies and improve the low success rate of surgical interventions.

[0004] In fact, the low surgical success rate is mostly attributed to the failure to localize the epileptogenic zone, and part of the reason for the failure to localize the epileptogenic zone is the incomplete image of the whole-brain seizure propagation pattern due to insufficient spatial sampling. In fact, the current standard for preoperative evaluation is to use implanted depth electrodes (stereoelectroencephalogram, SEEG) or subdural electrode grids. None of these methods allows exploration of the entire brain, and they are usually limited to those regions that are suspected to be part of the epileptogenic network based on non-invasive evaluation.

[0005] Computer-aided methods aimed at improving surgical planning avoid this situation of insufficient spatial sampling, assuming that the relevant targets that may be resected are located in the regions explored by electrode implantation. In this case, the behavior of other parts of the brain network is not taken into account. Some of these methods are based on the analysis of recorded signals using spectral or temporal features. Other methods among these are based on the analysis of functional networks derived from intracranial recordings.

[0006] Known methods for modeling the activity in the entire brain network (not just the explored subnetwork) require some manual adjustment of the model settings, usually based on the epileptogenic zone hypothesis specified by a clinical expert, and thus cannot be fully automated. SUMMARY OF THE INVENTION

[0007] According to a first aspect, the present invention relates to a method for determining the excitability and / or onset time of brain regions that are not observed as recruited or not observed as non-recruited during seizure activity in the brain of an epileptic patient, comprising the steps of:

[0008] Providing a computerized brain network that models the individual regions of the patient's brain and the connectivity between the regions for a training cohort of the epileptic patient's brain;

[0009] Providing, for the training cohort and for the patient's brain, a data set of observations of the state of the brain network during epileptic seizures, the observations defining regions in the brain network as

[0010] a first region that is observed as recruited during epileptic seizure activity in the patient's brain at the onset time,

[0011] a second region that is observed as non-recruited during the epileptic seizure activity of the patient's brain, and

[0012] a third region that is not observed as recruited or not observed as non-recruited during the seizure activity of the patient's brain;

[0013] Providing a dynamic model of the propagation of epileptic seizures in the brain network;

[0014] Providing a statistical model that defines the probability that the dynamic model generates the set of observations of the state of the brain network;

[0015] Using the statistical model and the data set of observations of the training cohort to train the dynamic model of the propagation of epileptic seizures; and

[0016] Inverting the trained dynamic model and using the statistical model to infer the onset time and / or excitability of the third region from the onset times observed for the first and second regions.

[0017] Preferably, - wherein, the initial movement time and excitability of the third region are inferred; - a computerized brain network is obtained from magnetic resonance neuroimaging and / or diffusion-weighted magnetic resonance imaging data; - a dataset of observations of the state of the brain network during a seizure is obtained by running a seizure onset detection algorithm on intracranial electroencephalogram signals and mapping these detected onset times to regions of the brain network; - the seizure onset detection algorithm employs time-frequency analysis of intracranial electroencephalogram signals; - mapping the onset times detected in intracranial electroencephalogram signals to brain regions is based on the physical distance between the electrode contacts from which the signals are recorded and the brain regions; - the dynamic model describes the evolution of the slow variables of individual regions of the network through an activation function, which is a function of the excitability of the region and the network effect, and wherein the regional onset time is defined as the time when the slow variable crosses a given threshold; - training of the dynamic model is equivalent to finding the optimal set of parameters of the activation function; the Hamiltonian Monte Carlo method is used to infer the parameters of the activation function, regional excitability, and onset time of the dynamic model; - for a brain network with n regions, the dynamic model defining the propagation of a seizure in the brain network is given by:

[0018] For i = 1, …, n

[0019] z i (0) = 0

[0020] where the function f q : R x [0; 1] → R + is the activation function, c i is the node excitability, and W = (ω ij ) is the connectivity matrix, which is normalized such that max i Σ j ω ij = 1, and H is the Heaviside step function representing the transition from the healthy state to the seizure state when the slow variable crosses the threshold z = 1; - the statistical model is a hierarchical model based on the principles of Bayesian inference, where the top-level parameters are the parameters of the activation function and the bottom-level parameters are the regional excitabilities; - the statistical model includes the following assumptions: the excitabilities of all regions for all training data have the same prior distribution; - the function f is a bilinear function followed by a power operation, and the parameters of the function are the values of the function f at four specified points in R x [0; 1]; - the statistical model such as q 11 , q 12 ~Normal(0, σ q )

[0021]

[0022] For k = 1, …, n seizures :

[0023] c k ~Normal(0, 1)

[0024] t k =P(c k , q, W k )

[0025]

[0026] t lim ~Normal(min(t k , nonseizing, t lim ), σ t )

[0027] wherein is the initial motion time of the observed seizure node, σ q = 30, σ t = 5s, t lim = 90s, and P(c k , q, W k ) represents a dynamic model that maps the excitability c k , the parameter q, and the connectome matrix W k to the initial motion time t k . BRIEF DESCRIPTION OF THE DRAWINGS

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

[0029] Figure 1 is a schematic overview of the seizure propagation inference problem according to the present invention;

[0030] Figure 2 shows the training and application phases of the method according to the present invention; and

[0031] Figure 3 shows the inference results obtained by the method according to the present invention for a single seizure with strong coupling and 21 observation regions;

[0032] Figure 4 shows the probability of seizure occurrence in the hidden nodes obtained from an embodiment of the method according to the present invention;

[0033] Figure 5 shows the probability that the inferred initial motion time is within T seconds of the true initial motion time in the hidden seizure region according to the method of the present invention;

[0034] Figure 6A , Figure 6B and Figure 6CShows the true and inferred regional excitability for the observed region according to the method of the present invention; and

[0035] Figure 7A 、 Figure 7B and Figure 7C Shows the true and inferred node excitability for the hidden region according to the method of the present invention. Detailed implementation

[0036] The present invention relates to a method for determining the excitability and / or onset time of brain regions that are not observed as recruited or non-recruited during seizure activity in the brain of an epilepsy patient. It relates to a method for attempting to infer the activity of the entire brain network during seizures from the activity and structural connectome of the observed brain subnetworks. The method is based on the assumption that the activity in any brain region can be classified as a normal state or a seizure state, and the change time (regional onset time) between these states is the only relevant feature for inferring the propagation pattern. The method is divided into two phases. In the training phase, the global parameters of the model are inferred from the training data, i.e., those parameters shared between different subjects and seizures. In the application phase, the trained model is applied to single seizure data, and the activity in the entire brain network is inferred for the seizure. Due to this data-driven method, only a few clearly interpretable constants need to be specified, so no manual adjustment is required at the global or patient-specific level.

[0037] As Figure 1 shown, the state of the brain network during the evolution of an epileptic seizure is only partially known. Some brain regions are observed to participate in seizure activity, with an onset time of t i . Other brain regions are observed, but they do not participate in seizure activity. There are also some regions that are not observed but are hidden, i.e., it is not known whether they participate in seizure activity. These are the brain regions that are not observed as recruited or non-recruited during seizure activity in the brain of an epilepsy patient.

[0038] As Figure 2 shown, the method of the present invention allows the use of the structural connectivity known from diffusion-weighted imaging to infer the state of hidden nodes. (Bottom) The two phases of inference. In the training phase, the hyperparameters of the model are learned from the training data. In the application phase, the model is applied to a single seizure.

[0039] To implement the method of the present invention, computerized brain networks are provided. These networks are for a training cohort of epileptic patients' brains and model the connectivity between individual regions of the patient's brain and between said regions. The computerized brain networks are obtained from magnetic resonance neuroimaging and / or diffusion-weighted magnetic resonance imaging data. In addition, a dataset of observations of the state of the brain network during seizures is provided for the training cohort and for the patient's brain, the observations defining regions in the brain network as: a first region that is observed to be recruited in the seizure activity of the patient's brain at the onset time; a second region that is observed to be non-recruited in the seizure activity of the patient's brain; and a third region that is not observed to be recruited or non-recruited in the seizure activity of the patient's brain.

[0040] These form the input data, thus constituting the connectome matrix. The computerized brain networks are obtained from magnetic resonance neuroimaging and / or diffusion-weighted magnetic resonance imaging data. For example, a single structural connectome matrix is obtained from the patient's T1 and diffusion-weighted MRI scans using existing known available software. Methods for connectome reconstruction are well known in the art.

[0041] These input data also include regional onset times. These are the information on whether and when a region enters the seizure state. To obtain such information data, the following steps are performed:

[0042] - First, intracranial seizure recordings are provided, and in each bipolar channel, the onset time is detected by evaluating when the signal power crosses a given threshold. Specifically, the whitened power in the [4, 60] Hz band and a log ratio threshold equal to 2 relative to the pre-seizure baseline are used. However, it should be noted that other embodiments are possible, either using different frequency bands or other methods, such as machine learning-based detection methods.

[0043] - Second, all bipolar channels are assigned to brain regions. For each channel, the nearest brain region is selected, determined by the volume segmentation provided by the known FreeSurfer TM software and the electrode contact positions determined by the brain CT scan of the implanted electrodes. If there are multiple regions with distances less than <2 near the contact, these regions are not assigned, and the information from the channel is not used.

[0044] - Third, regions without an assigned channel are marked as hidden. For a region with one assigned channel, the detected channel onset time is used as the regional onset time. If multiple channels are assigned to a region, the median of the detected onset times is taken.

[0045] According to the present invention, a dynamic model of seizure propagation and recruitment in a brain network is provided. The dynamic model describes the evolution of the slow variable of a single region of the network through an activation function, which is a function of the excitability of the region and the network effect, and wherein the onset time of the region is defined as the time when the slow variable crosses a given threshold.

[0046] For a network with n regions, the model is defined by Equation (1) as follows:

[0047]

[0048] In this equation, z i is the slow variable of region i, similar to the slow variable in the epilepsy model (Jirsa, v., Stacey, w., Quilichini, p., Ivanov, a., Bernard, c., July 2014, On the nature of seizure dynamics, Brain 137(8), 2110 - 2113). The function f q :R x[0;1]→R + is the activation function, parameterized by the parameter vector q. The function f q is an increasing function with respect to its first parameter. The parameter c i is the node excitability, and W=(ω ij ) is the connectivity matrix, which is normalized such that max i Σ j ω ij =1. Finally, H is the Heaviside step function representing the transition from the healthy state to the seizure state when the slow variable crosses the threshold z = 1. The onset time of region i is defined as t i =min{t|z i (t)≥1}. The system is completed with the initial condition z i (0)=0.

[0049] Since the function f q is positive, this model implies that each region will start to seize at some finite time. In fact, this is not the case, so a seizure time limit is introduced in the statistical model, and all regions with an onset time greater than this limit are considered non-seizure regions.

[0050] For known vectors of the excitability c, the parameter vector q, and the connectome matrix W, the model uniquely defines the vector of the onset time t. The following shorthand is used to represent this mapping, corresponding to Equation (2) below:

[0051] P(c,q,W)=t (2).

[0052] For the activation function f qParameterization of, using an exponential bilinear function as f q (c,y) parameterization. From the interpolation points [c a ,y a , [c a ,y b , [c b ,y a , [c b ,y b (where c a = -1, c b = 1, y a = 0, y b = 1), the four coefficients q = (q 11 ,q 12 ,q 21 ,q 22 ) are used for parameterization, and the function is given by

[0053]

[0054] To ensure that f increases with c, additional restrictions are placed on the coefficients: q 21 > q 11 and q 22 > q 12 .

[0055] The method according to the present invention further comprises the following steps: providing a statistical model that defines the probability of the set of observations of the state of the brain network generated by the dynamic model; and using the dataset of observations of the training cohort and the statistical model to train the dynamic model of the propagation of epileptic seizures. This is the training phase. Training the dynamic model is equivalent to finding the optimal set of parameters of the activation function.

[0056] The statistical model is, for example, a hierarchical model based on the principle of Bayesian inference, where the top-level parameters are the parameters of the activation function and the bottom-level parameters are the regional excitabilities. However, this is just an example. The statistical model preferably includes the following assumption: the excitabilities of all regions for all training data have the same prior distribution.

[0057] In the training phase, data from multiple seizures are used to infer the optimal value of the parameter vector q, which is shared across all seizures. The statistical model used for model training is established as follows: assume that for each seizure, the values of the excitability vector C k follow a standard normal distribution, and the activation function parameter q is regarded as a hyperparameter of the model. In addition, assume that the observed onset time is measured imprecisely with a standard deviation σ t . Finally, to handle regions that do not have seizures, a limit t lim, and each region that experiences an episode after this limit is considered to have no episode. The complete multiple-episode statistical model is as follows:

[0058] Input data: Connectivity matrix W k , episode set and no-episode node set, onset times of episode nodes

[0059] Parameters: σ q = 30, σ t = 5s, t lim = 90s

[0060] Model:

[0061] q 11 , q 12 ~ Normal(0, σ q )

[0062] q* 21 , q* 22 ~ HalfNormal(0, σ q )

[0063] q = (q 11 , q 12 , q 11 + q * 21 , q 12 + q * 22 )

[0064] For k = 1, …, n seizures :

[0065]

[0066] This statistical model defines the posterior probability distribution for the hyperparameter q and the episode-specific excitability parameter c k . Using the sampler implemented in Stan TM software, samples are drawn from this probability distribution using the Hamiltonian Monte Carlo method. Since at this stage we want to obtain a point estimate of the hyperparameter q, the median of all samples is taken for each component of the parameter vector.

[0067] The method according to the present invention further comprises the following steps: inverting the trained dynamic model and using the statistical model to infer the onset time and / or excitability of a third region - in fact, at least one, or even all, of the third regions - from the onset times observed for the first and second regions. This is the application stage. At this stage, the parameters of the activation function of the dynamic model and the region excitability and onset time are inferred by the Hamiltonian Monte Carlo method.

[0068] Regarding the application phase, once the parameters of the activation function f are learned from the training batches of seizure data, the following single-seizure model can be used to process additional seizures one by one: q Once the parameters of the activation function f are learned from the training batches of seizure data, the following single-seizure model can be used to process additional seizures one by one:

[0069] Input data: Connectivity matrix W, seizure set and non-seizure node set, onset time of seizure nodes Parameters of the activation function q

[0070] Parameters: σ q = 30, σ t = 5s, t lim = 90s

[0071] Model:

[0072]

[0073] This model is a simplification of the multiple-seizure model using the knowledge of the parameter vector q. Like the multiple-seizure model, samples from the posterior distribution of the single-seizure model can be drawn using the Hamiltonian Monte Carlo method. The results of single-seizure inference are the sampled excitability c and onset time t.

[0074] In this section, results are presented showing the performance of the method on synthetic data generated by the same model. These results obtained using the method according to the present invention demonstrate that if the assumptions of the method are met in a real epileptic network, it may be possible to infer the hidden state. In fact, the present invention allows predicting whether a hidden region will have a seizure. Furthermore, if a hidden region has a seizure, the present invention allows accurately predicting when it will occur. Finally, the present invention allows recovering the regional excitability.

[0075] Example 1: Test data

[0076] For testing purposes, synthetic data was generated by the same model used for inference. Connectivity matrices from 10 subjects were used, with 84 nodes using the Desikan-Killiany segmentation. Two sets of 10 seizures each were generated from these 10 connectivities, all seizures having three different sets of activation function parameters: one for no coupling (q = (-5.0, -5.0, -3.0, -3.0)), one for weak coupling (q = (-6.5, -3.0, -3.5, 12.0)), and one for strong coupling (q = (-11.2, 5.3, -6.2, 69.3)). Finally, the number of observed regions among the 84 regions in the segmentation was set to 21, 42, and 63. A total of (2 sets) x (3 coupling strengths) x (3 numbers of observed regions) = 18 sets, with 10 seizures per set. For all cases, the excitability c was randomly drawn from a standard normal distribution, and the observed nodes were randomly selected. Inference was performed separately for each of the 18 sets.

[0077] Example 2: Results

[0078] Figure 3 Partially shows exemplary results of the inference for a single seizure, while Figure 4 and 5 Figures 6A, 6B, 6C, 7A, 7B, and 7C show the cumulative results from all tests.

[0079] In Figure 3 , for a single seizure with strong coupling and 21 observation regions, an example of the inference results is partially shown. The left partial panel shows the true (dots) and inferred (plots) excitabilities of 84 brain regions. Solid dots mark the observation regions, while hollow dots mark the hidden regions. The light gray regions are those of the seizure, while the dark gray regions are those without seizure. The numbers in the left column are diagnostic values. The numbers in the right column are the inferred probability p(c>1), which is equivalent to the probability that the node is an epileptogenic node. The right panel shows the true (light gray / dark gray dots) and inferred (gray dots) onset times of the regions. The rules for solid dots / hollow dots and light gray dots / dark gray dots are the same as in the left panel. The results show that in some cases, the onset times of the hidden nodes are correctly inferred.

[0080] Figure 4 Shows how well the hidden regions can be predicted to have a seizure or not. It shows the inferred probability of seizure activity in the hidden nodes. In the grid, the rows represent three different coupling strengths, and the columns represent different numbers of observation regions. Each panel in the grid represents the results of two sets of inferences, with ten subjects in each set. The histograms calculate the probability that the hidden regions are involved in seizure activity. The light gray is the histogram of the seizure nodes, and the dark gray is the histogram of the non-seizure nodes. The inset text shows the classification metric when regions are classified as having a seizure or not based on the inferred probability with a threshold equal to 0.5. The results show that only in the case of strong coupling can the seizure nodes and non-seizure nodes be reliably distinguished. The number of observed nodes (nobs) does not necessarily improve the classification quality.

[0081] Figure 5 Further shows how well the accurate onset time can be estimated for the seizure regions. It shows the probability that the inferred onset time falls within T seconds of the true onset time of the hidden seizure regions. This probability is calculated for each hidden seizure region by counting the fraction of onset time samples in the window [t true –T; t true +T]. The results show that in the case of strong coupling, the onset time can be reliably inferred, and increasing the number of observed nodes further improves the accuracy. On the other hand, in the case of no coupling, the number of observed nodes does not change the prediction quality.

[0082] Figure 6A , Figure 6B and Figure 6C show how regional excitability can be recovered for unobserved nodes. They show the true and inferred regional excitability for the observed regions. In the grid, the rows represent three different coupling strengths and the columns represent different numbers of observed regions. Each panel shows the true (x-axis) and mean inferred excitability (y-axis) for the seizure (light grey) and non-seizure (dark grey) regions. Perfect inference would be all points lying along the diagonal. The inset numbers represent the R 2 coefficient of the linear regression. In the absence of coupling, the inference is very good and when strong coupling is introduced, some accuracy is lost due to network effects.

[0083] Figure 7A , Figure 7B and Figure 7C show the same for the hidden nodes. They show the true and inferred node excitability for the hidden regions. The layout is the same as Figure 6A , 6B and 6C. In the absence of coupling, it is not possible to infer the excitability of the hidden regions (row "None"). For strong coupling, there is a relationship between the true and inferred excitabilities (row "Strong") which increases with the number of observed regions.

[0084] Thus, according to the present invention: using the strong coupling effect, it is possible to infer whether a hidden node is in seizure, as Figure 4 shown. Using the strong coupling effect, it is also possible to predict the onset time of node seizure, as Figure 5 shown. For the observed nodes, the excitability can be inferred well, especially in the case of weak coupling effects, as Figures 6A to 6C shown. In the case of strong coupling and a large number of observed nodes, the excitability of the hidden nodes is inferred to some extent, as Figures 7A to 7C shown.

[0085] Finally, using partial observations of seizure evolution in the brain network and combining with the hypothesis of seizure propagation along structural connections (especially estimable from diffusion-weighted MRI), it can be inferred whether and when hidden regions are recruited during seizures. For this purpose, a data-driven dynamic model of seizure recruitment and propagation on a weighted network is used. The strong non-linearity caused by the discontinuous change of node state from normal to seizure state enriches this simple dynamic model. For the inversion of this model, a Bayesian inference framework is preferably adopted, and for example, the Stan software is used to infer the model parameters by the Hamiltonian Monte Carlo method. The results of computational experiments using synthetic data generated by the same model show that the quality of the inference results depends on the number of observed nodes and the connection strength. The accuracy of the inferred onset time improves with the number of observed nodes and the strength of the network effect, while in all cases studied, the ability of the model to infer the excitability of hidden nodes is limited. In real seizures, these results suggest that the state of hidden nodes during seizures is inferred from incomplete observations.

Claims

1. A method for determining the excitability and / or onset time of brain regions that are not observed as recruited or not observed as non-recruited during seizure activity in the brain of an epilepsy patient, comprising the following steps: Provide a computerized brain network that models the various regions of the brain and the connectivity between the regions for a training cohort of epileptic patients' brains; Provide a dataset of observations of the state of the brain network during seizures for the training cohort; Provide a dynamic model of the propagation of seizures in the brain network, the dynamic model describing the evolution of the slow variables of individual regions of the brain network through a parameterized activation function that is a function of the excitability of the region, and wherein the regional onset time is defined as the time when the slow variable crosses a given threshold; Provide a statistical model that defines the probability that the dynamic model generates the set of observations of the state of the brain network; Use the statistical model and the dataset of observations of the training cohort to train the dynamic model of seizure propagation in order to determine the optimal set of parameters for the activation function; Provide a dataset of observations of the state of the brain network during seizures for a patient's brain, the observations defining regions in the brain network as a first region that is observed to be recruited in the seizure activity of the patient's brain at the onset time, a second region that is observed to be non-recruited in the seizure activity of the patient's brain, and a third region that is not observed to be recruited or non-recruited in the seizure activity of the patient's brain; Invert the trained dynamic model and use the statistical model to infer the onset time and / or excitability of the third region from the onset times observed for the first and second regions; wherein the dataset of observations of the state of the brain network during seizures is obtained by running a seizure onset detection algorithm on intracranial electroencephalogram signals and by mapping these detected onset times to the regions of the brain network.

2. The method according to claim 1, wherein, Infer the onset time and excitability of the third region.

3. The method according to claim 1, wherein, Obtain the computerized brain network from magnetic resonance neuroimaging and / or diffusion-weighted magnetic resonance imaging data.

4. The method according to claim 1, wherein, The seizure onset detection algorithm employs time-frequency analysis of intracranial electroencephalogram signals.

5. The method according to any one of claims 1 to 4, wherein, Mapping the onset times detected in the intracranial electroencephalogram signals to the brain regions is based on the physical distance between the electrode contacts from which the signals are recorded and the brain regions.

6. The method according to any one of claims 1 to 4, wherein, Use the Hamiltonian Monte Carlo method to infer the parameters of the activation function of the dynamic model and the regional excitability and onset time.

7. The method according to any one of claims 1 to 4, wherein, The dynamic model is a data-driven dynamic model.

8. The method according to any one of claims 1 to 4, wherein, For a brain network with n regions, the dynamic model of the propagation of seizures in the brain network is defined by the following formula: For i = 1, …, n z i (0)＝0 where z i is the slow variable of region i, the function f q : R x [0; 1] → R + is the activation function, c i is the node excitability, and W = (ω ij ) is the connectivity matrix, which is normalized such that max i Σ j ω ij = 1, and H is the Heaviside step function representing the transition from the healthy state to the seizure state when the slow variable crosses the threshold z = 1.

9. The method according to any one of claims 1 to 4, wherein, The statistical model is a hierarchical model based on the principle of Bayesian inference, wherein the top-level parameters are the parameters of the activation function and the bottom-level parameters are the regional excitabilities.

10. The method according to any one of claims 1 to 4, wherein, The statistical model includes the assumption that the excitabilities of all regions for all training data have the same prior distribution.

11. The method according to any one of claims 1 to 4, wherein, The function f is a bilinear function followed by a power operation, and the parameters of the function are the values of the function f at four specified points in R x [0; 1].

12. The method according to any one of claims 1 to 4, wherein, The statistical model is q 11 ,q 12 ~Normal(0,σ q ) For k = 1, …, n seizures : c k ~Normal(0,1) t k = P(c k , q, W k ) t lim ~Normal(min(t k,nonseizing ,t lim ),σ t ) Among them, is the observed onset node initial motion time, σ q = 30, σ t = 5s, t lim = 90s, and P(c k , q, W k ) represents a dynamic model that maps the excitability c k , the parameter q, and the connectome matrix W k to the initial motion time t k .

Images