Epileptic seizure prediction method combining electroencephalogram and image data

Abstract

The invention relates to the technical field of biomedical signal processing and neural image analysis, in particular to an epileptic seizure prediction method combining electroencephalogram and image data, which comprises the following steps: acquiring diffusion tensor imaging data of a target object, calculating the diffusion tensor of whole brain voxels, and calculating the image data of the target object; an anisotropic transmission resistance tensor field representing the conduction performance of the white matter fiber bundle of the brain is constructed based on inverse transformation of the diffusion tensor; and based on the anisotropic transmission resistance tensor field, calculating an anisotropic geodesic distance between any two points in the cerebral cortex space, and constructing a Riemannian manifold geometric substrate containing anatomical conduction constraint. According to the method, by constructing the anisotropic transmission resistance field and the Riemannian manifold substrate, the minimum biophysical cost constrained by anatomy is calculated, so that the problems that correlation analysis lacking physical constraint is mostly adopted in the traditional technology, and due to the neglect of the white matter fiber conduction rule, the analysis accuracy is low are solved; therefore, the false alarm problem caused by false synchronization caused by the fact that the volume conductor effect cannot be eliminated is solved.

Description

Technical Field

[0001] This invention relates to the field of biomedical signal processing and neuroimaging analysis technology, and in particular to a method for predicting epileptic seizures by combining electroencephalogram (EEG) and imaging data. Background Technology

[0002] In recent years, with the rapid development of neural engineering and artificial intelligence technologies, epilepsy diagnosis and treatment technology is undergoing a paradigm shift from post-event detection to pre-event prediction. In the field of epilepsy seizure prediction, early techniques primarily relied on single-modal scalp EEG signal analysis. These methods focused on extracting time-domain, frequency-domain, or nonlinear dynamic characteristics of EEG signals, attempting to capture aura states that could predict seizures minutes to hours before they occur. However, traditional techniques mostly employ correlation analysis lacking physical constraints. Because they ignore white matter fiber conduction patterns, they fail to eliminate false alarms caused by spurious synchronization resulting from volume conductor effects. Summary of the Invention

[0003] To overcome the above shortcomings, this invention provides a method for predicting epileptic seizures that combines EEG and imaging data. This method aims to improve upon the problem that traditional techniques mostly rely on correlation analysis lacking physical constraints. Due to the neglect of white matter fiber conduction patterns, false alarms caused by the inability to eliminate spurious synchronization resulting from volume conductor effects are not eliminated.

[0004] This invention provides the following technical solution: a method for predicting epileptic seizures by combining electroencephalogram (EEG) and imaging data, comprising the following steps: S1. Obtain diffusion tensor imaging data of the target object, calculate the diffusion tensor of the whole brain voxels, and construct an anisotropic transmission resistance tensor field characterizing the conduction performance of the brain white matter fiber bundles based on the inverse transformation of the diffusion tensor. S2. Based on the anisotropic transmission resistance tensor field, calculate the anisotropic geodesic distance between any two points in the cerebral cortex space, and construct a Riemannian manifold geometric basis containing anatomical conduction constraints. S3. Acquire real-time multi-channel EEG signals of the target object, extract energy data of characteristic frequency bands, and combine them with the spatial position of the electrodes to map them into a non-equilibrium dynamic energy measure distribution defined on the Riemannian manifold geometric basis. S4. Construct an abnormal flow evolution model that connects the current dynamic energy measurement distribution with the preset attack template measurement distribution. This model decomposes the state transition into an advection transport process constrained by the anisotropic transport resistance tensor field and a reactive growth process constrained by the excitability of local brain regions. S5. Define an energy functional that includes a transmission cost term and a growth cost term. By solving the infimum of the energy functional, calculate the minimum biophysical cost required to evolve the current state into the eruption template state. S6. Establish the mapping relationship between the minimum biophysical cost and epilepsy susceptibility, generate a real-time epilepsy susceptibility index, and output an epileptic seizure prediction signal when the index meets the preset alarm conditions.

[0005] By adopting the above technical solution, an anisotropic transmission resistance field and a Riemannian manifold basis are constructed, and then the minimum biophysical cost under anatomical constraints is calculated. This improves the problem that traditional techniques mostly use correlation analysis that lacks physical constraints. Due to the neglect of the conduction law of white matter fibers, false alarms caused by the inability to eliminate the false synchronization caused by the volume conductor effect are not eliminated.

[0006] Preferably, the anisotropic transport resistance tensor field includes: Traverse each voxel of the diffusion tensor imaging data and extract the set of eigenvalues ​​and the corresponding set of eigenvectors of the diffusion tensor at that voxel location. Construct an identity matrix containing preset regularization parameters, and perform a weighted summation with the diffusion tensor to obtain the regularized diffusion tensor; Perform matrix inversion on the regularized diffusion tensor to generate a local transmission resistance tensor at the voxel location, wherein the eigenvector direction of the tensor is consistent with that of the original diffusion tensor, and the eigenvalues ​​are reciprocals of the eigenvalues ​​of the original diffusion tensor.

[0007] Preferably, the Riemannian manifold geometric basis comprises: A three-dimensional mesh map covering the cerebral cortex region is established, and the anisotropic transmission resistance tensor field generated in the previous step is assigned to the corresponding node in the mesh as the metric tensor of the mesh space. Based on the anisotropic process function equation, wavefront evolution simulation is performed on the three-dimensional mesh diagram using the fast travel method; Calculate the cumulative time integral value of the wavefront propagating from the source node to the target node, and define this integral value as the anisotropic geodesic distance between the two points.

[0008] Preferably, the non-equilibrium dynamic energy measure distribution includes: Short-time Fourier transform was performed on the real-time acquired multi-channel EEG signals to extract the instantaneous power value of the Gamma band in the EEG signals; A discrete positive measure consisting of a set of weighted Dirac functions is constructed, where the spatial coordinates of each Dirac function correspond to the physical spatial coordinates of the EEG electrode, and its weight corresponds to the instantaneous power value. A preset background noise constant with a uniform global distribution is superimposed on the discrete positive measure to form a dynamic energy measure distribution with non-zero support.

[0009] Preferably, the abnormal flow evolution model includes: A partial differential equation constraint describing the state evolution is established, which sets out the linear relationship between the energy density with respect to time, the divergence of the translation flux, and the source term. The translational flux is defined as the product of energy density and the translational velocity vector field; The source term is defined as the product of energy density and scalar reaction rate field.

[0010] Preferably, the energy functional includes: Define a transmission cost density function, which is a quadratic form of the advection velocity vector with respect to the anisotropic transmission drag tensor constructed in the preceding steps; Define a growth cost density function, which is the product of the square of the scalar response rate and the location-dependent epileptogenic damping coefficient; The integral of the sum of the transmission cost density function and the growth cost density function in the spatiotemporal domain is defined as the energy functional.

[0011] Preferably, the epileptogenic damping coefficient includes: Based on anatomical structure, the brain space is divided into several anatomical regions of interest, and different damping scalar values ​​are preset for each region of interest; A first damping scalar value is set in the hippocampus and amygdala region, and a second damping scalar value is set in the frontal and parietal cortex regions. The value of the first damping scalar is less than the value of the second damping scalar.

[0012] Preferably, the infimum of solving the energy functional includes: A distance matrix is ​​constructed based on the anisotropic geodesic distances calculated in the previous steps, and the mixed cost matrix is ​​calculated using analytical formulas in conjunction with the epileptogenic damping coefficients determined in the previous steps. The mixed cost matrix is ​​negative and divided by a preset entropy regularization parameter, and then an exponential operation is performed to generate the Gibbs kernel matrix. Initialize two dual vectors with preset constant elements, and use the Gibbs kernel matrix to update the energy measure distribution vector and the attack template measure distribution vector at the current moment by alternating division; When the dual vector converges iteratively, the inner product of the dual vector and the Gibbs kernel matrix is ​​calculated to obtain the minimum biophysical cost.

[0013] Preferably, the mapping relationship includes: Obtain the minimum biophysical cost obtained by solving the energy functional; Calculate the reciprocal of the sum of the minimum biophysical cost and a preset constant; The reciprocal is used as the real-time epilepsy susceptibility index at the current moment.

[0014] Preferably, the preset alarm conditions include: Maintain a sliding time window that stores historical epilepsy susceptibility indices; Calculate the mean and standard deviation of the historical data within the sliding time window; The sum of the mean and the standard deviation of the preset multiple is set as the dynamic adaptive threshold; Compare the current epilepsy susceptibility index with the dynamic adaptive threshold. If the current index is greater than the threshold, the alarm condition is determined to be met.

[0015] The present invention has the following beneficial effects: 1. In this invention, by constructing an anisotropic transmission resistance field and a Riemannian manifold substrate, the minimum biophysical cost under anatomical constraints is calculated, thereby improving the problem of false alarms caused by false synchronization due to the inability to eliminate volume conductor effects caused by the neglect of white matter fiber conduction law, which is mostly used in correlation analysis lacking physical constraints in traditional techniques.

[0016] 2. In this invention, by mapping real-time multi-channel EEG signals into a non-equilibrium dynamic energy metric distribution and decomposing state transitions into a reaction growth process constrained by the excitability of local brain regions, the total mass of the metric distribution is allowed to dynamically fluctuate over time to simulate the energy burst mechanism during abnormal neuronal discharge. This improves upon the problem that traditional geometric manifold or optimal transmission methods, which mostly employ equilibrium transmission models based on the mass conservation assumption, cannot describe the non-conservative generation characteristics of neural energy during epileptic seizures, thus making it difficult to accurately quantify the evolution process of local discharge intensity in the epileptogenic focus.

[0017] 3. In this invention, by constructing an abnormal flow evolution model and defining an energy functional that includes a transmission cost term and a growth cost term, the evolution of brain states is decoupled into advection transmission along nerve fibers and reactive growth generated locally, and the minimum biophysical cost is calculated. This improves upon the fact that most traditional multimodal fusion prediction techniques use simple feature-level splicing or general deep learning models lacking mechanistic constraints. Because they fail to distinguish between the spatial diffusion of signals and the generation of source terms from a biophysical perspective, the models suffer from poor interpretability and insufficient ability to locate epileptogenic mechanisms when facing complex epilepsy propagation networks.

[0018] 4. In this invention, the minimum biophysical cost required to evolve the current state into a seizure template state is calculated by solving the infimum of the energy functional, and a real-time epilepsy susceptibility index is generated accordingly. This quantifies the ease with which the brain overcomes anatomical resistance and physiological inhibitory barriers to enter a seizure state. This improves upon traditional clinical monitoring methods, which mostly use threshold detection of waveform amplitude or specific rhythms. Because alarms can only be triggered after the scalp EEG signal shows obvious clinical seizure characteristics, there is a time lag in issuing early warning signals during the pre-seizure prodromal period when the energy barrier is lowered but the waveform has not yet shown significant abnormalities. Attached Figure Description

[0019] Figure 1 This is a flowchart of a method for predicting epileptic seizures that combines electroencephalogram (EEG) and imaging data, as proposed in this invention. Figure 2 This is a flowchart of the anisotropic transmission resistance tensor field construction for an epileptic seizure prediction method combining EEG and imaging data proposed in this invention. Figure 3 This is a flowchart of the Riemannian manifold geometric basis construction for a method for predicting epileptic seizures that combines EEG and imaging data, as proposed in this invention. Figure 4 This is a flowchart illustrating the construction of a non-equilibrium dynamic energy measure distribution for a method of predicting epileptic seizures by combining EEG and imaging data, as proposed in this invention. Figure 5 This is a flowchart of the minimum biophysical cost calculation for a method for predicting epileptic seizures that combines EEG and imaging data, as proposed in this invention. Figure 6 This is a flowchart of the prediction signal output of a method for predicting epileptic seizures that combines EEG and imaging data, as proposed in this invention. Detailed Implementation

[0020] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0021] Example 1: In a first embodiment of the present invention, the present invention provides a method for predicting epileptic seizures by combining electroencephalogram (EEG) and imaging data, such as... Figures 1-6 As shown, it includes the following steps: S1. Obtain diffusion tensor imaging data of the target object, calculate the diffusion tensor of the whole brain voxels, and construct an anisotropic transmission resistance tensor field characterizing the conduction performance of the brain white matter fiber bundles based on the inverse transformation of the diffusion tensor. Furthermore, the anisotropic transport resistance tensor field includes: Traverse each voxel of the diffusion tensor imaging data and extract the set of eigenvalues ​​and the corresponding set of eigenvectors of the diffusion tensor at that voxel location. Construct an identity matrix containing preset regularization parameters, and sum it with the diffusion tensor using weighted summation to obtain the regularized diffusion tensor; Perform matrix inversion on the regularized diffusion tensor to generate a local transport resistance tensor at the voxel location, wherein the eigenvector direction of the tensor is consistent with that of the original diffusion tensor, and the eigenvalues ​​are reciprocals of the eigenvalues ​​of the original diffusion tensor.

[0022] Specifically, the core of S1 lies in transforming the physical quantity characterizing the diffusion ability of water molecules into a geometric quantity characterizing the cost of signal transmission, thereby providing a physical constraint base for subsequent optimal transmission calculations.

[0023] The data processing system receives diffusion tensor imaging data of the target object at its input. This data is typically four-dimensional image data that has undergone eddy current correction, head motion correction, and noise reduction. For each voxel location within the whole-brain space... Diffusion tensor imaging data provides The symmetric positive definite matrix, i.e., the diffusion tensor. Diffusion tensor The three-dimensional Gaussian diffusion distribution of water molecules at this voxel location is described, comprising three orthogonal eigenvectors and three corresponding non-negative eigenvalues. The magnitude of the eigenvalues ​​reflects the diffusion intensity along the direction of the corresponding eigenvector. In the white matter fiber bundle region, the eigenvalues ​​along the fiber bundle direction are much larger than those perpendicular to the fiber bundle direction, exhibiting significant anisotropy.

[0024] Directly inverting the original diffusion tensor carries the risk of numerical instability. In cerebrospinal fluid regions or background noise regions, the diffusion tensor may approach a singular matrix, or its eigenvalues ​​in certain directions may be close to zero. To ensure the convergence of subsequent calculations, this step introduces regularization processing on the original diffusion tensor.

[0025] Construct an identity matrix containing preset regularization parameters, and calculate the regularized diffusion tensor accordingly. : ;in This represents the regularized diffusion tensor; Represents the original diffusion tensor; express The identity matrix; This indicates the preset regularization parameters. The range of values ​​is In a preferred embodiment, a preset regularization parameter is used. Values The average diffusion coefficient of water molecules in brain tissue is typically within a certain range. Magnitude. Selection to The range ensures that numerical divergence or singularity problems caused by eigenvalues ​​approaching zero are avoided during matrix inversion operations, while also ensuring that the introduced deviation is much smaller than the actual diffuse signal intensity, thereby preserving the anisotropic conduction characteristics of white matter fiber bundles to the maximum extent.

[0026] After obtaining the regularized diffusion tensor, the system performs a matrix inversion operation to generate the transmission resistance tensor. From the perspective of physical conduction, a larger diffusion coefficient means better conductivity of the medium, and the corresponding transmission resistance should be smaller; conversely, a smaller diffusion coefficient means poorer conductivity, and the transmission resistance should be larger.

[0027] Based on the calculated voxel position Local transmission resistance tensor at point : ;in That is, the anisotropic transport resistance tensor field at position The tensor value at that location; This represents the matrix inversion operation.

[0028] The local transmission resistance tensor generated by the above matrix inversion operation It inherits the directional properties of the original diffusion tensor and performs a reciprocal mapping on its amplitude properties.

[0029] Let the original diffusion tensor be... Feature decomposition into ,in For including eigenvalues diagonal matrix, Let be the eigenvector matrix. Then the generated transmission resistance tensor eigenvalues Satisfying Relationship: ;in For the original diffusion tensor, the first Diffusion coefficient in each principal direction, The above regularization parameters are used.

[0030] This results in a technical constraint: in the principal direction of the nerve fiber bundle, the original diffusion coefficient... The calculated resistance characteristic value is relatively large. Smaller. This corresponds to a shorter metric length in Riemannian manifold geometry, meaning lower cost for signal propagation along this direction. In the direction perpendicular to the nerve fiber bundle, the original diffusion coefficient... The calculated resistance characteristic value is relatively small. Larger. In Riemannian manifold geometry, this corresponds to a longer metric length, meaning that the cost of signal propagation along this direction is high.

[0031] After traversing and calculating the locations of all voxels in the whole brain, the system outputs a three-dimensional tensor field data structure, namely the anisotropic transmission resistance tensor field. This tensor field serves as the input data for the construction of the Riemannian manifold geometric basis in subsequent steps. It directly defines the metric tensor in the optimal transmission path integral, forcing the subsequently calculated geodesic paths to preferentially follow the anatomical paths of the brain's white matter fiber tracts, rather than straight lines in Euclidean space.

[0032] S2. Based on the anisotropic transmission resistance tensor field, calculate the anisotropic geodesic distance between any two points in the cerebral cortex space, and construct the Riemannian manifold geometric basis containing anatomical conduction constraints. Furthermore, the geometric basis of Riemannian manifolds includes: A three-dimensional mesh map covering the cerebral cortex region is established, and the anisotropic transmission resistance tensor field generated in the previous step is assigned to the corresponding node in the mesh as the metric tensor of the mesh space. Based on the anisotropic process function equation, wavefront evolution simulation is performed on a three-dimensional mesh diagram using the fast marching method; Calculate the cumulative time integral value of the wavefront propagating from the source node to the target node, and define this integral value as the anisotropic geodesic distance between the two points.

[0033] Specifically, S2 receives the anisotropic transmission resistance tensor field output by S1 and uses numerical calculation methods to transform the local tensor properties of the physical layer into the global distance metric of the geometric layer, thereby constructing a Riemannian manifold space that reflects the connectivity characteristics of white matter fiber bundles in the brain.

[0034] The data processing system receives anatomical data of the cerebral cortex and anisotropic transport resistance tensor field generated in the preceding steps. First, a three-dimensional mesh is constructed based on the morphological data of the cerebral cortex. This mesh consists of a set of nodes and a set of edges, where the node set covers the cerebral cortex and white matter regions. The spatial resolution between nodes is consistent with or interpolated to match the voxel resolution of the diffusion tensor imaging data.

[0035] Subsequently, the system performs a metric tensor assignment operation. For each node in the 3D mesh graph, the system retrieves its corresponding anisotropic transport resistance tensor in physical space coordinates. This tensor is assigned to the corresponding node as the Riemannian metric tensor of that local space. Thus, the three-dimensional mesh is transformed into a discrete Riemannian manifold structure, where the length measurement standard at any point in its space is no longer the isotropic Euclidean distance, but rather controlled by the transmission resistance tensor at that point.

[0036] To calculate the shortest path between any two points on a manifold, this embodiment employs a wavefront evolution simulation method. Unlike uniform wavefront diffusion in isotropic media, in anisotropic transport resistance tensor fields, the propagation speed of the wavefront is highly correlated with its direction.

[0037] The propagation law of wavefronts is described based on anisotropic equations. Let... For wavefront from source node Propagate to any point in space The shortest time required, which is a time function that satisfies a nonlinear partial differential equation: ;in Represents the arrival time function At point The spatial gradient vector at that location; Point The original diffusion tensor at that point is the transmission resistance tensor. The inverse matrix; The expression represents the vector transpose operation; the value 1 represents the normalized velocity potential of wavefront propagation. This equation shows that the magnitude of the time gradient of wavefront propagation in a certain direction depends on the diffusion coefficient in that direction. The larger the diffusion coefficient, the smaller the time gradient, and the faster the wavefront propagates.

[0038] The system employs a fast traversal method to numerically solve the formula. The algorithm initializes the source node. Arrival time The arrival time of the first node is zero, and the arrival time of all other nodes is set to infinity. The algorithm maintains a narrowband set and updates the arrival time values ​​of neighboring nodes from the inside out according to the principle of prioritizing the minimum arrival time of the wavefront. During the update process, the cost of wavefront propagation along the grid edge is calculated using the local transmission resistance tensor of the grid nodes, until the wavefront covers all nodes in the target area.

[0039] Through the wavefront evolution simulation described above, the minimum cumulative time from the source node to the target node is obtained. Within the Riemannian geometric framework, this minimum cumulative time is equivalent to the geodesic distance between the two points.

[0040] Based on the definition of two points and Anisotropic geodesic distance between : ; in Indicates connection point and Arbitrary parameterized smooth path, parameters The value range is from 0 to 1; Representing a path At any moment The tangent vector, i.e., the propagation velocity vector; The term represents the anisotropic transport resistance tensor at a point on the path; the integral term represents the length of the path element weighted by the resistance tensor field. Minimization operator This means selecting the path with the smallest integral value from all possible paths connecting two points as the geodesic.

[0041] The system performs traversal calculations and outputs an anisotropic geodesic distance matrix or a distance query function. The values ​​in this distance matrix represent the geometric cost of signal transmission along the anatomical path of white matter fiber tracts between any two points in the cerebral cortex. If there is a direct and robust neural fiber connection between the two points, even if their Euclidean distance is large, the calculated anisotropic geodesic distance is small; if there is no direct anatomical connection between the two points, even if their Euclidean distance is small, the calculated distance is large because the wavefront needs to bypass a high-resistance region. This geometric basis data serves as the foundational metric space for constructing the abnormal flow evolution model and calculating the optimal transmission cost in subsequent step S3.

[0042] S3. Acquire real-time multi-channel EEG signals of the target object, extract energy data of characteristic frequency bands, and combine them with the spatial location of the electrodes to map them into a non-equilibrium dynamic energy measure distribution defined on the geometric basis of the Riemannian manifold. Furthermore, the distribution of non-equilibrium dynamic energy measures includes: Short-time Fourier transform was performed on the real-time acquired multi-channel EEG signals to extract the instantaneous power value of the Gamma band in the EEG signals; Construct a discrete positive measure consisting of a set of weighted Dirac functions, where the spatial coordinates of each Dirac function correspond to the physical spatial coordinates of the EEG electrodes, and its weights correspond to the instantaneous power values. A preset background noise constant with a uniform global distribution is superimposed on the discrete positive measure to form a dynamic energy measure distribution with non-zero support.

[0043] Specifically, S3 aims to transform discrete, one-dimensional time-series signals into mathematical measurement objects defined on a continuous geometric space, thereby establishing a mapping relationship between functional data and anatomical structures.

[0044] The data processing system receives real-time multi-channel EEG signals from the target object at its input. The system acquires the preset time window length and overlap rate, and performs a short-time Fourier transform on the EEG signals of each channel. This transform maps the time-domain signal to the time-frequency domain, generating spectral information that varies with time.

[0045] The spectrum data is extracted based on a preset frequency band range, specifically the frequency components of the Gamma band. The Gamma band is defined as signal components with frequencies between 30 Hz and 80 Hz. The signal energy within this band is calculated to obtain the instantaneous power value of each channel at the current moment.

[0046] Let the first Each brainwave channel at all times instantaneous power value Calculation formula: ; in Indicates the first Each channel at time and frequency The magnitude of the short-time Fourier transform coefficients at the specified location; This indicates the lower limit frequency of the Gamma band, with a value of 30 Hz. This represents the upper limit frequency of the Gamma band, with a value of 80 Hz; the integration operation represents the accumulation of the power spectral density within this frequency band.

[0047] The system retrieves the three-dimensional coordinate data of each EEG electrode in physical space. Based on the node distribution of the Riemannian manifold geometric basis, the electrode coordinates are mapped to corresponding point positions on the manifold. A discrete positive measure consisting of a set of weighted Dirac functions is constructed. In this measure, each support point of a Dirac function corresponds to the physical location of an EEG electrode, and the weight value at that location is equal to the calculated instantaneous power value. .

[0048] This process did not affect the instantaneous power value. Normalization is applied, meaning the sum of the power of each channel changes dynamically over time. This non-normalization characteristic allows the constructed measure to reflect fluctuations in the overall energy level of the brain, meaning that total mass is not conserved, thus meeting the computational requirements of non-equilibrium optimal transmission.

[0049] The measure consisting solely of the Dirac function is spatially sparse and has a value of zero at non-electrode locations. To satisfy the numerical stability requirements of the subsequent Sinkhorn algorithm in the logarithmic field, the measure must have a non-zero support property across the entire domain.

[0050] The system sets a preset background noise constant that is uniformly distributed across the entire domain. This constant is an extremely small positive real number used to simulate environmental thermal noise or basal neural activity. This background noise constant is then superimposed onto the aforementioned discrete positive measure to generate the final non-equilibrium dynamic energy measure distribution. .

[0051] Non-equilibrium dynamic energy measure distribution The mathematical expression formula is as follows: ; in Indicates time In spatial location The density value measured at the location; Indicates the total number of channels for brainwave signals; Indicates the first The spatial coordinate vectors of each electrode on the Riemannian manifold basis; Dirac The function takes the value of infinity when the independent variable is zero and has an integral value of 1; here it is used to represent the energy of a point source. For the first The instantaneous power value of each channel is used as the quality weight of the point source; This represents the preset background noise constant, which exists throughout the entire manifold space. It is evenly distributed on the surface.

[0052] After the above processing, the system outputs a continuously varying time-varying non-equilibrium dynamic energy measure distribution sequence. Each time frame in this sequence... Both contain two types of information: spatial topology information determined by the electrode locations, and functional strength information determined by the power in the Gamma band. This metric distribution data serves as the input state for the subsequent abnormal flow evolution model in S4, i.e., the source distribution for the optimal transmission problem.

[0053] S4. Construct an abnormal flow evolution model that connects the current dynamic energy measurement distribution with the preset attack template measurement distribution. This model decomposes the state transition into an advection transport process constrained by an anisotropic transport resistance tensor field and a reactive growth process constrained by the excitability of local brain regions. Furthermore, the abnormal flow evolution model includes: A partial differential equation constraint describing the state evolution is established, which sets out the linear relationship between the energy density with respect to time, the divergence of the translation flux, and the source term. The translational flux is defined as the product of the energy density and the translational velocity vector field; The source term is defined as the product of energy density and scalar reaction rate field.

[0054] Specifically, S4 establishes a dynamic equation describing how the brain electrical energy state evolves from the current observed distribution to the preset attack template distribution. This equation constitutes the mathematical constraints for subsequent calculations of biophysical costs.

[0055] As the target state input to the abnormal flow evolution model, the pre-defined attack template measure distribution is used. This is not a general mathematical model, but rather a baseline state distribution built based on historical pathological data or a standard case database of the target object. The template construction process includes the following specific operations: If the target subject has historical long-term video EEG monitoring data, the data processing system executes the following processing flow: First, it reviews the target subject's historical monitoring records and filters out clinically confirmed typical epileptic seizure events. At least three independent seizure events are selected, and a time window of 0.5 to 2.0 seconds after the onset of each seizure is extracted as an analysis segment. Second, for each analysis segment, the system uses the same short-time Fourier transform algorithm as in step S3 to extract the average power spectral density of each channel in the Gamma band. Finally, combining the physical spatial coordinates of the electrodes, the extracted power data is mapped to a discrete measure distribution of a single seizure, and the discrete measure distributions of all selected events are spatially aligned and arithmetically averaged to generate the final preset seizure template measure distribution. .

[0056] If the target subject lacks historical seizure data, the data processing system executes the following steps: First, it accesses a pre-defined standardized epilepsy EEG database. This database stores population average seizure data categorized by epilepsy type and lesion side. Second, the system matches the corresponding population average data based on the target subject's clinical diagnosis information. Finally, it maps the population average data to the target subject's current electrode layout coordinate system using a spatial interpolation algorithm, generating an initial preset seizure template metric distribution. This template is iteratively updated as the target subject's actual seizure data is subsequently collected.

[0057] The input terminal of the data processing system receives the current dynamic energy measurement distribution output in step S3. and the aforementioned pre-constructed pre-seizure template measure distribution Define a normalized virtual evolution time interval, denoted as . Within this time interval, a time-varying variable is introduced. Continuously changing intermediate state density field .

[0058] The system sets the boundary conditions for the evolution process: at the beginning of the virtual evolution time. At that time, the intermediate state density field Equal to the current dynamic energy measure distribution At the end of virtual evolution time At that time, the intermediate state density field Equal to the episodic template measure distribution .

[0059] To describe the state Based on the continuous variation patterns in geometric space and time dimensions, a generalized continuity equation is established as a constraint for the partial differential equation. This equation decomposes the change in the density field into the superposition of two types of physical processes: one is the flow of matter in space, i.e., advection transport; the other is the generation or destruction of matter in situ, i.e. reactive growth.

[0060] The constraints set for the partial differential equations are: ;in Represents the energy density field virtual time The partial derivative of the density, i.e., the local rate of change of the density; This represents the divergence operator, used to calculate the flux source and sink distribution of a vector field; Represents the translational flux vector field; Let represent the scalar source term field. This equation shows that the change in local energy density is determined by the inflow and outflow flux divergence and the locally generated source terms, reflecting the generalized form of the law of conservation of energy in open systems.

[0061] The system provides a specific definition for the translational flux term. Translational flux. It describes the spatial transport process of neural energy along the cerebral cortex and white matter fiber bundles.

[0062] According to the definition of horizontal flow: ;in Instantaneous energy density; This is the advection velocity vector field, belonging to the tangent vector field on the geometric basis of the Riemannian manifold. This definition establishes that the spatial movement of energy is achieved through the density distribution driven by the velocity field. In the subsequent cost calculation, this velocity field... It will be directly constrained by the anisotropic transmission resistance tensor field constructed in S1, that is, if the direction of the velocity vector deviates from the main direction of the nerve fiber bundle, it will result in extremely high resistance cost.

[0063] The system defines source terms in detail. Source Term It describes the non-conservative growth or decay of neural energy in local brain regions, corresponding to the abnormal neuronal discharge or inhibition process during an epileptic seizure.

[0064] Define the source term according to the formula: ;in Instantaneous energy density; For scalar reaction rate field. The value represents the relative growth rate of energy per unit density. When When it is greater than zero, it indicates energy generation, corresponding to the burst of energy during an attack; when... A value less than zero indicates energy dissipation. This definition establishes that changes in energy intensity are driven by a reaction rate field. This reaction rate field is used in subsequent cost calculations. It will be constrained by the excitability parameters of local brain regions.

[0065] Based on the above definitions, a complete framework for a metamorphic flow evolution model is constructed. This model outputs a set of constraint equations containing the variables to be solved, including the intermediate state density field. Advection velocity vector field and scalar reaction rate field This set of constraint equations will serve as the boundary of the feasible region for energy functional optimization in S5, ensuring that the minimum biophysical cost calculated subsequently strictly follows the laws of physical conservation and birth-death.

[0066] S5. Define an energy functional that includes a transmission cost term and a growth cost term. By solving the infimum of the energy functional, calculate the minimum biophysical cost required to evolve the current state into the eruption template state. Furthermore, the energy functionals include: Define a transmission cost density function, which is a quadratic form of the advection velocity vector with respect to the anisotropic transmission drag tensor constructed in the preceding steps; Define a growth cost density function, which is the product of the square of the scalar response rate and the location-dependent epileptogenic damping coefficient; The integral of the sum of the transmission cost density function and the growth cost density function in the spatiotemporal domain is defined as an energy functional.

[0067] The epileptogenic damping coefficient includes: Based on anatomical structure, the brain space is divided into several anatomical regions of interest, and different damping scalar values ​​are preset for each region of interest; A first damping scalar value is set in the hippocampus and amygdala region, and a second damping scalar value is set in the frontal and parietal cortex regions. The value of the first damping scalar is less than the value of the second damping scalar.

[0068] Solving for the infimum of the energy functional includes: A distance matrix is ​​constructed based on the anisotropic geodesic distances calculated in the previous steps, and the mixed cost matrix is ​​calculated using analytical formulas in conjunction with the epileptogenic damping coefficients determined in the previous steps. The mixed cost matrix is ​​negative and divided by a preset entropy regularization parameter, and then an exponential operation is performed to generate the Gibbs kernel matrix. Initialize two dual vectors with preset constant elements, and use the Gibbs kernel matrix to update the energy measure distribution vector and the eruption template measure distribution vector at the current moment by alternating division; When the dual vector converges iteratively, the inner product of the dual vector and the Gibbs kernel matrix is ​​calculated to obtain the minimum biophysical cost.

[0069] Specifically, S5 integrates the geometric constraints, functional states, and dynamic models constructed in the preceding steps into a single scalar value—the minimum biophysical cost—by solving a variational problem. This value quantifies the anatomical resistance and physiological barriers that the brain must overcome to evolve from its current observed state to the template state for an attack.

[0070] The data processing system first establishes an energy functional to evaluate the total cost of any evolutionary path connecting the current energy measure distribution and the eruption template measure distribution. This energy functional... It consists of two parts: the transmission cost term and the growth cost term.

[0071] Define the transmission cost density function This function characterizes the kinetic energy loss during spatial movement, and its mathematical form is a quadratic form of the anisotropic transport drag tensor constructed from the advection velocity vector with respect to the preceding S1. Specifically: ;in For virtual time and location The advection velocity vector at that location; For position Anisotropic transport resistance tensor at the location; This represents the transpose of a vector. Because... It is generated based on the inverse transformation of the diffusion tensor, when the velocity vector... When the velocity vector is parallel to the main direction of the nerve fiber bundle, the quadratic value is smaller; when the velocity vector is perpendicular to the main direction of the nerve fiber bundle, the quadratic value is larger.

[0072] Define the growth cost density function This function characterizes the potential energy consumed by locally generated energy, and its mathematical form is the product of the square of the scalar reaction rate and the location-dependent epileptogenic damping coefficient. Specifically: ;in For virtual time and location scalar reaction rate at the point; For position The epileptogenic damping coefficient at the location.

[0073] The integral of the sum of the transport cost density function and the growth cost density function in the spatiotemporal domain is defined as the total energy functional. : ;in Represents the spatial domain of the cerebral cortex; Energy density distribution in the intermediate state; integration interval For virtual evolution time.

[0074] When calculating the growth cost density function, the system uses the epileptogenic damping coefficient based on anatomical structure. Perform spatial non-uniform assignment.

[0075] The system retrieves a pre-stored brain anatomical atlas and divides the cerebral cortex into several regions of interest (ROIs). For each ROI, the system assigns a specific damping scalar value. In this embodiment, the system identifies the hippocampus and amygdala regions and sets a first damping scalar value for them. The system identifies the frontal and parietal cortical regions and assigns a second damping scalar value to them. .

[0076] The numerical settings satisfy the following relation: This numerical relationship establishes physical constraints: in epilepsy-susceptible areas such as the hippocampus, the cost per unit reaction rate is lower, meaning that energy is easily generated explosively in these areas; while in areas such as the frontal lobe, the cost per unit reaction rate is higher, meaning that the inhibition of energy generation is stronger.

[0077] Solving the energy functional The infimum is equivalent to solving the Wasserstein-Fisher-Rao distance in the unbalanced optimal transport problem. Since directly solving the infinite-dimensional variational problem is computationally complex, this embodiment employs the Sinkhorn iterative algorithm based on entropy regularization for numerical solution.

[0078] First, calculate the mixed cost matrix between discrete grid node pairs. This matrix incorporates the anisotropic geodesic distances calculated in step S2. and the aforementioned epileptogenic damping coefficient Based on the analytic properties of the Wasserstein-Fisher-Rao metric, for any two points... and Mixed cost elements The calculation shows: ;in The anisotropic geodesic distance between two points; The mean of the epileptogenic damping coefficients at the two points; The function is used to truncate the distance to prevent the independent variable of the cosine function from going out of bounds; It is the natural logarithm. This formula maps the transmission cost and growth cost to a unified coupling cost between each pair of points.

[0079] For the mixed cost matrix Perform entropy regularization transformation to generate the Gibbs kernel matrix. The calculation shows: ;in This is a preset entropy regularization parameter, and its value is a positive real number; It is an exponential function.

[0080] The system initializes two dual vectors. and The initial values ​​of its elements are all set to preset constants, such as 1. Then, the Gibbs kernel matrix is ​​used... The energy measure distribution vector at the current moment and the distribution vector of the episodic template measure Perform alternating division updates. Update rules: ;in Indicates the number of iterations; For S3 output Discrete vector form; for It is a discrete vector form; division is an element-wise division.

[0081] The iteration terminates when the rate of change of the dual vector falls below a preset convergence threshold. Calculate the minimum biophysical cost. : ;in This is the approximate lower bound value of the energy functional obtained by solving. This value is output to subsequent steps to generate the epilepsy susceptibility index.

[0082] S6. Establish a mapping relationship between minimum biophysical cost and epilepsy susceptibility, generate a real-time epilepsy susceptibility index, and output an epileptic seizure prediction signal when the index meets the preset alarm conditions. Furthermore, the mapping relationships include: Obtain the minimum biophysical cost obtained by solving the energy functional; Calculate the reciprocal of the sum of the minimum biophysical cost and a preset constant; The reciprocal is used as the real-time epilepsy susceptibility index at the current moment.

[0083] Preset alarm conditions include: Maintain a sliding time window that stores historical epilepsy susceptibility indices; Calculate the mean and standard deviation of the historical data within the sliding time window; Set the sum of the standard deviations of the mean and a preset multiple as the dynamic adaptive threshold; Compare the current epilepsy susceptibility index with the dynamic adaptive threshold. If the current index is greater than the threshold, the alarm condition is determined to be met.

[0084] Specifically, S6, as the end of the data processing flow, is responsible for transforming the abstract physical quantities calculated in the preceding steps into intuitive clinical monitoring indicators and executing automatic judgments based on statistical laws.

[0085] The data processing system receives the minimum biophysical cost value output by S5 at its input. This value, in a physical sense, represents the difficulty the system faces in overcoming anatomical resistance and physiological inhibitory barriers to enter an episodic state. To align with clinical monitoring practices, where higher values ​​generally indicate greater risk, an inverse proportional mapping function is constructed.

[0086] The minimum biophysical cost obtained by solving the energy functional is denoted as . The system introduces a preset constant, denoted as . This preset constant is a non-zero positive real number, and its value is set according to the floating-point precision of the calculation system. It is used to prevent abnormal numerical cases where the denominator is zero.

[0087] Based on the calculation of the real-time epilepsy susceptibility index at the current moment : ;in Indicates time Epilepsy susceptibility index; Indicates time The minimum biophysical cost; This represents a preset constant.

[0088] This calculation process establishes a negative correlation between the two. When the brain is in a stable state, the calculated biophysical cost decreases due to transmission resistance and growth damping. A larger value corresponds to a higher susceptibility index. Smaller, at a low baseline. When the brain enters the pre-ictal phase, the anisotropic transmission cost of functional connectivity decreases and local excitability increases, leading to... The susceptibility index decreases rapidly. It exhibits a non-linear upward trend.

[0089] To accommodate the baseline EEG fluctuations of different subjects and the same subject at different times, the system abandons the fixed threshold decision method and adopts dynamic statistical decision logic based on sliding time windows.

[0090] Maintain a first-in, first-out (FIFO) sliding time window queue in memory. This queue stores the epilepsy susceptibility index sequence calculated at historical time points. Let the length of the sliding time window be... The data set contained in the queue is .

[0091] Real-time calculation of the mean of the data within the sliding time window with standard deviation : ; in The number of sampling points included in the sliding time window; This represents the susceptibility index for historical moments.

[0092] The system sets a preset multiplier parameter, denoted as . This parameter determines the sensitivity of the detection algorithm, typically ranging from 3 to 5, corresponding to the 3-sigma criterion in statistics or a more stringent anomaly detection standard. The current dynamic adaptive threshold is calculated using the formula. : ;in For a moment Dynamic alarm threshold; The baseline mean; The standard deviation of the fluctuation; This is a preset multiple. This threshold fluctuates in real time according to the background fluctuations in the subject's brain state.

[0093] Perform a comparison logic operation to compare the real-time epilepsy susceptibility index calculated at the current moment. With the calculated dynamic adaptive threshold Compare them.

[0094] Judgment Logic: If The system determines that the current brain state falls within the background fluctuation range, does not trigger an alarm, and continues to update the sliding time window data. If The system determines that the current brain state shows a statistically significant abnormality, and that this abnormality is caused by a significant reduction in biophysical costs, thus meeting the alarm conditions.

[0095] When alarm conditions are met, the system outputs a seizure prediction signal. This signal is either a digital interrupt signal or message data containing the predicted probability of occurrence. This prediction signal is transmitted to an external display device or neuromodulation device to alert medical personnel or trigger closed-loop electrical stimulation intervention.

[0096] Example 2: In clinical monitoring scenarios, existing prediction techniques based on simple signal processing or general deep learning face severe challenges: due to the lack of physical constraints on brain anatomy, existing techniques struggle to effectively distinguish between the "volume conductor effect" caused by physiological artifacts such as chewing and blinking and the synchronization of actual neural conduction, leading to numerous false alarms caused by spurious connections. Furthermore, for seizures originating in deep brain regions such as the hippocampus, conventional scalp EEG analysis exhibits significant detection lag due to signal attenuation and cannot adapt to individual anatomical differences caused by pathological changes, making it difficult to provide accurate and robust early warnings based on the patient's unique biophysical characteristics. To address these issues, this invention provides a method for predicting epileptic seizures that combines EEG and imaging data, the structure of which is as follows... Figure 1 As shown. The specific implementation process of this method is as follows: An anisotropic transport resistance tensor field is constructed by inverse transforming diffusion tensor imaging data, forcing subsequent signal propagation models to follow the anatomical paths of brain white matter fiber tracts, rather than straight paths in Euclidean space. Based on this resistance field, anisotropic geodesic distances are calculated to construct a Riemannian manifold geometric basis, which can distinguish physically adjacent but anatomically disconnected brain regions, eliminating spurious synchronization interference caused by volume conductor effects.

[0097] Mapping EEG signals to a non-equilibrium dynamic energy metric distribution preserves the energy bursts and mass non-conservation characteristics caused by abnormal neuronal discharges during epileptic seizures. The abnormal flow evolution model decouples the transition process of brain states into anatomically restricted advection transport and locally excitability-restricted reactive growth, enabling independent quantification of the signal diffusion process along nerve fibers and the local burst process of the epileptogenic focus.

[0098] By solving the lower bound of the energy functional that includes transmission and growth costs, the obtained minimum biophysical cost numerically characterizes the energy barrier required for the brain to transition from the current state to a seizure state. The reduction in this cost directly reflects the loss of brain network stability and the increase in seizure tendency, providing a quantitative basis based on biophysical mechanisms for the generated epilepsy susceptibility index. This solves the problems of weak generalization ability and high false alarm rate caused by neglecting anatomical physical constraints in traditional purely data-driven methods.

[0099] Finally, it should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for predicting epileptic seizures by combining electroencephalogram (EEG) and imaging data, characterized in that, Includes the following steps: S1. Obtain diffusion tensor imaging data of the target object, calculate the diffusion tensor of the whole brain voxels, and construct an anisotropic transmission resistance tensor field characterizing the conduction performance of the brain white matter fiber bundles based on the inverse transformation of the diffusion tensor. S2. Based on the anisotropic transmission resistance tensor field, calculate the anisotropic geodesic distance between any two points in the cerebral cortex space, and construct a Riemannian manifold geometric basis containing anatomical conduction constraints. S3. Acquire real-time multi-channel EEG signals of the target object, extract energy data of characteristic frequency bands, and combine them with the spatial position of the electrodes to map them into a non-equilibrium dynamic energy measure distribution defined on the Riemannian manifold geometric basis. S4. Construct an abnormal flow evolution model that connects the current dynamic energy measurement distribution with the preset attack template measurement distribution. This model decomposes the state transition into an advection transport process constrained by the anisotropic transport resistance tensor field and a reactive growth process constrained by the excitability of local brain regions. S5. Define an energy functional that includes a transmission cost term and a growth cost term. By solving the infimum of the energy functional, calculate the minimum biophysical cost required to evolve the current state into the eruption template state. S6. Establish the mapping relationship between the minimum biophysical cost and epilepsy susceptibility, generate a real-time epilepsy susceptibility index, and output an epileptic seizure prediction signal when the index meets the preset alarm conditions.

2. The method for predicting epileptic seizures by combining EEG and imaging data according to claim 1, characterized in that, The anisotropic transport resistance tensor field includes: Traverse each voxel of the diffusion tensor imaging data and extract the set of eigenvalues ​​and the corresponding set of eigenvectors of the diffusion tensor at that voxel location. Construct an identity matrix containing preset regularization parameters, and perform a weighted summation with the diffusion tensor to obtain the regularized diffusion tensor; Perform matrix inversion on the regularized diffusion tensor to generate a local transmission resistance tensor at the voxel location, wherein the eigenvector direction of the tensor is consistent with that of the original diffusion tensor, and the eigenvalues ​​are reciprocals of the eigenvalues ​​of the original diffusion tensor.

3. The method for predicting epileptic seizures by combining EEG and imaging data according to claim 1, characterized in that, The Riemannian manifold geometry basis includes: A three-dimensional mesh map covering the cerebral cortex region is established, and the anisotropic transmission resistance tensor field generated in the previous step is assigned to the corresponding node in the mesh as the metric tensor of the mesh space. Based on the anisotropic process function equation, wavefront evolution simulation is performed on the three-dimensional mesh diagram using the fast travel method; Calculate the cumulative time integral value of the wavefront propagating from the source node to the target node, and define this integral value as the anisotropic geodesic distance between the two points.

4. The method for predicting epileptic seizures by combining EEG and imaging data according to claim 1, characterized in that, The non-equilibrium dynamic energy measure distribution includes: Short-time Fourier transform was performed on the real-time acquired multi-channel EEG signals to extract the instantaneous power value of the Gamma band in the EEG signals; A discrete positive measure consisting of a set of weighted Dirac functions is constructed, where the spatial coordinates of each Dirac function correspond to the physical spatial coordinates of the EEG electrode, and its weight corresponds to the instantaneous power value. A preset background noise constant with a uniform global distribution is superimposed on the discrete positive measure to form a dynamic energy measure distribution with non-zero support.

5. The method for predicting epileptic seizures by combining EEG and imaging data according to claim 1, characterized in that, The abnormal flow evolution model includes: A partial differential equation constraint describing the state evolution is established, which sets out the linear relationship between the energy density with respect to time, the divergence of the translation flux, and the source term. The translational flux is defined as the product of energy density and the translational velocity vector field; The source term is defined as the product of energy density and scalar reaction rate field.

6. The method for predicting epileptic seizures by combining EEG and imaging data according to claim 1, characterized in that, The energy functional includes: Define a transmission cost density function, which is a quadratic form of the advection velocity vector with respect to the anisotropic transmission drag tensor constructed in the preceding steps; Define a growth cost density function, which is the product of the square of the scalar response rate and the location-dependent epileptogenic damping coefficient; The integral of the sum of the transmission cost density function and the growth cost density function in the spatiotemporal domain is defined as the energy functional.

7. The method for predicting epileptic seizures by combining EEG and imaging data according to claim 6, characterized in that, The epileptogenic damping coefficient includes: Based on anatomical structure, the brain space is divided into several anatomical regions of interest, and different damping scalar values ​​are preset for each region of interest; A first damping scalar value is set in the hippocampus and amygdala region, and a second damping scalar value is set in the frontal and parietal cortex regions. The value of the first damping scalar is less than the value of the second damping scalar.

8. The method for predicting epileptic seizures by combining EEG and imaging data according to claim 1, characterized in that, The infimum of the energy functional solution includes: A distance matrix is ​​constructed based on the anisotropic geodesic distances calculated in the previous steps, and the mixed cost matrix is ​​calculated using analytical formulas in conjunction with the epileptogenic damping coefficients determined in the previous steps. The mixed cost matrix is ​​negative and divided by a preset entropy regularization parameter, and then an exponential operation is performed to generate the Gibbs kernel matrix. Initialize two dual vectors with preset constant elements, and use the Gibbs kernel matrix to update the energy measure distribution vector and the attack template measure distribution vector at the current moment by alternating division; When the dual vector converges iteratively, the inner product of the dual vector and the Gibbs kernel matrix is ​​calculated to obtain the minimum biophysical cost.

9. The method for predicting epileptic seizures by combining EEG and imaging data according to claim 1, characterized in that, The mapping relationship includes: Obtain the minimum biophysical cost obtained by solving the energy functional; Calculate the reciprocal of the sum of the minimum biophysical cost and a preset constant; The reciprocal is used as the real-time epilepsy susceptibility index at the current moment.

10. The method for predicting epileptic seizures by combining EEG and imaging data according to claim 1, characterized in that, The preset alarm conditions include: Maintain a sliding time window that stores historical epilepsy susceptibility indices; Calculate the mean and standard deviation of the historical data within the sliding time window; The sum of the mean and the standard deviation of the preset multiple is set as the dynamic adaptive threshold; Compare the current epilepsy susceptibility index with the dynamic adaptive threshold. If the current index is greater than the threshold, the alarm condition is determined to be met.

Images