A cluster analysis method for epileptic spikes based on magnetoencephalographic topographic coding

Abstract

The present invention discloses a cluster analysis method for epileptic spikes based on magnetoencephalogram (MEG) topography coding, which comprises the following steps: 1) calculating the global energy spectrum of each sample of multiple spike samples during the interictal period of the same patient; 2) performing dipole fitting on each sample and deleting samples with a GOF value less than a threshold; 3) extracting the MEG topography corresponding to the maximum moment of each global energy spectrum; 4) performing perceptual hash coding on the magnetometer topography and gradiometer topography in each extracted MEG topography to obtain corresponding magnetometer topography coding features and gradiometer topography coding features; 5) using the two features to respectively calculate a first distance matrix and a second distance matrix between samples and perform fusion; 6) using the fused distance matrix for clustering and deleting outliers within the class; 7) retaining classes whose number of spike samples within the class is greater than a set threshold H; and then superimposing and averaging the samples within the class to obtain spikes after superimposition and averaging.

Description

Technical Field

[0001] The present invention relates to the technical field of magnetoencephalography and epilepsy analysis, and in particular to a spike clustering method based on spatial coding features of magnetoencephalogram topography. Background Art

[0002] Noninvasive magnetoencephalography (MEG), characterized by high temporal resolution, is a key modern neurological imaging technique. The signals detected by MEG primarily originate from postsynaptic potentials. Because electromagnetic field transmission is unaffected by the medium used, source imaging using MEG signals offers greater localization accuracy than electroencephalography (EEG). Furthermore, its noninvasive nature makes MEG promising for clinical neuroscience applications.

[0003] The signals recorded by magnetoencephalography and electroencephalography (EEG) originate from the same neurophysiological activity, and both can be used to record and observe abnormal discharges in epilepsy. Numerous experiments and clinical observations have shown that patients with epilepsy experience abnormal epileptiform discharges (IEDs) between attacks. The brain region that causes these discharges is called the stimulatory focus. Studies have shown that the epileptogenic focus is closely related to the epileptogenic focus. Detecting and tracing IEDs from MEG can pinpoint the patient's stimulatory focus, providing an effective reference for surgical resection of epilepsy patients.

[0004] Averaging IEDs with the same spatiotemporal characteristics can reduce noise, improve the signal-to-noise ratio (SNR), and thus enhance localization accuracy. In clinical practice, experienced physicians are often required to visually identify IEDs. Previous studies have shown that clustering algorithms can automatically identify IEDs with the same spatiotemporal characteristics. However, traditional clustering algorithms directly construct a distance matrix for all channel time-series IED samples. Considering the complexity of MEG IED waveforms, in patients with focal epilepsy, the epileptogenic focus is often confined to a specific brain region, with observed IED signals occurring in only a few MEG channels. Directly using signals from all channels to calculate the distance matrix increases errors. When signal noise is high, the distance matrix cannot accurately calculate the distances between IEDs. Furthermore, although automated IED detection algorithms have improved IED identification efficiency in recent years, the results often suffer from high false positive rates and inaccurate detection. Traditional clustering algorithms have low tolerance for clustering IED samples with false positives, resulting in poor clustering results.

[0005] Therefore, we need to design new clustering analysis strategies to reduce the algorithm's sensitivity to false-positive IEDs samples and improve clustering accuracy. Summary of the Invention

[0006] To address the technical problems existing in the prior art, the present invention provides a method for clustering spike waves in epileptic IEDs based on spatial encoding of magnetoencephalographic topography. This method uses the specific pattern of spike dipole morphology between seizures for the same patient as a priori conditions to perform perceptual hash coding on the topography. The perceptual hash coding is then used to calculate a distance matrix for subsequent cluster analysis.

[0007] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0008] 1) For spike wave samples from a single patient's interictal period, calculate the global field power (GFP) corresponding to each spike wave sample and align the spike wave samples according to the time point corresponding to the maximum value in the global field power (GFP). That is, the time point corresponding to the maximum value in the global field power (GFP) of each sample is set to the same moment, thereby aligning the samples.

[0009] 2) Perform dipole fitting on each spike sample, calculate the goodness of fit (GOF), dipole spatial position coordinates, and dipole directional parameters for each spike sample, and use the calculated GOF value to delete spike samples with a GOF less than 0.3. Fitting each spike sample yields a dipole with spatial and directional parameters and goodness of fit.

[0010] 3) For each sample, extract the magnetoencephalogram (MEG) topography corresponding to the GFP maximum point; MEG topography includes magnetometer topography and gradiometer topography;

[0011] 4) performing perceptual hash coding on the obtained magnetometer topography map and gradiometer topography map to obtain corresponding 64-dimensional coding features, namely, the magnetometer topography coding features and the gradiometer topography coding features;

[0012] 5) The distance matrix between spike samples is calculated using the magnetometer topography encoding features and the gradiometer topography encoding features of IEDs respectively, and the distance matrices calculated by the two are fused to construct a distance matrix for clustering, and the fused distance matrix is ​​clustered using a hierarchical clustering algorithm;

[0013] 6) The number of cluster categories is determined using the silhouette coefficient method, that is, the number of categories corresponding to the local maximum of the silhouette coefficient is used as the number of clusters;

[0014] 7) Clustering results often contain outliers and points with small GOF values. First, the local outlier factor (LOF) algorithm is used to detect and delete outliers within the cluster. Then, spike samples with a GOF value less than 0.5 within the cluster are removed.

[0015] 8) Setting a threshold H, retaining the class with spike wave samples greater than H, and superimposing and averaging the samples within the class to obtain spike waves after superimposition and averaging. In the present invention, the threshold H=4.

[0016] Traditional clustering algorithms cluster based on the temporal waveforms of spikes. However, when clinicians manually group (cluster) the selected spikes, they extract the topographic map at the spike GFP maximum as the grouping basis. Because spikes from the same individual exhibit several specific dipole morphologies on the topographic map, and these specific morphologies recur, the present invention uses the doctor-grouping method to encode the topographic map at the GFP maximum and perform subsequent clustering. Because the signal-to-noise ratio of individual spike samples is low, stacking and averaging the samples in each category after clustering can improve the spike signal-to-noise ratio. The stacked and averaged results are used to more accurately trace the stimulation focus.

[0017] At present, the more popular spike wave screening method is an automated spike wave detection algorithm based on pattern recognition or deep learning. This type of algorithm will detect a large number of false positive samples. The topography of the false positive samples usually does not show a specific dipole morphology, and usually cannot be clustered into a class with a larger sample size. To address this problem, the present invention uses clustering on the results of the automated detection algorithm, retains the class with a larger sample size, and removes the class with a small sample size (which can be understood as a class mixed with false positive samples) to remove the false positive samples for subsequent analysis. If a class with a larger sample size shows a specific dipole morphology, this class will be retained and the class will be superimposed and averaged within the class. Therefore, the results obtained by the method of the present invention can greatly improve the accuracy of the analysis.

[0018] Furthermore, in step 1), the global energy spectrum describes the standard deviation of the magnetic field intensity / gradient value detected by the magnetoencephalography sensor at each moment:

[0019]

[0020] Where N is the number of sensors in the MEG magnetometer / gradiometer; i is the sensor number; u i is the magnetic field intensity / magnetic field gradient value collected by the i-th detector at a certain moment, The average value of the magnetic field intensity / magnetic field gradient detected by all sensors. The time point corresponding to the maximum value of GFP during the spike time is calculated using formula (1). The present invention sets two groups of sensors, with one magnetometer and two gradiometers set at each collection position.

[0021] Furthermore, as described in step 2), the samples with GOF less than 0.3 in the spike wave are deleted. The dipole has three direction parameters and three spatial position coordinates.

[0022] Furthermore, the topographic map described in step 3) is divided into a magnetometer topographic map and a gradiometer topographic map. There are two gradiometer sensors at the same location, and the gradient values ​​of the corresponding positions of the two gradiometers are calculated as shown in formula (2):

[0023]

[0024] Among them, g i1 and g i2 Indicates the magnetic field gradient value collected by two gradiometers at the same position; g i is the gradient value of the corresponding position calculated by formula (2).

[0025] Furthermore, the perceptual hash coding algorithm process described in step 4) is as follows:

[0026] 1. Downsample and resize the topographic maps (magnetometer and gradiometer) corresponding to the time point of the extracted GFP maximum to remove details and high-frequency information. For example, if the original image size is 64×64, downsampling it to 48×48 will preserve the basic image structure but blur the details.

[0027] 2. Perform discrete cosine transform (DCT) on the scaled image as shown in formula (3).

[0028]

[0029] Among them, x, y = 0, 1, ..., M-1, M represents the length and width of the image (only square images are considered here), x and y represent the coordinates of a pixel in the image, f(x, y) represents the coordinate value; u and v represent the coordinates in the image transformation domain, F(u, v) represents the value of the coordinate (u, v) in the transformation domain, After the topographic map is transformed, most of its information is compressed into a few transformation coefficients.

[0030] 3. Extract the low-frequency part of the topographic map transform domain.

[0031] 4. Calculate the mean of the intercepted area, and set the feature values ​​above the mean to 1 among the extracted 64-dimensional features, and set the feature values ​​below or equal to the mean to 0.

[0032] Furthermore, in step 5), the distance matrix is ​​calculated using the Hamming distance of formula (4) for the 64-dimensional encoding features of the IEDs of all samples.

[0033]

[0034] Among them, l i ,s i Represents the i-th eigenvalue of the encoded feature; It is an exclusive OR operation.

[0035] For the distance matrix calculated from the magnetometer topography map and the distance matrix calculated from the gradiometer topography map, the total distance matrix is ​​calculated using formula (5):

[0036] Matrix D =aMatrix D_Mag +(1-a)Matrix D_Grad (5)

[0037] Among them, Matrix D Represents the final calculated distance matrix; Matrix D_Mag Matrix represents the distance matrix of the magnetometer encoding features; D_Grad Represents the distance matrix calculated by the gradiometer encoding feature; α is the weight coefficient, ranging from 0 to 1. Using the final distance matrix, a hierarchical clustering algorithm is used to cluster the samples.

[0038] Furthermore, the number of categories of the clustering algorithm is determined using the Silhouette Coefficient method described in step 6), and the number of categories corresponding to the local maximum value in the Silhouette Coefficient curve is selected as the final number of categories of the hierarchical clustering algorithm.

[0039] Furthermore, a higher GOF threshold (e.g., 0.5) is set compared to the threshold in step 2) and samples with a GOF value below this threshold (GOF∈[0,1]) are deleted. For the main cluster obtained by clustering, the local outlier factor algorithm is used to remove outliers within the main cluster. The local outlier factor algorithm requires the construction of a distance matrix, which is calculated using the dipole direction and position parameters obtained from the dipole fitting in step 2).

[0040] Furthermore, clustering results are generated according to the above steps. The main class in the clustering results will be selected, and for the main class samples, the waveforms of each channel are superimposed and averaged for subsequent traceability analysis.

[0041] The present invention provides a clustering analysis method using the information encoding features of spike wave topography patterns, which can cluster spike waves generated by interictal discharges collected by magnetoencephalography. The algorithm has the following advantages:

[0042] 1) Taking advantage of the specific topographic patterns of interictal discharges in epileptic patients, this study used magnetoencephalographic topographic information to cluster spikes for the first time.

[0043] 2) The algorithm uses the parameters generated by dipole fitting, such as goodness of fit, dipole direction, position and other information, combined with the outlier detection algorithm to eliminate abnormal samples;

[0044] 3) The results generated by this cluster analysis method can be used for tracing the source of epileptic lesions and have a stronger signal-to-noise ratio than single samples. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 It is a flow chart of the clustering algorithm of the present invention.

[0046] Figure 2 It is a flowchart of perceptual hash coding magnetoencephalography topography.

[0047] Figure 3 It is the result of algorithm clustering.

[0048] Figure 4 yes Figure 3 In-class samples and average results of category 1 in the clustering results;

[0049] (a) In-class samples of category 1, (b) averaged results.

[0050] Figure 5 yes Figure 3 In-class samples and average results of category 2 in the clustering results;

[0051] (a) In-class samples of category 2, (b) averaged results.

[0052] Figure 6 yes Figure 3 In-class samples and average results of category 3 in the clustering results;

[0053] (a) In-class samples of category 3, (b) averaged results.

[0054] Figure 7 yes Figure 3 In-class samples and average results of category 4 in the clustering results;

[0055] (a) In-class samples of category 4, (b) averaged results.

[0056] Figure 8 yes Figure 3 In-class samples and average results of category 5 in the clustering results;

[0057] (a) In-class samples of category 5, (b) averaged results. DETAILED DESCRIPTION

[0058] In order to make the process and implementation scheme of the present invention more intuitive, the technical solution of the present invention is further described below through specific examples and in conjunction with the accompanying drawings.

[0059] It should be clear that the clustering analysis algorithm provided by the present invention can have different performances in different instances, and based on various changes in the framework of the present invention, such as the magnetoencephalography device model, feature encoding method, clustering algorithm and weight coefficients involved described in the present invention, there may also be other alternatives or different combinations, but they do not deviate from the scope of the present invention in essence. In addition, the present invention does not describe some well-known algorithms and technical details. The examples and illustrations provided by the present invention are only responsible for intuitively explaining the algorithm flow and are not used to limit the present invention. The specific methods of the present invention are described in detail, such as Figure 1 shown.

[0060] For a 306-channel SQUID magnetoencephalography device, the number of magnetometer sensors is 102 and the number of gradiometer sensors is 204; the system sampling frequency is 1000 Hz. The spike waves generated by the interictal discharges of epilepsy in the sample are detected by an automated spike wave detection algorithm. There is a certain proportion of false positive samples or spike wave samples screened by doctors.

[0061] The time point corresponding to the GFP maximum typically represents the peak of the spike wave. This point typically has a high signal-to-noise ratio and is suitable for tracing the epileptogenic focus. Before clustering, the spike wave samples were aligned at the time point corresponding to the GFP maximum as the sample midpoint to facilitate subsequent cluster analysis.

[0062] The GFP maximum point of each sample is dipole fitted to obtain the goodness of fit GOF. The larger this indicator is, the better the fitting effect is. Therefore, a smaller threshold is first used to delete samples with a goodness of fit less than 0.3. Here, too large a threshold will result in too many deleted samples, making the distance matrix constructed during clustering more susceptible to noise samples. Single dipole fitting will obtain the 6 degree of freedom parameters of the dipole in space. For patients with focal epilepsy, their dipoles are usually concentrated near the lesion area or in symmetrical positions, and appear repeatedly. This prior will be used for subsequent abnormal point detection and removal operations.

[0063] Extract the topographic map corresponding to the time point of the maximum GFP value of each sample. The topographic map is extracted using the Python third-party package MNE to obtain the magnetometer topographic map and the gradiometer topographic map respectively. The topographic map is perceptually hashed as follows Figure 2As shown in Figure 2. The topographic map shape is resampled to 48×48 dimensions, and redundant information of the details is removed without losing the structural information of the topographic map. The resampled image is subjected to discrete cosine transform as shown in formula (3) to obtain the image frequency domain features. The discrete cosine transform compresses the image information in the upper left corner of the frequency domain, so the frequency domain features of size 8×8 in the upper left corner of the frequency domain are intercepted. This part of the features retains most of the original image information and can be used for further encoding operations. The mean of the intercepted features is calculated, and each feature is traversed. The features greater than the mean are set to 1, and vice versa. Finally, the topographic map of each sample is encoded as a 64-dimensional binary feature.

[0064] Clustering algorithms usually use a suitable distance metric to construct a distance matrix. Since the encoded features are binary features, the Hamming distance shown in formula (4) is used as the metric. The magnetometer distance matrix Matrix is ​​generated based on the encoded features of the magnetometer and gradiometer topographic maps. D_Mag and the gradiometer distance matrix Matrix D_Grad In order to comprehensively consider the contribution of magnetometer and gradiometer topography to the clustering results, the two distance matrices are fused to construct a total distance matrix as the final distance metric matrix, as shown in formula (5). In actual clinical analysis, magnetometer topography usually presents a specific dipole morphology, which is often used as an important basis for determining whether it is an epilepsy-related spike wave. Therefore, in this algorithm, α = 0.6, which makes the magnetometer distance matrix contribute more to the total distance matrix.

[0065] The resulting total distance matrix is ​​clustered using a hierarchical clustering algorithm. Hierarchical clustering ensures unique clustering results for large sample sizes. Clustering is determined using the silhouette coefficient method. By calculating the silhouette coefficient curve, the number of clusters corresponding to the local maximum or the greatest change in the silhouette coefficient is selected as the final number of clusters. This method limits the independent variable k of the silhouette coefficient curve to a range of [5, 25] and adjusts it based on the clustering results.

[0066] Under normal circumstances, clustering results are prone to outlier samples. Outlier samples are manifested in that, on the one hand, outlier samples are clustered individually or a very small number of outliers are clustered, and on the other hand, a small number of samples in the class are significantly different from other samples in the class, but belong to the same cluster. In view of the above analysis, the results obtained by the clustering algorithm of the present invention will delete clusters with less than 4 clustered samples, and use the LOF algorithm to remove outlier samples in the class for each retained cluster. The LOF algorithm also requires samples to construct a distance matrix. The distance matrix constructed here is constructed using the dipole position and direction coordinate information obtained by dipole fitting. First, the sample coordinate values ​​in the cluster are constructed using Euclidean distance to construct a distance matrix, and then the LOF algorithm is used to find the candidate abnormal samples. Secondly, the direction parameters of the samples in the class are used to construct a distance matrix using cosine distance, and the LOF algorithm is used to find the candidate abnormal samples. Finally, the union of the two results is taken as the final abnormal sample to be deleted. The parameter LOF-k in the LOF algorithm is selected as 4.

[0067] Finally, the samples with GOF values ​​less than 0.5 in the retained samples are deleted to obtain the final clustering results, as shown in Figure 3 shown. Figure 4-8 Indicates the results of categories 1-5 in the final clustering results, Figure 4 (a) shows five intra-class samples, where the first and second columns represent the magnetometer topography and gradiometer topography of the samples, respectively, and the third column represents the corresponding time series waveform of the samples; Figure 4 (b) The first row is the topographic map after the samples in the cluster are superimposed and averaged, and the second row is the temporal waveform after the samples are superimposed and averaged; Figures 5 to 8 Same thing.

[0068] In summary, the method of the present invention can achieve better clustering results for cluster analysis of epileptic spikes in MEG. This method is applicable to MEG spike datasets screened by physicians, as well as spike datasets containing a certain number of false-positive samples. The stacking and averaging of the clustering results of the present invention can improve the signal-to-noise ratio of spike signals, thereby improving the accuracy of tracing epileptogenic foci.

Claims

1. A method for cluster analysis of epileptic spikes based on magnetoencephalographic topography encoding, comprising the following steps: 1) for multiple spike wave samples between epileptic seizures of the same patient, calculating the global energy spectrum of each spike wave sample, and aligning the spike wave samples according to the time point corresponding to the maximum value in the global energy spectrum; 2) performing dipole fitting on each spike wave sample to obtain a dipole having spatial position parameters, direction parameters, and a goodness of fit (GOF); deleting spike wave samples whose GOF values ​​are less than a set first GOF threshold; 3) extracting the magnetoencephalogram topography corresponding to each maximum moment of the global energy spectrum; 4) performing perceptual hash coding on the magnetometer topography map in each MEG topography map extracted in step 3) to obtain a corresponding magnetometer topography coding feature; performing perceptual hash coding on the gradiometer topography map in each MEG topography map extracted in step 3) to obtain a corresponding gradiometer topography coding feature; 5) calculating a first distance matrix between spike samples using the magnetometer topography encoding features, calculating a second distance matrix between spike samples using the gradiometer topography encoding features, and fusing the obtained first and second distance matrices; 6) Clustering the fused distance matrix and deleting outliers within the class and spike samples with GOF values ​​within the class less than the set second GOF threshold; 7) For each class processed in step 6), if the number of spike wave samples in the class is greater than the set threshold H, the class is retained; then the spike wave samples in each retained class are superimposed and averaged to obtain the superimposed and averaged spike wave.

2. The method according to claim 1, characterized in that The fused distance matrix is ​​clustered using a hierarchical clustering algorithm. For the main cluster obtained by clustering, the local outlier factor algorithm is used to remove outliers in the main cluster.

3. The method according to claim 2, characterized in that The distance matrix required by the local anomaly factor algorithm is calculated using the direction parameters and spatial position parameters of the dipole.

4. The method according to claim 2, characterized in that The silhouette coefficient method is used to determine the number of cluster categories.

5. The method according to claim 1, characterized in that In step 5), use the formula Calculate the first distance matrix and the second distance matrix; wherein, l i ,s i represents the i-th eigenvalue of the eigenvector, wherein the feature is a magnetometer topographic map encoding feature or a gradiometer topographic map encoding feature; Is an exclusive OR operation; using the formula Matrix D =aMatrix D_Mag +(1-a)Matrix D_Grad The first distance matrix and the second distance matrix are fused; Matrix D Represents the distance matrix after fusion; Matrix D_Mag Represents the first distance matrix calculated based on the encoded features of the magnetometer topographic map; Matrix D_Grad represents the second distance matrix calculated based on the coded features of the gradiometer topographic map; α is the weight coefficient.

6. The method according to claim 1, characterized in that The second GOF threshold is greater than the first GOF threshold.

7. The method according to claim 6, characterized in that The second GOF threshold is 0.5; the first GOF threshold is 0.

3.

8. The method according to any one of claims 1 to 7, characterized in that: The perceptual hash coding process is as follows: a) extracting the magnetometer topography map and the gradiometer topography map corresponding to the time point at the maximum value of the global energy spectrum and removing details and high-frequency information therein; b) performing a discrete cosine transform on the image processed in step a); c) intercepting the low-frequency portion of the transform domain of the magnetometer topography map after the discrete cosine transform and calculating its mean to obtain the corresponding magnetometer topography map encoding feature; The low-frequency part of the transform domain of the gradiometer topography map after discrete cosine transform is intercepted and its mean is calculated to obtain the corresponding gradiometer topography encoding feature.

9. A server, characterized in that: The method comprises a memory and a processor, wherein the memory stores a computer program, the computer program is configured to be executed by the processor, and the computer program comprises instructions for executing the steps of the method according to any one of claims 1 to 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 8 are implemented.

Images