Quantitative analysis method for multi-parameter transcranial magnetic stimulation

Abstract

The invention relates to a multi-parameter transcranial magnetic stimulation quantitative analysis method, which comprises the following steps of: acquiring neuron activity index data, and performing standardization processing on the neuron activity index data to generate a first data set F1; identifying the acquired data by using a DBSCAN clustering algorithm and removing noise points to obtain a second data set F2 only containing core points and boundary points; and S3, calculating a distance metric and a link distance based on the second data set F2 obtained in the step S2, applying a hierarchical clustering algorithm to obtain a clustering tree diagram, and determining a classification threshold and performing classification according to the form of the clustering tree diagram and the requirements of a five-level quantitative classification system. According to the quantitative analysis method, the parameters such as the stimulation intensity, the frequency and the pulse mode are quantitatively graded, so that a clinician can more accurately select a proper treatment scheme according to the specific conditions such as the age, the gender and the disease state of a patient.

Description

Technical Field

[0001] The present invention relates to the field of medical technology, and more specifically, to a quantitative analysis method for multi-parameter transcranial magnetic stimulation (TMS). Background Art

[0002] Transcranial magnetic stimulation (TMS), as a non-invasive neuromodulation technique, has attracted extensive attention in the fields of neuroscience and clinical rehabilitation. By applying external stimuli to the brain, it can regulate brain nerve activities, thereby improving the symptoms of various neurological diseases or enhancing the cognitive abilities of healthy people, etc. TMS generates a magnetic field outside the scalp and induces an induced current in the skull, thus affecting the activities of neurons. Although certain achievements have been made in clinical applications, there are still significant uncertainties in the treatment effects. This uncertainty mainly stems from individual differences and the lack of precise quantification of stimulus intensity and effects. Factors such as the brain anatomical structure, neurophysiological characteristics, and disease states of different individuals vary, and these differences will significantly affect the treatment effects of TMS. For example, differences in cerebral cortex thickness, neuronal excitability, and neurotransmitter levels may all lead to different responses of different individuals to the same stimulus parameters. Therefore, in order to improve the accuracy and effectiveness of treatment, it is crucial to establish a scientific and reasonable quantitative grading system. Summary of the Invention

[0003] In view of the above technical problems, the present invention proposes a quantitative analysis method for multi-parameter transcranial magnetic stimulation (TMS), specifically including: Step 1, obtaining neuron activity index data and performing standardization processing on it to generate a first data set F1; Step 2, using the DBSCAN clustering algorithm to remove noise points from the first data set F1 to obtain a second data set F2 containing only core points and boundary points; Step 3, based on the second data set F2 obtained in Step 2, calculating the distance metric and the linkage distance, and applying the hierarchical clustering algorithm to obtain a clustering dendrogram, and combining the morphology of the clustering dendrogram, determining the grading threshold and performing grading according to the requirements of the five-level quantitative grading system.

[0004] Further, Step 2 specifically includes:[[]]

[0005] Step 21, determining the preliminary range of clustering algorithm parameters, specifically including:[[]]

[0006] Step 211, determining the preliminary range of the neighborhood radius ∈, specifically including;[[]]

[0007] Calculating the average reachable distance ARD(∈) of data points under different ∈ values:[[]]

[0008] Wherein,[[]]

[0009] N∈ (x i ) represents the set of data points within the ∈-neighborhood of the data point x i . |N ∈ (x i )| represents the number of data points within the ∈-neighborhood of the data point x i , and d(x i , x j ) represents the Euclidean distance between the data point x i and x j : where N represents the number of data points in the dataset, x ik and x jk are their k-th eigenvalues respectively, and D is the dimension of the data;

[0010] Plot the curve of ARD(∈) versus ∈, find the elbow position of the curve, and determine the preliminary range of ∈;

[0011] Step 212, determine the initial range of the MinPts value, 1*D ≤ MinPts ≤ 10*D;

[0012] Step 22, clustering classification and noise removal, specifically including:

[0013] Step 221, for each data point in the first dataset F1, calculate the number of data points within the neighborhood of ∈ determined in Step 2123. If the number is greater than or equal to MinPts, mark this point as a core point and create a new cluster, and add this point and all data points within its ∈-neighborhood to the cluster;

[0014] Step 222, for each core point in the cluster, repeat the above process, and add the unmarked data points within the ∈-neighborhood of this point to the cluster until the cluster no longer expands;

[0015] Step 223, if the number of data points within the ∈-neighborhood of a data point is less than MinPts, temporarily mark this point as a noise point;

[0016] Step 224, after completing the clustering, remove the data marked as noise points from the first dataset F1;

[0017] Step 23, determine the final values of the neighborhood radius ∈ and the MinPts value, specifically including;

[0018] Step 2321, for each ∈ value, select different MinPts values for experimentation;

[0019] Step 2322, calculate the silhouette coefficient of the corresponding clustering result

[0020] Wherein, a(i) is the average distance from data point i to other data points within its affiliated cluster, and b(i) is the minimum average distance from data point i to data points in other clusters;

[0021] Step 2323, select the combination of ∈ and MinPts that maximizes the silhouette coefficient s(i) as the final neighborhood radius ∈ and the value of MinPts.

[0022] Step 24, perform the clustering classification and noise removal in Step 22 according to the finally determined neighborhood radius ∈ and the value of MinPts, and obtain the second data set F2 containing only core points and boundary points.

[0023] Further, the specific steps of Step 3 are as follows:

[0024] Step 31, calculate the distance metric and the linkage distance, specifically including:

[0025] Step 311, calculate the distance metric. For two data points i and j in the second data set F2, calculate the Euclidean distance between them:

[0026] Where x ik and x jk are their k-th eigenvalues respectively, and D is the dimension of the data;

[0027] Step 312, calculate the average linkage distance, specifically including:

[0028] When calculating the distance between clusters, the average linkage distance method is adopted:

[0029] Where |A| and |B| are the numbers of data points in the new clusters A and B respectively, and d(i, j) is the Euclidean distance between data points i and j;

[0030] Step 32, clustering process and dendrogram generation, specifically including:

[0031] Step 321, initialization and distance matrix calculation: Initially, regard each data point in the second data set F2 as a separate cluster, calculate the Euclidean distance between every two clusters, and form an M×M distance matrix, where M is the number of data points in the second data set F;

[0032] Step 322, cluster merging and dendrogram updating, specifically including:

[0033] Step 3221, find the two clusters A and B with the minimum distance in the distance matrix, merge the clusters A and B into a new cluster, update the distance matrix, delete the rows and columns related to A and B, and calculate the distance between the new cluster and other clusters;

[0034] Step 3222: Repeat Step 3221. When all data points are merged into one cluster or the preset stop condition is reached, the merging process ends, and a clustering dendrogram is gradually constructed, recording the inter-cluster distance and the corresponding clustering situation at each merge.

[0035] Step 33: Determination of hierarchical threshold and grade division, specifically including:

[0036] Step 331: Determine the hierarchical threshold based on the generated clustering dendrogram above, specifically including: using the method of equally dividing the height of the dendrogram to determine the hierarchical threshold. Assume the range of the inter-cluster distance in the dendrogram is [D min , D max , divide it into five equal segments, and the threshold points are D1, D2, D3, D4 respectively.

[0037] Step 332: According to the determined hierarchical threshold points, divide the data set into five grades, corresponding to weak-intensity stimulation, low-intensity stimulation, medium-intensity stimulation, high-intensity stimulation, and ultra-high-intensity stimulation respectively.

[0038] Furthermore, the eigenvalue in Step 3 includes: relative change rate of membrane potential, relative change rate of action potential amplitude, relative change rate of frequency, and relative reduction rate of discharge threshold.

[0039] The present invention proposes a quantitative analysis method for multi-parameter transcranial magnetic stimulation, which combines the advantages of hierarchical clustering algorithm and DBSCAN clustering algorithm to construct a quantitative grading system, and can provide clear guidance for parameter selection in TMS treatment. By quantitatively grading parameters such as stimulation intensity, frequency, and pulse pattern (for TMS), clinicians can more accurately select a suitable treatment plan according to the specific conditions of patients, such as age, gender, disease status, etc. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] The drawings described herein are used to provide a further understanding of the present disclosure, and constitute a part of the present disclosure. The schematic embodiments and descriptions thereof are used to explain the present disclosure and do not constitute an improper limitation to the present disclosure. In the drawings:

[0041] Figure 1 shows a schematic flow chart of the quantitative analysis method for multi-parameter transcranial magnetic stimulation of the present invention;

[0042] Figure 2 shows a schematic diagram of the relative change rate of membrane potential varying with the stimulation intensity range in an embodiment of the quantitative analysis method for multi-parameter transcranial magnetic stimulation of the present invention;

[0043] Figure 3Shows a schematic diagram of the relative change rate of membrane potential varying with the stimulation frequency range in an embodiment of the quantitative analysis method for multi-parameter transcranial magnetic stimulation of the present invention;

[0044] Figure 4 Shows a schematic diagram of the relative change rate of action potential amplitude varying with the intensity range in an embodiment of the quantitative analysis method for multi-parameter transcranial magnetic stimulation of the present invention;

[0045] Figure 5 Shows a schematic diagram of the relative change rate of action potential amplitude varying with the stimulation frequency range in an embodiment of the quantitative analysis method for multi-parameter transcranial magnetic stimulation of the present invention;

[0046] Figure 6 Shows a schematic diagram of the relative change rate of frequency varying with the intensity range in an embodiment of the quantitative analysis method for multi-parameter transcranial magnetic stimulation of the present invention;

[0047] Figure 7 Shows a schematic diagram of the relative change rate of frequency varying with the stimulation frequency range in an embodiment of the quantitative analysis method for multi-parameter transcranial magnetic stimulation of the present invention;

[0048] Figure 8 Shows a schematic diagram of the relative reduction rate of discharge threshold varying with the intensity range in an embodiment of the quantitative analysis method for multi-parameter transcranial magnetic stimulation of the present invention;

[0049] Figure 9 Shows a schematic diagram of the relative reduction rate of discharge threshold varying with the stimulation frequency range in an embodiment of the quantitative analysis method for multi-parameter transcranial magnetic stimulation of the present invention. Detailed implementation manners

[0050] To make the objectives, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are some but not all of the embodiments of the present invention.

[0051] Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the claimed present invention, but merely represents selected embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts fall within the scope of protection of the present invention.

[0052] The data of transcranial magnetic stimulation intensity quantification grading involves multiple neuron activity indexes. In the present invention, the relative change rate of membrane potential, the relative change rate of action potential amplitude, the relative change rate of frequency, and the relative reduction rate of neuron firing threshold are selected as analysis indexes. These indexes reflect the response of neurons to stimulation from different perspectives, and they are interrelated and jointly affect the judgment of stimulation intensity levels. For example, the change of membrane potential may be related to the generation and propagation of action potentials, and the change of firing threshold will affect the excitability of neurons, and thus be associated with frequency changes. Therefore, when selecting a clustering method, an algorithm that can handle multivariate data needs to be considered to make full use of the information contained in these indexes and accurately divide the stimulation intensity levels.

[0053] In addition, since the indexes are not completely independent of each other and there is a certain correlation. For example, the change of membrane potential may lead to corresponding changes in action potential amplitude and frequency, and the decrease of firing threshold may make neurons more likely to generate action potentials, thus affecting indexes such as frequency. This correlation makes the structure of the data more complex, and traditional univariate analysis methods or simple clustering algorithms may not be able to effectively capture the internal patterns in the data. Therefore, the selected clustering method should be able to handle the correlation between variables, avoid bias in the clustering results caused by repeated consideration of relevant information, and at the same time be able to uncover the relationship between the unique neural activity patterns represented by different index combinations and the stimulation intensity.

[0054] As the stimulation intensity changes, the changes of neuron activity indexes do not show a simple linear relationship. At low-intensity stimulation, the indexes may change relatively slowly, while when the stimulation intensity reaches a certain level, the indexes may change sharply, showing a non-linear growth or fluctuation trend. For example, the relative change rate of membrane potential may only increase slightly under low-intensity stimulation, but rise rapidly under high-intensity stimulation. This non-linear distribution makes the distribution of data points in the feature space uneven, and traditional clustering methods based on linear models are difficult to accurately describe and divide the data under different stimulation intensity levels, and a clustering algorithm that can adapt to non-linear data distribution needs to be selected.

[0055] In addition, due to the difference in the density of data in different stimulation intensity regions. In some stimulation intensity ranges, there may be more data points clustering, indicating that the neuron activity indexes change relatively concentrated in this range and have a high density; while in other intensity ranges, the data points may be relatively sparse and the changes are relatively scattered. For example, there may be more data points near medium-intensity stimulation because this intensity range is often used in treatment research and the experimental data is relatively rich, while the data points in the weak-intensity or ultra-intensity stimulation regions may be fewer. This density non-uniformity requires the clustering method to be able to identify the data patterns in different density regions, avoid mismerging or ignoring the data in low-density regions, and ensure accurate division of stimulation intensity levels.

[0056] The DBSCAN algorithm is highly robust to noise and outliers, which is an important advantage in the clustering of transcranial magnetic stimulation data. By defining core points, border points, and noise points, it can effectively identify and exclude noisy data. In a dataset, even if there are a small number of outliers caused by experimental interference or individual special responses, the DBSCAN algorithm can correctly label them as noise points without affecting the clustering process of core points and border points, thus ensuring that the clustering results can accurately reflect the relationship between the true neuron activity patterns and the stimulation intensity represented by most data points. For example, for those abnormal data points that significantly deviate from the normal data distribution, the DBSCAN algorithm will not force them into a certain cluster but handle them separately, ensuring the accuracy and stability of clustering.

[0057] Given the non - linear distribution and density inhomogeneity of transcranial magnetic stimulation data, the advantages of the DBSCAN algorithm are fully demonstrated. It can discover clusters of any shape, without being restricted by the assumptions about the shape of clusters (such as spherical clusters) in traditional clustering methods. When the neuron activity index shows a complex - shaped distribution with the change of stimulation intensity, the DBSCAN algorithm can accurately identify the natural data cluster structure. At the same time, its density - based clustering method can well adapt to the density inhomogeneity of data, forming tight clusters in high - density regions and reasonably dividing clusters or identifying isolated data points in low - density regions. For example, in the high - density region of medium - intensity stimulation data, the DBSCAN algorithm can effectively identify the cluster composed of core points and border points, and in the low - density regions of weak - intensity or super - intensity stimulation, it can also correctly handle the attribution of data points, providing strong support for accurately classifying the stimulation intensity levels.

[0058] In addition, when dealing with multi - variable data, the DBSCAN algorithm mainly clusters based on the density of data points in the feature space and makes relatively insufficient use of the correlations between variables. Although it can consider the comprehensive density of multi - variable data to a certain extent, it may not be able to fully capture the complex non - linear correlations between variables. To make up for this deficiency, before applying the DBSCAN algorithm, dimensionality reduction processing (such as principal component analysis, etc.) can be performed on multi - variable data, converting the original correlated variables into a few uncorrelated principal components, reducing the data dimension while retaining the main information, so that the DBSCAN algorithm can more effectively process the dimensionality - reduced data and improve the clustering effect. Or combine it with other algorithms that can better handle multi - variable correlations (such as hierarchical clustering algorithm) to give play to their respective advantages and achieve more accurate quantification and classification of stimulation intensity.

[0059] On the other hand, hierarchical clustering algorithms can simultaneously consider the distance relationships between data points with multiple variables. It measures similarity by calculating the distances between data points in a multi-variable feature space (such as Euclidean distance, etc.), thereby enabling a comprehensive analysis of information on multiple metrics such as membrane potential, action potential, frequency, and firing threshold. Regarding the correlations between metrics, hierarchical clustering algorithms will naturally take into account the co-variation between variables when calculating distances, and will not treat each variable separately, thus being able to better capture the inherent patterns in the data. For example, when the membrane potential and the amplitude of the action potential change simultaneously, the hierarchical clustering algorithm can perform clustering based on their combined distances in the multi-variable space, avoiding ignoring the association between them due to considering variables separately.

[0060] Although the distance-based calculation of the hierarchical clustering algorithm itself may have certain limitations when dealing with highly non-linear data, it can, to a certain extent, adapt to the complexity of the data distribution. Its agglomerative or divisive clustering process can gradually reveal the hierarchical structure of the data. In the case of uneven data density, it can reflect the differences in different density regions through clustering mergers or splits at different levels. For example, in regions with higher density, data points are more likely to be merged together in the early clustering stage, while data points in regions with lower density are merged in subsequent stages, thus adapting to density unevenness to a certain extent. However, for extremely non-linear and complex density change situations, the standalone hierarchical clustering algorithm may need further improvement or combination with other methods.

[0061] Hierarchical clustering algorithms are sensitive to noise and outliers, which is one of their main drawbacks. Outliers can greatly affect the calculation of distances between data points, resulting in distortion of the structure of the clustering dendrogram, and further affecting the determination of the hierarchical threshold and the final hierarchical result. To overcome this problem, in this invention, it is chosen to preprocess the data before applying the hierarchical clustering algorithm, such as using data cleaning techniques to identify and remove obvious outliers, or using robust distance metrics (such as median distance, etc.) to reduce the impact of outliers. In addition, combining other algorithms that are robust to noise (such as the DBSCAN algorithm) can further improve the overall clustering effect. First, use the DBSCAN algorithm to remove noise points, and then apply the hierarchical clustering algorithm on the clean data for more refined classification.

[0062] Therefore, the hierarchical clustering algorithm and the DBSCAN algorithm have complementary characteristics. The hierarchical clustering algorithm is good at presenting the hierarchical structure of data, can provide intuitive information about the progressive relationship between different stimulus intensity levels, and has certain advantages in dealing with multivariate correlations, but is sensitive to noise. The DBSCAN algorithm is robust to noise, can discover clusters of arbitrary shapes and adapt to density changes, but is relatively weak in dealing with multivariate correlations and presenting hierarchical structures. By combining the two, the present invention can first use the DBSCAN algorithm to remove noise points to obtain relatively pure data, and then apply the hierarchical clustering algorithm to these data, which can not only avoid the interference of noise on hierarchical clustering, but also give full play to the advantages of hierarchical clustering in revealing the data hierarchical structure and dealing with multivariate correlations, so as to achieve complementary advantages and improve the accuracy of clustering and grading.

[0063] As Figure 1 shown, a quantitative analysis method for multi-parameter transcranial magnetic stimulation provided by a specific embodiment of the present invention specifically includes:

[0064] Step 1, obtain neuron activity index data and perform normalization processing on it to generate the first data set F1.

[0065] In one embodiment, the present invention uses SimNIBS software to perform macroscopic field simulation of transcranial electromagnetic fields. For transcranial magnetic stimulation (TMS), select a suitable coil type, such as a circular coil, figure-eight coil, etc. The present invention uses a MagVenture_MC-B70 Magstim 70mm coil. Based on the brain anatomical structure and functional regions, determine the coil placement according to the 10-10 system, select the dorsolateral prefrontal cortex of F3, and place the coil parallel to the scalp surface. The present invention selects the Aberra_L5_model neuron cell model, which is constructed based on the detailed morphology and electrophysiological characteristics of actual cerebral cortex neurons, can accurately simulate the real response of neurons and effectively respond to TMS stimulation.

[0066] Use the created Aberra_L5_model neuron model for TMS stimulation simulation experiments. Among them, for the frequency parameter, set multiple different frequency values, starting from 0.1 Hz, increasing by 0.1 Hz each time, such as 0.1 Hz, 0.2 Hz, 0.3 Hz, etc. When the frequency is greater than 1 Hz, increase by 1 Hz each time, such as 2 Hz, 3 Hz, 4 Hz, etc. For the intensity parameter, start from 10% RMT (relative intensity), increase by 10% RMT each time until 150% RMT, and increase by 75%, 85%, 95%, 105%, 115%, 125% within the common range. In each group of experiments, keep other parameters consistent, such as the pulse width is fixed at 0.2 ms, etc.

[0067] After the TMS stimulation starts, the changes in the neuronal membrane potential are monitored in real time. Under different TMS parameter settings of frequency and intensity, the membrane potential shows different changing trends. By plotting the curve of the membrane potential changing with time and comparing the curve differences between different frequency and intensity groups, the influence rules of TMS frequency and intensity on the changes in the neuronal membrane potential are analyzed.

[0068] Under TMS stimulation, the firing threshold of neurons may change. For example, when the TMS stimulation intensity increases, the firing threshold of neurons may decrease, which means that neurons are more likely to generate action potentials, that is, the excitability of neurons increases. By measuring the changes in the neuronal firing threshold, the influence of TMS on the excitability of neurons can be understood.

[0069] The changes in the action potential amplitude may be related to factors such as the ion channel function of neurons, the intracellular ion concentration, and the electrical properties of the cell membrane. Under different stimulation conditions, the changes in the action potential amplitude can reflect the functional state of neurons and their adaptability to stimuli.

[0070] Obtain data such as the relative change rate of the membrane potential, the relative change rate of the membrane potential frequency, the relative change rate of the action potential amplitude, and the relative reduction rate of the firing threshold under the above TMS stimulation. These indicators reflect the response of neurons to stimuli from different angles, and they are interrelated and jointly affect the judgment of the stimulus intensity level. For example, the change in the membrane potential may be related to the generation and propagation of action potentials, and the change in the firing threshold will affect the excitability of neurons, and thus be related to the frequency change.

[0071] In order to fully analyze the action mechanism of TMS and provide a theoretical basis for optimizing the treatment plan, the present invention selects data such as the relative change rate of the membrane potential, the relative change rate of the membrane potential frequency, the relative change rate of the action potential amplitude, and the relative reduction rate of the firing threshold under TMS stimulation as the target data for analysis.

[0072] It should be noted that the above indicators are not completely independent of each other and there is a certain correlation. For example, the change in the membrane potential may cause corresponding changes in the action potential amplitude and frequency, and the decrease in the firing threshold may make neurons more likely to generate action potentials, thus affecting indicators such as frequency. This correlation makes the structure of the data more complex, and traditional univariate analysis methods or simple clustering algorithms may not be able to effectively capture the internal patterns in the data. Therefore, the selected clustering method should be able to handle the correlation between variables, avoid bias in the clustering results due to repeated consideration of relevant information, and at the same time be able to uncover the relationship between the unique neural activity patterns represented by different combinations of indicators and the stimulus intensity.

[0073] In one embodiment, the neuron activity metrics (such as membrane potential, action potential amplitude, frequency, and firing threshold changes) are normalized, and the normalized data is output as a feature vector to obtain the first dataset F1.

[0074] In one embodiment, the normalization can be performed by methods such as mean normalization, Z - score normalization, and min - max normalization. The present invention does not make specific limitations.

[0075] Step 2: Use the DBSCAN clustering algorithm to remove noise points from the first dataset F1, obtaining a second dataset F2 that only contains core points and border points.

[0076] Through comparison tests, the inventors found that the DBSCAN algorithm has strong robustness to noise and outliers in transcranial magnetic stimulation data clustering. By defining core points, border points, and noise points, it can effectively identify and exclude noisy data. In the dataset, even if there are a small number of outliers caused by experimental interference or individual special reactions, the DBSCAN algorithm can correctly label them as noise points without affecting the clustering process of core points and border points, thus ensuring that the clustering result can accurately reflect the relationship between the true neuron activity pattern represented by most data points and the stimulation intensity.

[0077] In a preferred embodiment, for the neuron activity metric data obtained in step 1, the present invention selects to use the DBSCAN clustering algorithm to identify and remove noise points from the obtained data, obtaining a dataset that only contains core points and border points.

[0078] The DBSCAN clustering algorithm identifies and removes noise points from the obtained data, obtaining a dataset that only contains core points and border points, specifically including:

[0079] Step 21: Determine the preliminary range of the clustering algorithm parameters, specifically including:

[0080] Step 211: Determine the preliminary range of the neighborhood radius ∈, specifically including;

[0081] Calculate the average reachability distance ARD(∈) of data points under different ∈ values:

[0082] Where

[0083] N ∈ (x i ) represents the set of data points within the ∈ - neighborhood of data point x i , |N ∈ (x i )| represents the number of data points within the ∈ - neighborhood of data point x iThe number of data points within the ∈-neighborhood, d(x i , x j ) represents the Euclidean distance between data points x i and x j : where N represents the number of data points in the dataset, x ik and x jk are their k-th eigenvalues respectively, and D is the dimension of the data.

[0084] Plot the curve of ARD(∈) versus ∈, find the elbow position of the curve, and use the ∈ value corresponding to this elbow position as the initial ∈ value. The elbow position refers to the turning point where the curve changes from a rapid descent to a slow descent, which can be found by observing the change in the curve slope. Different ∈ values, as the neighborhood radius parameter in the DBSCAN algorithm, will affect the neighborhood range of data points and thus have different effects on the clustering result.

[0085] Step 212, Determine the initial range of the MinPts value, 1*D ≤ MinPts ≤ 10*D.

[0086] Step 22, Clustering classification and noise removal, specifically including:

[0087] Step 221, For each data point in the first dataset F1, calculate the number of data points within the neighborhood of ∈ determined in step 2123. If the number is greater than or equal to MinPts, mark this point as a core point and create a new cluster, and add this point and all data points within its ∈-neighborhood to the cluster;

[0088] Step 222, For each core point in the cluster, repeat the above process, adding the unmarked data points within the ∈-neighborhood of this point to the cluster until the cluster no longer expands;

[0089] Step 223, If the number of data points within the ∈-neighborhood of a data point is less than MinPts, temporarily mark this point as a noise point;

[0090] Step 224, After completing the clustering, remove the data marked as noise points from the first dataset F1;

[0091] Step 23, Determine the final values of the neighborhood radius ∈ and the MinPts value, specifically including;

[0092] Step 2321, For each ∈ value, select different MinPts values for experimentation;

[0093] Step 2322, Calculate the silhouette coefficient of the corresponding clustering result

[0094] Among them, a(i) is the average distance from data point i to other data points within its affiliated cluster, and b(i) is the minimum average distance from data point i to data points in other clusters;

[0095] Step 2323, select the combination of ∈ and MinPts that maximizes the silhouette coefficient s(i) as the final neighborhood radius ∈ and the value of MinPts.

[0096] Step 24, perform the clustering classification and noise removal in Step 22 according to the finally determined neighborhood radius ∈ and the value of MinPts, and obtain the second data set F2 that only contains core points and border points.

[0097] After clustering by the DBSCAN algorithm, data points are divided into core points, border points, and noise points. A core point is a point whose number of data points within its ∈-neighborhood is greater than or equal to MinPts; a border point is a point within the neighborhood of a core point, but the number of data points within its own neighborhood is less than MinPts. The algorithm mainly focuses on core points and border points to define clusters. Noise points do not belong to any meaningful clusters, so the data set after removing noise points only contains core points and border points.

[0098] Step 3, based on the second data set F2 obtained in Step 2, calculate the distance metric and the linkage distance, and apply the hierarchical clustering algorithm to obtain a clustering dendrogram. Combining the morphology of the clustering dendrogram, according to the requirements of the five-level quantitative grading system, determine the grading threshold and perform grading, specifically including:

[0099] Step 31, calculate the distance metric and the linkage distance, specifically including:

[0100] Step 311, calculate the distance metric. For two data points i and j in the second data set F2, calculate the Euclidean distance between them:

[0101] where, x ik and x jk are their k-th eigenvalue respectively, and D is the dimension of the data. In a preferred embodiment, the dimension D is 4. Specifically, the 4-dimensional eigenvalue includes the relative change rate of membrane potential, the relative change rate of action potential amplitude, the relative change rate of frequency, and the relative reduction rate of discharge threshold.

[0102] Step 312, calculate the average linkage distance, specifically including:

[0103] In the initial stage of hierarchical clustering, each data point in the second data set F2 is regarded as a separate cluster. When calculating the distance between clusters, the average linkage distance method can be used:

[0104] where, |A| and |B| are the number of data points in the new clusters A and B respectively, d(xi , x j ) is the data point x i and x j The Euclidean distance between them, where the "new cluster" is formed by gradually merging the initial single data point clusters as the algorithm runs during the hierarchical clustering process.

[0105] Step 32, Clustering process and dendrogram generation, specifically including:

[0106] Step 321, Initialization and distance matrix calculation. Initially, each data point in the second data set F2 is regarded as a separate cluster, and the Euclidean distance between pairwise clusters is calculated to form an M×M distance matrix. Here, M is the number of data points in the second data set F2, which is different from the total number N of data points in the first data set F1 before being processed by the DBSCAN algorithm because the number of data points may decrease after removing noise points by the DBSCAN algorithm.

[0107] Step 322, Cluster merging and dendrogram update, specifically including:

[0108] Step 3221, Find the two clusters A and B with the smallest distance in the distance matrix, merge the cluster A and cluster B into a new cluster, update the distance matrix, delete the rows and columns related to A and B, and calculate the distance between the new cluster and other clusters;

[0109] Step 3222, Repeat step 3221. When all data points are merged into one cluster, or when a pre-set stop condition is reached, the merging process ends, and the inter-cluster distance and the corresponding clustering situation at each merge are recorded. The specific recording method can be to use a list, where each element is a tuple containing the numbers of the two merged clusters and their distance, thereby gradually constructing a clustering dendrogram.

[0110] In one embodiment, the pre-set stop condition can be that the number of clusters reaches a specified value, or the inter-cluster distance exceeds a certain threshold, etc.

[0111] Step 33, Hierarchical threshold determination and level division, specifically including:

[0112] Step 331, Determine the hierarchical threshold based on the generated clustering dendrogram, specifically including: Using the method of equally dividing the height of the dendrogram to determine the hierarchical threshold: Assume the range of the inter-cluster distance in the dendrogram is [D min , D max , divide it into five equal segments, and the threshold points are D1, D2, D3, D4;

[0113] Step 332: Divide the dataset into five levels according to the determined grading threshold points, corresponding to weak-intensity stimulation, low-intensity stimulation, medium-intensity stimulation, high-intensity stimulation, and ultra-high-intensity stimulation respectively.

[0114] Specifically, the specific definitions of these five levels are as follows:

[0115] 1) Level 1: Weak-intensity stimulation

[0116] Definition: The stimulation intensity is low, and it has almost no observable impact on various physiological indicators of neurons. There are no obvious changes in indicators such as the firing threshold and frequency of neurons, membrane potential changes, and action potential amplitude.

[0117] Applicable situations: Generally used for setting control groups in basic research, or when conducting initial exploratory stimulation on certain special sensitive populations to evaluate their basic tolerance to stimulation.

[0118] Under weak-intensity stimulation, the stimulation intensity range (RMT%) is less than 75%, and the stimulation frequency range (Hz) is less than 0.5.

[0119] 2) Level 2: Low-intensity stimulation

[0120] Definition: At this level, the changes in indicators such as the firing threshold and frequency of neurons, membrane potential changes, and action potential amplitude are relatively small. It usually corresponds to a stimulation intensity lower than a certain threshold, which can be calculated based on simulation results. Specifically, the firing frequency of neurons may only increase slightly, the amplitude of membrane potential changes is small, the change in action potential amplitude is not significant, and the degree of neuronal synchronization increases slightly but not significantly.

[0121] Applicable situations: Suitable for mild neuromodulation or exploratory treatment. For example, in some preliminary studies, when wanting to understand the basic impact of TMS on a certain area of the brain, low-intensity stimulation can be used. Or when treating some patients with relatively mild conditions and relatively stable symptoms, low-intensity stimulation can be used as a preliminary intervention method to observe the patients' reactions and avoid adverse reactions that may be caused by excessive stimulation. It is also commonly used in the early treatment stage of some patients with poor tolerance to stimulation to provide a reference for adjusting subsequent treatment plans.

[0122] Under low-intensity stimulation, the stimulation intensity range (RMT%) is greater than or equal to 75% and less than 95%, and the stimulation frequency range (Hz) is greater than or equal to 0.5 and less than 1.

[0123] 3) Level 3: Medium-intensity stimulation

[0124] Definition: The changes in neuron response indicators are relatively obvious. There is a moderate increase or change in indicators such as neuron firing frequency, membrane potential change, and action potential amplitude. The degree of neuron synchronization may increase to a certain extent. This corresponds to a certain range of stimulus intensity, which is also determined according to the grading standard. Generally, the neuron firing frequency increases significantly compared to the baseline level, the change amplitude of the membrane potential is moderate, and the action potential amplitude also changes to a certain extent. The enhanced synchronous activity between neurons makes the transmission and integration of nerve signals more efficient.

[0125] Applicable situations: It is used for the treatment of moderate nerve dysfunction. For example, in some patients with moderate cognitive impairment or motor dysfunction, medium-intensity stimulation can promote neuroplasticity changes by regulating neuron activity and improve the symptoms of patients. For some patients with chronic neurological diseases who still have certain functional impairments during the stable stage of the disease, medium-intensity stimulation can be used as a routine treatment option to maintain or further improve nerve function and improve the quality of life of patients.

[0126] Under medium-intensity stimulation, the stimulus intensity range (RMT%) is greater than or equal to 95% and less than 110%, and the stimulus frequency range (Hz) is greater than or equal to 1 and less than 5.

[0127] 4) Grade 4: High-intensity stimulation

[0128] Definition: It will cause relatively significant changes in indicators such as neuron firing frequency, membrane potential change, and action potential amplitude. The degree of neuron synchronization increases to a large extent, and it may induce neuroplasticity changes similar to long-term potentiation (LTP) or long-term depression (LTD). This corresponds to a stimulus intensity higher than a certain conventional threshold but still within the clinically acceptable and relatively safe range. At this time, the neuron firing frequency increases significantly, the membrane potential changes significantly, the change in action potential amplitude is relatively prominent, and the neurons are highly synchronized. This intensity of stimulation can have a strong remodeling effect on the brain neural circuit.

[0129] Applicable situations: It is used for the treatment of relatively severe neurological diseases or to promote the further improvement of nerve function. For example, in some patients with moderate to severe depression, certain stages of Parkinson's disease, or stroke sequelae with obvious functional impairments who respond poorly to conventional treatments, high-intensity stimulation can be used as an intensive treatment method to strongly regulate neuron activity in order to achieve more significant nerve function recovery and reconstruction effects. However, during the use process, it is necessary to closely monitor various physiological indicators and reactions of patients to ensure the safety and effectiveness of the treatment, and adjust the stimulation parameters in a timely manner according to the specific situation of patients.

[0130] Under high-intensity stimulation, the stimulation intensity range (RMT%) is greater than or equal to 110% and less than 120%, and the stimulation frequency range (Hz) is greater than or equal to 5 and less than 15.

[0131] 5) Level 5: Ultra-high-intensity stimulation

[0132] Definition: The stimulation intensity exceeds the conventional use by a large margin and is usually only applied with caution in specific situations. This intensity of stimulation will cause very significant changes in indicators such as the discharge frequency of neurons, the change in membrane potential, and the amplitude of action potentials, and the degree of neuronal synchronization is extremely high. Although it is still within the safe range under strict control and specific conditions, there are certain risks and it requires great caution.

[0133] Applicable situations: Generally used for special scientific research exploration or extremely special clinical cases. In scientific research, for example, in some frontier neuroscience research, in order to deeply explore the potential response mechanisms or special neurophysiological phenomena of the brain under extreme stimulation conditions, ultra-high-intensity stimulation may be used under strict experimental design and safety guarantees. In the clinical field, it is only considered for a very small number of special and difficult cases that have been comprehensively evaluated, and for which all other treatment methods have been tried and are ineffective, and only on the premise of full demonstration by a multi-disciplinary expert team and strict monitoring. At the same time, before application, relevant risks and possible consequences must be fully informed to the patient and their family members to obtain informed consent.

[0134] Under ultra-high-intensity stimulation, the stimulation intensity range (RMT%) is greater than or equal to 120%, and the stimulation frequency range (Hz) is greater than or equal to 15.

[0135] As Figures 2 - 9 It can be seen from the quantization grading results shown that by grading the experimental results through the quantization grading system proposed by the present invention, under TMS stimulation of different intensities and frequencies, the changes in neuronal activity indicators fully conform to the expectations of the quantization grading system. Specifically, when the experimental results are graded according to the first-level standard, it is expected that the actual changes in neuroelectrophysiological indicators also show small changes. For example, the membrane potential amplitude increases slightly and the frequency changes little. When the experimental results are graded according to the second-level standard, it is expected that the changes in the membrane potential amplitude and frequency also conform to the expectations of a moderate impact. When the experimental results are graded according to the third-level standard, it is expected that the membrane potential amplitude increases significantly and the frequency increases significantly, conforming to the expectations of a strong impact.

[0136] In summary, the embodiments of the present invention propose a quantitative analysis method for multi-parameter transcranial magnetic stimulation, which combines the advantages of hierarchical clustering algorithm and DBSCAN clustering algorithm. First, the DBSCAN algorithm is used to remove noise points to obtain relatively pure data, and then the hierarchical clustering algorithm is applied to these data, which can not only avoid the interference of noise on hierarchical clustering, but also give full play to the advantages of hierarchical clustering in revealing the data hierarchical structure and processing multi-variable correlations, so as to achieve complementary advantages and improve the accuracy of clustering and grading.

[0137] Those skilled in the art should understand that the embodiments of the present invention can be provided as a method, a system, or a computer program product. Therefore, the present invention can take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present invention can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk memory, CD-ROM, optical memory, etc.) that contain computer-usable program code.

[0138] The present invention is described with reference to the flowcharts and / or block diagrams of methods, apparatuses (systems), and computer program products according to the embodiments of the present invention. It should be understood that each flow and / or block in the flowchart and / or block diagram, and the combination of flows and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing devices generate means for implementing the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0139] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer-readable memory generate a manufactured article including instruction means, and the instruction means implements the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0140] These computer program instructions can also be loaded onto a computer or other programmable data processing device, so that a series of operation steps are executed on the computer or other programmable device to generate a computer-implemented process, and thus the instructions executed on the computer or other programmable device provide steps for implementing the functions specified in Figure 1 one flow or multiple flows and / or blocks Figure 1 one block or multiple blocks.

[0141] The embodiments of the present invention have been described above in conjunction with the accompanying drawings. However, the present invention is not limited to the above specific embodiments. The above specific embodiments are merely illustrative and not restrictive. Under the inspiration of the present invention, those of ordinary skill in the art can also make many forms without departing from the spirit and scope protected by the claims of the present invention. All of these are within the protection scope of the present invention.

[0142] It should be noted that in this text, relational terms such as "first" and "second" are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the term "comprising", "including" or any other variant thereof is intended to cover non-exclusive inclusion, so that a process, method, article or device comprising a series of elements not only includes those elements, but also includes other elements not expressly listed, or further includes elements inherent to such process, method, article or device. Without further limitation, an element defined by the statement "comprising a..." does not exclude the existence of additional identical elements in the process, method, article or device comprising the element.

Claims

1. A multi-parameter quantitative analysis method for transcranial magnetic stimulation, characterized in that: include: Step 1, obtaining neuron activity index data and performing standardization processing on the data to generate a first data set F1; Step 2, using the DBSCAN clustering algorithm to remove noise points from the first data set F1 to obtain a second data set F2 containing only core points and boundary points; Step 3, based on the second data set F2 obtained in step 2, calculate the distance metric and the link distance, and apply the hierarchical clustering algorithm to obtain the clustering dendrogram, combine the morphology of the clustering dendrogram, and determine the grading threshold and perform grading according to the requirements of the five-level quantitative grading system.

2. The multi-parameter transcranial magnetic stimulation quantitative analysis method according to claim 1, characterized in that: The step 2 specifically includes: Step 21, determining the preliminary range of clustering algorithm parameters, specifically including: Step 211, determining the preliminary range of the neighborhood radius ∈, specifically including: Calculate the average reachable distance ARD(∈) of data points under different ∈ values: in, N ∈ (x i ) represents the data point x i The set of data points in the ∈-neighborhood of |N ∈ (x i )| represents the data point x i The number of data points in the ∈-neighborhood, d(x i ,x j ) represents the data point x i and x j The Euclidean distance between: Where N represents the number of data points in the data set, x ik and x jk are their kth eigenvalues ​​respectively, and D is the dimension of the data; Draw a curve of ARD(∈) varying with ∈, find the elbow position of the curve, and determine the preliminary range of ∈; Step 212, determine the initial range of the MinPts value, 1*D≤MinPts≤10*D; Step 22, clustering classification and noise removal, specifically includes: Step 221, for each data point in the first data set F1, calculate the number of data points in the neighborhood ∈ determined in step 2123, if the number is greater than or equal to MinPts, mark the point as a core point, and create a new cluster, adding the point and all data points in its ∈-neighborhood to the cluster; Step 222, for each core point in the cluster, repeat the above process, adding unlabeled data points in the ∈-neighborhood of the point to the cluster until the cluster no longer expands; Step 223, if the number of data points in the ∈-neighborhood of a data point is less than MinPts, the point is temporarily marked as a noise point; Step 224, after clustering is completed, the data marked as noise points are removed from the first data set F1; Step 23, determining the final value of the neighborhood radius ∈ and the MinPts value, specifically including; Step 2321, for each ∈ value, select different MinPts values ​​for testing; Step 2322, calculate the silhouette coefficient of the corresponding clustering result Where a(i) is the average distance from data point i to other data points in its cluster, and b(i) is the minimum average distance from data point i to data points in other clusters; Step 2323, select the ∈ and MinPts combination that maximizes the silhouette coefficient s(i) as the final neighborhood radius ∈ and MinPts value. Step 24, performing clustering classification and noise removal in step 22 according to the finally determined neighborhood radius ∈ and MinPts value, and obtaining a second data set F2 containing only core points and boundary points.

3. The multi-parameter transcranial magnetic stimulation quantitative analysis method according to claim 2, characterized in that: The step 3 specifically includes: Step 31, calculating the distance metric and the link distance, specifically includes: Step 311, calculate the distance metric, for two data points x in the second data set F2 i and x j , calculate the Euclidean distance between them: Among them, x ik and x jk are their kth eigenvalues ​​respectively, and D is the dimension of the data; Step 312, calculating the average link distance, specifically includes: The average link distance method is used to calculate the distance between clusters: Among them, |A| and |B| are the number of data points in the new clusters A and B respectively, and d(x i ,x j ) is the data point x i and x j The Euclidean distance between Step 32, clustering process and dendrogram generation, specifically includes: Step 321, initialization and distance matrix calculation: Initially, each data point in the second data set F2 is regarded as a separate cluster, and the Euclidean distance between each cluster is calculated to form an M×M distance matrix, where M is the number of data points in the second data set F; Step 322, cluster merging and tree diagram updating, specifically includes: Step 3221, find the two clusters A and B with the smallest distance in the distance matrix, merge cluster A and cluster B into a new cluster, update the distance matrix, delete the rows and columns related to A and B, and calculate the distance between the new cluster and other clusters; Step 3222, repeating step 3221, when all data points are merged into one cluster, or when a preset stop condition is reached, the merging process ends, and a clustering dendrogram is gradually constructed, recording the inter-cluster distance and the corresponding clustering situation at each merging; Step 33, determining the classification threshold and classifying the levels, specifically includes: Step 331, based on the generated clustering dendrogram, determines the classification threshold, specifically including: using the method of equally dividing the height of the dendrogram to determine the classification threshold: assuming that the range of the distance between clusters in the dendrogram is [D min ,D max ], and divide it into five equal segments, with threshold points D1, D2, D3, and D4 respectively; Step 332, according to the determined classification threshold point, the data set is divided into five levels, corresponding to weak intensity stimulation, low intensity stimulation, medium intensity stimulation, high intensity stimulation and super intensity stimulation respectively.

4. The multi-parameter transcranial magnetic stimulation quantitative analysis method according to claim 3, characterized in that: The characteristic values ​​in step 3 include: relative rate of change of membrane potential, relative rate of change of action potential amplitude, relative rate of change of frequency and relative rate of decrease of discharge threshold.

Images