Methods for determining the time of origin and excitability of brain regions

Abstract

The present invention relates to a method for determining the onset time and excitability of brain regions that are not observed in the brains of epileptic patients, whether or not they experience seizure activity. The method according to the present invention comprises the steps of: providing a mechanical model of epileptic seizure propagation in a brain network; providing a statistical model that defines the probability of producing a set of observations of brain network states by the mechanical model; training the mechanical model of epileptic seizure propagation using the statistical model and a dataset of observations of a training cohort; and inverting the trained mechanical model and estimating the onset time and excitability of a third region from onset times observed in first and second regions using the statistical model.

Description

[Technical Field]

[0001] FIELD OF THE INVENTION The present invention relates to a method for determining the onset time and excitability of brain regions not observed in the brains of epileptic patients, both with and without the onset of seizure activity. [Background technology]

[0002] Background of the Invention A possible treatment for patients with intractable epilepsy is surgical intervention aimed at removing one or more suspected epileptogenic zones, i.e., the brain regions responsible for seizure initiation. However, the success rate of these surgical interventions is only 60-70%.

[0003] However, such success rates may be increased by utilizing the method disclosed in document WO2018 / 015778A1, which discloses a method for modulating epileptogenesis in the brain of an epileptic patient, comprising the steps of: providing a virtual brain; providing models of epileptogenic and propagation regions and loading the models into the virtual brain to create a virtual epileptic brain; acquiring brain data of an epileptic patient; identifying in the data the location of at least one possible epileptogenic region; fitting the virtual epileptic brain to the data acquired from the epileptic patient to parameterize the at least one possible epileptogenic region in the virtual epileptic brain as an epileptic region; and simulating in the virtual epileptic brain the effects of network modulation that mimic clinical interventions on the patient's brain. Mimicking clinical interventions on the patient's brain in turn allows for the definition of improved strategies for surgery and increases the low success rates of surgical interventions.

[0004] In fact, the low success rate of surgical procedures is due in most cases to the inability to localize the epileptogenic zone, which in part results from incomplete imaging of seizure propagation patterns throughout the brain resulting from insufficient spatial sampling. In practice, the current standard for preoperative evaluation uses either implanted depth electrodes (stereotactic intracranial electroencephalography, SEEG) or subdural grid electrodes. Neither of these methods allows for the exploration of the entire brain, and they are typically limited to regions suspected of being part of the epileptogenic network based on noninvasive assessment.

[0005] Computer-assisted methods aimed at improving surgical planning eliminate this insufficient spatial sampling by assuming that the relevant potential ablation targets lie within the region explored by the electrode implantation, without taking into account the behavior of the rest of the brain network. Some of these methods are based on the analysis of recorded signals using spectral or temporal features. Others are based on the analysis of functional networks derived from intracranial recordings.

[0006] Known methods that model the activity of whole-brain networks and not the activity of the examined sub-networks require some manual adjustment of the model settings, which are often hypotheses made by clinical experts about the epileptogenic zone and therefore cannot be used in a fully automated manner. Summary of the Invention

[0007] Summary of the Invention According to a first aspect, the present invention provides a method for manufacturing a semiconductor device comprising: providing a computerized brain network that models various brain regions and the connectivity between said regions for a training cohort of brains of epilepsy patients and for the patient brains; providing a dataset of observations of brain network states during epileptic seizures for the brains of a training cohort and patients, the observations classifying regions in the brain network as: At onset time, the first area observed to manifest epileptic seizure activity in the patient's brain; a second region of the patient's brain observed to be free of said epileptic seizure activity; and a third region of the patient's brain that is not observed with or without the onset of said seizure activity; providing a dataset, providing a dynamical model of epileptic seizure propagation within brain networks; providing a statistical model that defines the probability of producing the set of brain network state observations by the dynamical model; training a dynamical model of epileptic seizure propagation using the statistical model and a dataset of observations of the training cohort; inverting the trained mechanical model and utilizing a statistical model to estimate the onset time and / or excitability of a third region from the onset times observed in the first and second regions; The present invention relates to a method for determining the onset time and / or excitability of brain regions not observed in the brain of an epileptic patient, whether or not seizure activity occurs, comprising:

[0008] Preferentially, - wherein the onset time and excitability of a third region are estimated; - a computerized brain network is obtained from magnetic resonance neuroimaging and / or diffusion-weighted magnetic resonance imaging data; - a dataset of observations of brain network states during epileptic seizures is obtained by running a seizure onset detection algorithm on intracranial EEG signals and mapping these detected onset times to regions of the brain network; - the seizure onset detection algorithm utilizes a time-frequency analysis of the intracranial EEG signals; - the mapping of onset times detected in the intracranial EEG signals to brain regions is based on the physical distance between the electrode contacts recording the signal and the brain region; - a dynamical model describes the spread of a single region of the network to a slow variable by an activation function which is a function of the excitability of the region and the network effect, the onset time of the region being defined as the time when the slow variable exceeds a given threshold; - training the dynamical model amounts to finding an optimal set of parameters of the activation function; the parameters of the activation function of the dynamical model and the region excitability and onset times are estimated by a Hamiltonian Monte Carlo method; The dynamical model of epileptic seizure propagation in a brain network is defined by the following equation for a brain network with n regions:

number

number

number

number

number

[0009] Other features and aspects of the present invention will become apparent from the following description and accompanying drawings. [Brief explanation of the drawings]

[0010] [Figure 1] 1 is a diagrammatic overview of the problem of seizure propagation estimation according to the present invention. [Figure 2] 1 illustrates the training and application phases of the method of the present invention. [Figure 3] The estimated results obtained by the method of the present invention for a single seizure with strong coupling and 21 observation regions are shown. [Figure 4]1 shows the probability of having a seizure at a hidden node obtained by an implementation of the method of the present invention. [Figure 5] The method of the present invention shows the probability that the estimated onset time will be within T seconds of the true onset time in the hidden seizure region. [Figure 6A] The method of the present invention provides true and estimated regional excitability of the observed region. [Figure 6B] The method of the present invention provides true and estimated regional excitability of the observed region. [Figure 6C] The method of the present invention provides true and estimated regional excitability of the observed region. [Figure 7A] The method of the present invention shows the true and estimated node excitability of the hidden region. [Figure 7B] The method of the present invention shows the true and estimated node excitability of the hidden region. [Figure 7C] The method of the present invention shows the true and estimated node excitability of the hidden region. DETAILED DESCRIPTION OF THE INVENTION

[0011] Detailed Description of the Invention The present invention relates to a method for determining the onset time and / or excitability of brain regions not observed in epileptic patients, both during and without seizure activity. It attempts to estimate the activity of the entire brain network during a seizure from the observed activity of brain subnetworks and from the structural connectome. This method is based on the assumption that activity in any brain region can be classified as either a normal or seizure state, and that the transition time between these states (region onset time) is the only relevant feature for estimating propagation patterns. The method consists of two phases: in the training phase, global parameters of the model, i.e., parameters shared across different subjects and across seizures, are estimated from the training data. In the application phase, the trained model is applied to single seizure data, and the activity of the entire brain network is estimated for that seizure. This data-driven approach eliminates the need for manual global or patient-specific tuning, requiring only a few clearly interpretable constants.

[0012] As shown in Figure 1, the brain network state during the spread of an epileptic seizure is only partially understood. Some brain regions are involved in the onset of the seizure at time t i Some brain regions are observed to be involved in seizure activity in the brains of epilepsy patients. Other brain regions are observed not to be involved in seizure activity. And still other regions are hidden, i.e., their involvement in seizure activity is unknown. These are brain regions that are not observed in the brains of epilepsy patients, whether or not seizure activity occurs.

[0013] As shown in Figure 2, our method can estimate the state of hidden nodes by exploiting the structural connectivity seen from diffusion-weighted images. Two phases of estimation (bottom): In the training phase, the hyperparameters of the model are learned from training data. In the application phase, the model is applied to a single seizure.

[0014] In an implementation of the method of the present invention, computerized brain networks are provided. These networks model various brain regions and the connectivity between said regions for a training cohort of brains of epilepsy patients and the patient brains. The computerized brain networks are obtained from magnetic resonance neuroimaging and / or diffusion-weighted magnetic resonance imaging data. Similarly, datasets of observations of brain network states during epileptic seizures are provided for the training cohort and the patient brains, the observations defining regions within the brain network at an onset time as follows: a first region observed when the patient's brain exhibits epileptic seizure activity; a second region observed when the patient's brain does not exhibit said epileptic seizure activity; and a third region not observed whether the patient's brain exhibits said seizure activity or not.

[0015] These thus form the input data that comprises the connectome matrix. The computerized brain network is obtained from magnetic resonance neuroimaging and / or diffusion-weighted magnetic resonance imaging data. Individual structural connectome matrices are obtained from patient T1 and diffusion-weighted MRI scans, for example, using existing known and available software. Methods for connectome reconstruction are generally known in the art.

[0016] These input data also include the onset time of the region, which is information about whether and when the region enters the seizure state. To obtain such information data, the following steps are performed: First, intracranial seizure recordings are provided and, in each bipolar channel, the onset time is detected by evaluating the time when the signal power exceeds a given threshold. Specifically, a whitened power in the band [4,60] Hz and a log ratio threshold of 2 relative to the pre-ictal baseline are used. Note, however, that other implementations are possible, either by using different frequency bands or by using other approaches, such as machine learning-based detection methods. Second, all bipolar channels are assigned to brain regions. For each channel, the closest brain region is selected, as determined by the volumetric parcellation provided by the known FreeSurfer™ software and the location of the electrode contacts determined from a CT scan of the brain with the implanted electrodes. If there are multiple regions close to the contact and with a distance ratio of less than 2, the region is not assigned and information from the channel is not used. - Third, regions with no assigned channel are marked as hidden. For regions with one assigned channel, the origin time of the detected channel is used as the origin time of the region. If multiple channels are assigned to a region, the median of the detected origin times is taken.

[0017] According to the present invention, a dynamical model of brain network development and epileptic seizure propagation is provided, which describes the spread of a single region of the network to a delay variable with an activation function that is a function of the excitability of the region and network effects, and the onset time of a region is defined as the time when the delay variable crosses a given threshold.

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

number

[0019] In this formula, z i is a delay variable for region i, similar to the delay variable in the Epileptor 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 :

number

number

number

[0020] Function f q Since is positive, the model suggests that each region will have a finite time for seizure onset. In reality, this is incorrect, and the statistical model takes this into account by introducing a time constraint for seizures, and considering all regions with onset times that have a larger constraint than non-seizures.

[0021] Given a known vector of excitability c, a parameter vector q, and a connectome matrix W, the model uniquely defines a vector of origin times t. It uses the following shorthand notation for this mapping, which corresponds to equation (2) below: P(c,q,W)=t (2)

[0022] Activation function f q For the parameterization of , the exponential function of bilinearity is f q (c,y) parameterization. The interpolation point [c a ,y a ], [c a ,y b ], [c b ,y a ], [c b ,y b ](ca =-1, c b =1, y a =0, y b =1) in which the four coefficients q=(q 11 ,q 12 ,q 21 ,q 22 ), the function becomes

number

[0023] To ensure that f increases c, an additional constraint on the coefficients: q 21 >q 11 and q 22 >q 12 is set.

[0024] The method according to the invention further comprises the steps of providing a statistical model defining the probability of producing said set of observations of brain network states by said mechanical model; and training a mechanical model of epileptic seizure propagation using the statistical model and a dataset of observations of a training cohort. This is the training phase. Training the mechanical model is equivalent to finding an optimal set of parameters of the activation function.

[0025] The statistical model is a hierarchical model, constructed using Bayesian principles, with the top parameters being activation function parameters and the bottom parameters being regional excitability. However, this is just one example. This statistical model includes the assumption that the excitability of all regions across all training data has the same prior distribution.

[0026] In the training phase, data from multiple seizures are used to estimate the optimal value of a parameter vector q that is shared among all seizures. A statistical model for model training is constructed in the following way: for each seizure, an excitability vector c kIt is assumed that the values ​​of follow a standard normal distribution, and the activation function parameter q is treated as a hyperparameter of the model. Furthermore, the observed starting times have a standard deviation σ t Finally, to address the non-seizure region, we set the limit t lim is set and each region that has a seizure after this limit is treated as a non-seizure. The complete multi-seizure statistical model is shown below: Input data: connectome matrix W k , a set of seizure nodes and a set of non-seizure nodes, the onset time of the seizure nodes t(-) k,seizing . Parameter: σ q = 30, σ t =5s,t lim =90s Model:

number

[0027] This statistical model is based on the hyperparameter q and the seizure-specific excitability parameter c k We define a posterior probability distribution for q. Samples from this probability distribution are drawn using a Hamiltonian Monte Carlo method using a sampler implemented in the Stan™ software. In this phase, we want to obtain point estimates for the hyperparameter q, so the median of all samples is taken for each component of the parameter vector.

[0028] The method further includes the steps of inverting the trained dynamical model and utilizing a statistical model to estimate the onset time and / or excitability of a third region, or indeed at least one or all of the third regions, from the onset times observed in the first and second regions. This is the application phase. In this phase, activation function parameters of the dynamical model, as well as region excitability and onset times, are estimated using a Hamiltonian Monte Carlo method.

[0029] For the application phase, the activation function f q Once the parameters of are learned from a training batch of seizure data, additional seizures can be processed one by one using the single seizure model as follows: Input data: connectome matrix W, a set of seizure nodes and a set of non-seizure nodes, and the onset time of the seizure nodes t(-) k,seizing , the parameter q of the activation function. Parameter: σ q = 30, σ t =5s,t lim =90s Model:

number

[0030] This model is a simplification of the multiple seizure model using knowledge of the parameter vector q. Samples from the posterior distribution of the single seizure model can be drawn using Hamiltonian Monte Carlo methods, just like the multiple seizure model. The results of the single seizure estimation are the sampled excitability c and onset time t.

[0031] In this section, we present results showing the performance of our method on synthetic data generated by the same model. These results obtained using our method suggest that if the assumptions of our method are compatible with real epileptic networks, it may be possible to estimate hidden states. Indeed, our method can predict whether a hidden region will undergo a seizure. In addition, if a hidden region will undergo a seizure, our method can accurately predict when. Finally, our method can restore regional excitability.

[0032] Example 1: Test Data For testing purposes, synthetic data were generated using the same model used for estimation. We used a connectome matrix from 10 subjects, containing 84 nodes, using Dessicane-Killiani parcellation. From these 10 connectomes, two groups containing 10 bouts were created, all with three different parameter sets for the activation function: one with no coupling (q = (-5.0, -5.0, -3.0, -3.0)), one with weak coupling (q = (-6.5, -3.0, -3.5, 12.0)), and one with strong coupling (q = (-11.2, 5.3, -6.2, 69.3)). Finally, the number of observation regions was set to 21, 42, and 63 from the 84 regions of the parcellation. In total, this gave (2 groups) × (3 coupling strengths) × (3 number of observation regions) = 18 groups, each with 10 bouts. In all examples, the excitability c is randomly drawn from a standard normal distribution and the observed node is selected randomly. Estimation is performed separately for each of the 18 groups.

[0033] Example 2: Results Exemplary results for the estimation of a single seizure are partially shown in Figure 3, and the cumulative results for all tests are presented in Figures 4, 5, 6A, 6B, 6C, 7A, 7B and 7C.

[0034] In Figure 3, an example of the estimation results is partially shown for a single seizure with strong coupling and 21 observed regions. The left partial panel shows the true (dots) and estimated (plots) excitability of 84 brain regions. Black dots mark observed regions, and white dots mark hidden regions. Light grey regions are regions with seizures, and dark grey are regions without seizures. The numbers in the left column indicate

number

[0035] Figure 4 shows how well hidden regions can predict whether they will undergo seizures. It shows the estimated probability of seizure activity at hidden nodes. In this grid, the rows represent three different coupling strengths, and the columns represent different numbers of observed regions. Each panel in the grid represents the results of estimation in two groups, each with 10 subjects. The histogram counts the probability that a hidden region is involved in seizure activity. In light gray are histograms of seizure nodes, and in dark gray are histograms of non-seizure nodes. The inset shows the classification measure when regions are classified as seizure or non-seizure based on the estimated probability at a threshold of 0.5. The results show that seizure and non-seizure nodes can be reliably distinguished with only strong coupling. The number of observed nodes (nobs) does not necessarily increase the classification quality.

[0036] Figure 5 further illustrates how well the exact onset time can be estimated in the seizure region. It shows the probability that the estimated onset time in the hidden seizure region will be within T seconds of the true onset time. This probability is calculated as the fraction of samples with onset times within the window.

number

[0037] Figures 6A, 6B, and 6C show how regional excitability can be restored in unobserved nodes. They show the true and estimated regional excitability in 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 estimated excitability (y-axis) in seizure (light gray) and non-seizure (dark gray) regions. A perfect estimate would have all dots along the diagonal. The inserted numbers represent the R of the linear regression. 2 The coefficients are shown. The estimation is very good without coupling, and when strong coupling is introduced, some precision is lost due to network effects.

[0038] Figures 7A, 7B, and 7C show the same thing for hidden nodes. They show the true and estimated node excitability of hidden regions. The layout is the same as Figures 6A, 6B, and 6C. The excitability of hidden regions cannot be estimated in the absence of coupling (no laterals). In the case of strong coupling, there is a specific relationship between the true and estimated excitability (lateral strength) that increases with the number of observed regions.

[0039] Therefore, according to the present invention, with a strong coupling effect, it is possible to estimate whether a hidden node will cause a seizure, as shown in Figure 4. When there is a strong coupling effect, the onset time of the seizure node can also be predicted, as shown in Figure 5. For observed nodes, the excitability can be well estimated, especially with a weak coupling effect, as shown in Figures 6A-6C. With strong coupling and a large number of observed nodes, the excitability of hidden nodes can be estimated to a certain degree, as shown in Figures 7A-7C.

[0040] Finally, partial observation of seizure spread within a brain network, supplemented by the assumption that epileptic seizures propagate along structural connections that can be inferred from diffusion-weighted MRI, can be used to estimate whether and when hidden regions will manifest in seizures. To this end, a data-driven dynamical model of seizure onset and propagation throughout the entire enhanced network is employed. This simple dynamical model is enhanced by strong nonlinearities caused by intermittent changes in node states from normal to seizure states. To invert this model, a Bayesian inference framework is preferentially employed, and its parameters are estimated using a Hamiltonian Monte Carlo method, for example, using Stan software. Computational experiments on synthetic data generated by this model show that the quality of the estimation results depends on the number of observed nodes and the strength of connections. While the accuracy of the estimated onset time increases with the number of observed nodes and the strength of network effects, the model's ability to estimate the excitability of hidden nodes remains limited in all studies. These results show that in real epileptic seizures, the state of hidden nodes during the seizure can be estimated from incomplete observations.

Claims

1. 1. A computer-implemented method for determining the onset time and / or excitability of unobserved brain regions involved or not involved in seizure activity in the brain of an epilepsy patient, comprising: providing a computerized brain network that models various brain regions and connectivity between said regions for a training cohort of brains of epilepsy patients; providing a dataset of observations of brain network states during epileptic seizures for the training cohort; providing a dynamical model of epileptic seizure propagation in the brain network, the dynamical model describing the evolution of a slow variable of a single region of the brain network by a parameterized activation function that is a function of the excitability of the region, and a region onset time is defined as the time at which the slow variable crosses a given threshold; providing a statistical model that defines the probability of generating the dataset of observations of the brain network states according to the dynamical model; training the dynamical model using the statistical model and the dataset of observations of the training cohort to determine an optimal set of parameters of the activation function; providing a dataset of observations of brain network states during epileptic seizures for the brains of epilepsy patients, said observations classifying regions in said brain networks as: a first region of said patient's brain observed to be involved in epileptic seizure activity at time of onset; a second region of said patient's brain observed not to be involved in said epileptic seizure activity; and a third region of the patient's brain that is not observed to be involved or uninvolved in the seizure activity; and applying the trained dynamic model inversely to the brain of the epileptic patient as a single seizure model, and estimating the onset time and excitability of a third region from the onset time observed in the first region using the statistical model; Including, a dataset of observations of brain network states during epileptic seizures 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 dynamical model of epileptic seizure propagation in the brain network is defined, for a brain network with n regions, by the following equation: [Equation 1] where z i is the delay variable of domain i, and the function [Equation 2] is the activation function parameterized by the parameter vector q, c i is the node excitability, and W=(ω ij ) is [Equation 3] H is a Heaviside step function representing the switch from normal to seizure state when the delay variable exceeds the threshold z=1; The activation function f q is a bilinear function with exponentiation, and the parameters of the activation function are [Equation 4] are the values ​​of the activation function f q at four specified points in To parameterize the activation function fq, a bilinear exponential function fq(c,y) is used, parameterized by four coefficients q=(q11, q12, q21, q22) among the interpolation points [c a , y a ], [c a , y b ], [c b , y a ], [c b , y b ] (c a =-1, c b =1, y a =0, ​​y b =1), and the activation function is given by [Equation 5] where q 21 >q 11 and q 22 >q 12 ; In the training step, estimating an optimal value of a parameter vector q that is shared among all seizures in the training cohort; The statistical model: [Equation 6] where k is a node, t(-)k,sizing is the onset time observed at a seizure node, σ q =30, σ t =5s, and limit t lim =90s; P(c k ,q,W k ) represents the dynamical model that maps excitability c k , parameter q, and connectome matrix W k onto the onset time t k ; the statistical model defines the posterior probability distribution of hyperparameter q and seizure-specific excitability parameter c k ; The single seizure model comprises: [Equation 7] where W is the connectome matrix, q is the determined optimal set of parameters of the activation function, t(-)k,sizing is the onset time of the seizure node, σ q =30, σ t =5s, t lim =90s, Thereby, the onset time t and excitability c of the third brain region of the epilepsy patient are estimated. method.

2. The method of claim 1 , wherein the computerized brain network is obtained from magnetic resonance neuroimaging and / or diffusion-weighted magnetic resonance imaging data.

3. 10. The method of claim 1, wherein the seizure onset detection algorithm utilizes time-frequency analysis of intracranial electroencephalogram signals.

4. 4. The method of claim 1 or claim 3, wherein the mapping of the time of origin detected in the intracranial EEG signal to the brain region is based on the physical distance between the electrode contacts and the brain region recording the signal.

5. The method according to any one of claims 1 to 4, wherein the parameters of the activation function of the dynamical model for the third region and the excitability and onset time are estimated by a Hamiltonian Monte Carlo method.

6. The method of any one of claims 1 to 5, wherein the dynamic model is data-driven.

7. The method according to any one of claims 1 to 6, wherein the statistical model is a hierarchical model constructed according to the principles of Bayesian estimation, in which the top-level parameters are parameters of the activation function and the bottom-level parameters are regional excitability.

8. The method of any one of claims 1 to 7, wherein the statistical model includes the assumption that the excitability of all regions for all training data has the same prior distribution.

Images