Multi-parameter brain stimulation effect simulation system and method based on transcranial ultrasonic simulation

Abstract

The invention relates to a multi-parameter brain stimulation effect simulation system and method based on transcranial ultrasound simulation, and the system is based on a TUSX (Transcranial Ultrasound Simulation Toolbox) toolbox, simulates the sound pressure distribution and effect of transcranial ultrasound simulation data on the brain under different parameter combinations through a high-precision three-dimensional model, and then carries out the classification and grading of simulation results through an artificial intelligence algorithm. Therefore, reference can be provided for actual clinical research.

Description

Technical Field

[0001] The present invention relates to the field of medical technology, and more specifically, to a multi-parameter brain stimulation effect simulation system and method based on transcranial ultrasound simulation. Background Art

[0002] Transcranial focused ultrasound is an emerging non-invasive brain stimulation technology that has attracted widespread attention due to its high spatial resolution and deep penetration ability. Under different stimulation parameters (such as power intensity, pulse repetition frequency, duty cycle, etc.), the effects of transcranial focused ultrasound on specific areas of the cerebral cortex (such as the M1 area) may be different. However, due to the diversity of parameter combinations, it is time-consuming and costly to study the effect of each parameter combination in clinical practice. Therefore, a method for systematic simulation in a computer is urgently needed to predict the brain response effects of different parameter combinations. Summary of the invention

[0003] In view of the above technical problems, the present invention proposes a multi-parameter brain stimulation effect simulation system and method based on transcranial focused ultrasound (tFUS) technology.

[0004] The present invention specifically relates to a multi-parameter brain stimulation effect simulation system based on transcranial ultrasound simulation, comprising:

[0005] An image acquisition module obtains high-resolution head MRI or CT images from a medical database;

[0006] A model generation module generates a three-dimensional model including skull and brain structures based on the head MRI or CT image;

[0007] Simulation calculation module, to perform multiple preset parameter combinations A i The sound field simulation under the condition of i The corresponding sound field simulation result data S i , wherein each set of parameter combinations A i Including power intensity, pulse repetition frequency, duty cycle and focal depth, sound field simulation result data S i Including each parameter combination A i The sound pressure distribution and power density obtained by simulation calculation below;

[0008] The data extraction module extracts the i The corresponding sound field simulation result data S i Extract each corresponding set of feature quantities, each set of feature quantities includes each set of parameter combination A i The sound field simulation result data S corresponding to the following i Multiple features x extracted fromij Each set of feature quantities includes at least three features of the maximum sound pressure, average sound pressure value, effective range of the sound field, sound field volume, sound field shape parameters, power density peak, temperature rise estimation, and energy absorption rate; each set of feature quantities constitutes each data point F i =(x i1 , x i2 , x ij ,…,x iM ), all data points F i Form a data set F, that is, F={F1,F2,…,F i …,F N}, where N is the number of parameter combinations, M is the number of extracted features, i is 1 to N, and j is 1 to M;

[0009] A data processing module, which performs clustering processing on the data set F, divides the clustering results into response levels, and outputs the division results;

[0010] The data processing module performs clustering processing on the data set F, specifically including:

[0011] Step 1: Density calculation. For each data point F i , calculate its local density ρ i :

[0012]

[0013] in, is the kernel function, which adopts Gaussian kernel function, and ε is the preset distance threshold;

[0014] Step 2, calculate the density weight:

[0015]

[0016] Among them, ω i is the density weight, reflecting the data point F i The density of the area;

[0017] Step 3, initialize cluster centers: select K initial cluster centers, denoted as H1, H2, …, H k ,…,H K , where k is 1 to K;

[0018] Step 4, clustering iteration: perform the clustering iteration process until the cluster center does not change or the change is less than the set threshold, including:

[0019] Step 41, cluster assignment: assign each data point F i Assign to the cluster C that is closest to it k :

[0020]

[0021] Step 42, cluster center update: update the cluster center according to the weighted average of the data points in the cluster:

[0022]

[0023] Step 43, convergence judgment: if the change of the cluster center is less than the set threshold δ, the algorithm converges and stops iterating; otherwise, continue to iterate from step 41 to step 43;

[0024] Step 5: Output the clustering results. When the clustering stops iterating, output the cluster C to which each data point belongs. k And the cluster C k Center H k .

[0025] Furthermore, each group of characteristic quantities includes maximum sound pressure, effective range of sound field and total power density.

[0026] Furthermore, it also includes a data standardization module, which performs data standardization processing on each feature extracted by the data extraction module, including:

[0027] For each extracted feature x ij , calculate its mean μ ij and standard deviation σ ij ;

[0028] Normalization is performed using the following formula:

[0029] Furthermore, the step of selecting K initial cluster centers specifically includes:

[0030] Step 31, randomly select the first cluster center: randomly select a data point from the data set as the first cluster center H1;

[0031] Step 32, calculate the distance: for each data point F in the data set i , calculate its distance D(F i ,H k ):

[0032]

[0033] Among them, H k is the center of the kth cluster, D(F i ,H k ) is the weighted distance;

[0034] Step 33, select the next cluster center: according to the probability distribution The probability of selecting the next cluster center is proportional to the square of the distance;

[0035] Step 34, repeat steps 42 and 43 until K clusters C1, C2, ..., C k ,…,C K and their corresponding cluster centers H1, H2,…, H k ,…,H K .

[0036] Furthermore, the clustering results are divided into response levels, specifically including:

[0037] Step 1, calculate the data point F i The weight coefficients of each set of feature quantities in include:

[0038]

[0039] Among them, w i is the weight coefficient of the i-th group of features, is the variance of the feature quantity of the i-th group;

[0040] Step 2: Calculate F for each data point i In the cluster C k The composite response score SC i , specifically including:

[0041]

[0042] Among them, SC i is the data point F i In the cluster C k The overall response score in SC i,j is the data point F i In the cluster C k The standardized score of the jth feature in ;

[0043] Step 3: Based on the comprehensive response score SC i Get each data point F i The corresponding response level.

[0044] Furthermore, we output the cluster C to which each data point belongs. k and cluster center H k Previously, it also included dynamically adjusting the number of clusters K, including:

[0045] Step 51', cluster merging: If the distance between two cluster centers is ‖H k′ -H k″‖ is smaller than the set merging threshold d merge , then merge these two clusters;

[0046] Step 52', cluster splitting: If the local density variance of the data points within a cluster is Exceeding the set split threshold The cluster is split and a new cluster center is added.

[0047] On the other hand, the present invention also proposes a multi-parameter brain stimulation effect simulation method based on transcranial ultrasound simulation, the method comprising:

[0048] An image acquisition step, obtaining high-resolution head MRI or CT images from a medical database;

[0049] A model generation step, generating a three-dimensional model including a skull and brain structure based on the head MRI or CT image;

[0050] Simulate the calculation steps and perform multiple preset parameter combinations A i The sound field simulation under the condition of i The corresponding sound field simulation result data S i , wherein each set of parameter combinations A i Including power intensity, pulse repetition frequency, duty cycle and focal depth, sound field simulation result data S i Including each parameter combination A i The sound pressure distribution and power density obtained by simulation calculation below;

[0051] Data extraction step, from each parameter combination A i The corresponding sound field simulation result data S i Extract the corresponding set of features from the above, which constitute each set of parameter combination A. i The sound field simulation result data S corresponding to the following i The data point value extracted from: F i =(x i1 , x i2 , x ij ,…,x iM ), all data points F i Form a data set, forming a data set F, that is, F = {F1, F2, ..., F i …,F N}, where M is the number of extracted features, N is the number of parameter combinations, i is 1 to N, and j is 1 to M;

[0052] A data processing step, after clustering the data set F, divides the clustering results into response levels and outputs the division results; wherein the clustering of the data set F in the data processing step specifically includes:

[0053] Step 1: Density calculation. For each data point F i , calculate its local density ρ i :

[0054]

[0055] in, is the kernel function, which adopts Gaussian kernel function, and ε is the preset distance threshold;

[0056] Step 2, calculate the density weight:

[0057]

[0058] Among them, ω i is the density weight, reflecting the data point F i The density of the area;

[0059] Step 3, initialize cluster centers: select K initial cluster centers, denoted as H1, H2, …, H k ,…,H K , where k is 1 to K;

[0060] Step 4, clustering iteration: perform the clustering iteration process until the cluster center does not change or the change is less than the set threshold, including:

[0061] Step 41, cluster assignment: assign each data point F i Assign to the cluster C that is closest to it k :

[0062]

[0063] Step 42, cluster center update: update the cluster center according to the weighted average of the data points in the cluster:

[0064]

[0065] Step 43, convergence judgment: if the change of the cluster center is less than the set threshold δ, the algorithm converges and stops iterating; otherwise, continue to iterate from step 41 to step 43;

[0066] Step 5: Output the clustering results. When the clustering stops iterating, output the cluster C to which each data point belongs. k And the cluster C k Center H k.

[0067] The present invention proposes a multi-parameter brain stimulation effect simulation system and method based on transcranial ultrasound simulation, which utilizes the TUSX toolbox and the k-Wave acoustic simulation package to realize individualized sound field simulation based on high-resolution medical images for the first time, and can provide more accurate simulation capabilities for transcranial ultrasound stimulation.

[0068] In addition, the present invention also proposes an adaptive clustering algorithm based on weighted density and spatial features (ACA-WDSF). The algorithm combines the local density, spatial position characteristics and weighting mechanism of data points, and has the following characteristics: 1. Combining density and spatial features: considering the local density and spatial distribution characteristics of data points, and enhancing the recognition ability of nonlinear structures and complex shapes. 2. Multi-dimensional feature fusion: integrating multi-dimensional features such as sound pressure, power density, sound field volume, spatial position and biological effects to comprehensively evaluate brain response. 3. Dynamic cluster number adjustment: automatically determining the optimal number of clusters through an adaptive mechanism, avoiding the complexity and uncertainty of manually setting the number of clusters. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] The drawings described herein are used to provide a further understanding of the present disclosure and constitute a part of the present disclosure. The illustrative embodiments of the present disclosure and their descriptions are used to explain the present disclosure and do not constitute an improper limitation on the present disclosure. In the drawings:

[0070] Figure 1 A schematic structural diagram of an embodiment of a multi-parameter brain stimulation effect simulation system of the present invention is shown;

[0071] Figure 2 A schematic flow chart of an embodiment of the multi-parameter brain stimulation effect simulation method of the present invention is shown. DETAILED DESCRIPTION

[0072] In order to make the purpose, 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 in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments.

[0073] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0074] like Figure 1 As shown, a schematic diagram of the structure of a multi-parameter brain stimulation effect simulation system based on transcranial focused ultrasound (tFUS) technology provided by a specific embodiment of the present invention is shown. The system may include: an image acquisition module, a model generation module, a simulation calculation module, a data extraction module, and a data processing module.

[0075] The following is a detailed introduction to the implementation and function of each module:

[0076] An image acquisition module is used to acquire high-resolution head MRI or CT images from a medical database. In the present invention, the head MRI or CT image may also be head MRI or CT image data directly acquired by an MRI scanning device or a CT scanning device for training a target model, and no specific limitation is made here.

[0077] A model generation module generates a three-dimensional model including skull and brain structures based on the acquired head MRI or CT image. In one embodiment of the present invention, the head MRI or CT image can be imported into the TUSX toolbox, and a three-dimensional model including skull and brain structures is generated through a preprocessing module of the TUSX toolbox.

[0078] In a specific embodiment, in the TUSX toolbox, the position and orientation of the ultrasonic transducer can be set using a selected combination of parameters, the coordinates of the M1 region where the center of the transducer is aligned can be determined, and the focus depth and the resolution of the sound field simulation grid can be set. In addition, the TUSX toolbox can be used in conjunction with the k-Wave acoustic simulation package to perform numerical simulations of three-dimensional sound fields. Each parameter combination corresponds to a simulation, recording the propagation path of the ultrasonic wave in the M1 region, the sound pressure distribution, and its power density map.

[0079] In this embodiment, a plurality of preset parameter combinations A are performed. i The sound field simulation under the condition of i They include: power intensity, unit W / cm 2 , pulse repetition frequency, unit Hz, duty cycle, unit %, focal depth, unit mm, each parameter combination A i Generate a sound field simulation result S i , that is, i is 1~N.

[0080] In another embodiment, it further includes using the interpolation and smoothing functions provided by TUSX to refine the three-dimensional model to ensure the accuracy of the simulation; for MRI data, a uniform skull medium model is used; for CT data, a uniform or heterogeneous skull medium model can be selected. Through the interpolation and smoothing processing, the accuracy of the simulated data can be effectively improved, and different models are used for MRI data and CT data, which can better ensure the accuracy of different types of data.

[0081] In addition, based on the spatial and physical characteristics of transcranial focused ultrasound (tFUS) data, relevant features reflecting brain responses can be extracted to analyze transcranial ultrasound stimulation. These features can generally be divided into the following categories:

[0082] 1. Spatial distribution characteristics:

[0083] Sound field volume (equivalent sound pressure area) V P :The sound pressure value exceeds the threshold P th The volume of space is:

[0084]

[0085] Among them, δ(·) is the indicator function, when P i ≥P th When δ=1, otherwise δ=0; V voxel The volume of each voxel (m 3 ).

[0086] Sound field shape parameters: such as length-to-width ratio and volume-to-surface area ratio, which reflect the geometric characteristics of the sound field.

[0087] 2. Strength characteristics:

[0088] Maximum sound pressure value P max : The highest sound pressure value in the entire sound field.

[0089] Average sound pressure The average value of the sound pressure value is:

[0090]

[0091] Where N is the total number of sound field grid points, P i is the sound pressure value at the i-th grid point.

[0092] Peak power density W max : The highest power density value in the entire sound field.

[0093] 3. Spatial location characteristics:

[0094] Sound field center coordinates (x c ,y c ,z c ): The center of mass position of the sound field, the formula is:

[0095]

[0096] Among them, x i ,y i ,z i is the spatial coordinate of the i-th grid point, W i is the corresponding power density.

[0097] Focus offset: simulates the distance between the focus and the target area (such as M1 area) to evaluate the focusing accuracy.

[0098] 4. Biological effect characteristics:

[0099] Energy absorption rate: The energy absorption rate is calculated based on the sound field data and tissue characteristics, reflecting the absorption of sound energy in the tissue.

[0100] Temperature rise estimation: Using the biothermal model, the possible temperature rise is estimated based on the power density and tissue thermal conductivity characteristics. The formula is:

[0101]

[0102] Where Q is the absorbed heat (J), ρ is the tissue density (kg / m 3 ), c p is the specific heat capacity (J / (kg·K)), V is the heated volume (m 3 ).

[0103] These features reflect the effects of tFUS on the brain under different parameters, and provide rich information for subsequent cluster analysis and response classification. In the present invention, in order to facilitate subsequent simulation and prediction analysis, some of the features are extracted through a data extraction module.

[0104] The data extraction module extracts the i The corresponding sound field simulation result data S i Extract the corresponding set of features from the above, which constitute each set of parameter combination A. i The sound field simulation result data S corresponding to the following i The data point value extracted from: F i =(x i1 , x i2 , x ij ,…,x iM), all data points F i Form a data set, forming a data set F, that is, F = {F1, F2, ..., F i …,F N}, where N is the number of parameter combinations, M is the number of extracted features, i is 1~N, and j is 1~M.

[0105] In one embodiment, each set of feature quantities includes at least three features of maximum sound pressure, average sound pressure value, effective range of sound field, sound field volume, sound field shape parameter, power density peak, temperature rise estimation, and energy absorption rate. In this embodiment, the present invention can comprehensively evaluate brain response by integrating multi-dimensional features such as sound pressure, power density, sound field volume, spatial position, and biological effects.

[0106] In another preferred embodiment, the extracted feature quantities include maximum sound pressure, effective range of sound field and total power density, wherein the maximum sound pressure is the highest sound pressure value in the sound field, which is used to describe the intensity characteristics of the sound field, the effective area of ​​the sound field is the area where the sound pressure exceeds a certain threshold, reflecting the coverage of ultrasound in the target area, and the total power density represents the sum of all energy in the sound field, reflecting the energy delivery of ultrasound to the target area. In this embodiment, by selecting only the three feature quantities of maximum sound pressure, effective range of sound field and total power density, not only can the overall amount of calculation be reduced, but also the adaptability and robustness of transcranial focused ultrasound data analysis can be guaranteed to the greatest extent.

[0107] In another embodiment, the brain stimulation effect simulation system based on transcranial focused ultrasound technology in the present invention may further include a data standardization module, which performs a data normalization on each feature point F extracted by the data extraction module. i Perform data standardization, including: for each extracted feature x ij , calculate its mean μ ij and standard deviation σ ij ; Normalized by the following formula: The standardized data has a mean of 0 and a standard deviation of 1, which ensures that different features have the same weight in the clustering process.

[0108] The data processing module is used to perform clustering processing on the data set F, divide the clustering results into response levels, and finally output the division results.

[0109] Since the distribution density of the extracted feature quantity data in the feature space varies greatly, there are high-density and low-density areas, and the density differences in different areas are significant. In order to facilitate the automatic classification of the simulation results, the present invention selects an unsupervised learning algorithm to perform clustering preprocessing on the data. Commonly used unsupervised learning algorithms include: K-means clustering (K-Means), which is usually suitable for situations where the data points are relatively dense and the cluster shape is spherical, and hierarchical clustering (HierarchicalClustering), which is usually suitable for situations with complex data structures. In the present invention, since the data points are relatively uniform and we hope to obtain a predetermined number of cluster categories, the K-means clustering method is selected.

[0110] In order to overcome the shortcomings of traditional clustering algorithms in processing transcranial focused ultrasound (tFUS) simulation data, the present invention designs an adaptive clustering algorithm based on weighted density and spatial features (ACA-WDSF). The algorithm combines the local density, spatial location characteristics and weighting mechanism of data points. In a preferred embodiment of the present invention, clustering processing is performed on the data set F, specifically including:

[0111] Step 1: Density calculation. For each data point F i , calculate its local density ρ i :

[0112]

[0113] in, is the kernel function, which adopts Gaussian kernel function, and ε is the preset distance threshold;

[0114] In one embodiment, a Gaussian kernel function may be selected as the kernel function for density estimation, and the formula is:

[0115]

[0116] Here, σ is the kernel bandwidth, which determines the sensitivity of the density calculation. A smaller σ value will make the density estimation more localized, which is suitable for fine-grained analysis in high-density areas; a larger σ value is suitable for global density estimation in low-density areas.

[0117] In a specific embodiment, how to calculate the preset distance threshold ε specifically includes: the data distribution of the three characteristic quantities of maximum sound pressure, effective range of sound field and total power density obtained by the data extraction module, and respectively calculating the mean μ and standard deviation σ of each characteristic quantity, and setting the distance threshold ε to a specific multiple of the characteristic quantity standard deviation σ. In a preferred embodiment, ε can be taken between 0.5σ and 2σ. Generally speaking, if the data distribution is relatively concentrated and the standard deviation is small, then ε can be taken within a smaller multiple of the standard deviation range; conversely, if the data distribution is more dispersed and the standard deviation is larger, then ε can be taken within a larger multiple of the standard deviation range.

[0118] Step 2, calculate the density weight coefficient:

[0119]

[0120] Among them, ω i is the density weight coefficient, reflecting the data point F i The higher the density of the area, the greater the weight, the greater the proportion of the distance measurement, which makes the data points in the high-density area have a stronger influence on the cluster center.

[0121] In another embodiment, a variance-based method can also be used to calculate the density weight ω i , specifically:

[0122]

[0123] in, is the variance of the i-th feature. The features with larger variance have smaller weights to prevent them from having too much influence on the clustering results; the features with smaller variance have larger weights to enhance their influence on the clustering results.

[0124] Step 3, initialize cluster centers: select K initial cluster centers, denoted as H1, H2, …, H k ,…,H K , where k is 1~K.

[0125] In one embodiment, the selecting K initial cluster centers specifically includes:

[0126] Step 31, randomly select the first cluster center: randomly select a data point from the data set as the first cluster center H1;

[0127] Step 32, calculate the distance: for each data point F in the data set i , calculate its distance D(F i ,H k ):

[0128]

[0129] Among them, H k is the center of the kth cluster, D(F i ,H k ) is the weighted distance;

[0130] Step 33, select the next cluster center: according to the probability distribution The next cluster center is selected with probability proportional to the square of the distance.

[0131] Specifically: First, for each data point F in the data set i (For example, the characteristic vector F composed of the maximum sound pressure, the effective range of the sound field and the total power density i ), calculate its weighted distance D(F i ,H k ). For example, if there are two selected cluster centers H1 and H2, for the data point F i , its weighted distance D(F i ,H k )=min{ω i ‖F i -H1‖,ω i ‖F i -H2‖}, where ‖·‖ represents the Euclidean norm of the distance calculation.

[0132] In the present invention, based on the traditional Euclidean distance, the density weight coefficient w is introduced i , thus forming a density-weighted distance metric. This method calculates the local density of each data point and gives higher weights to data points in high-density areas, so that the clustering process pays more attention to these areas and improves the accuracy of clustering.

[0133] Secondly, according to the formula Calculate F for each data point i The probability of being selected as the next cluster center includes: first calculating the denominator ∑ Fi∈F D(F i ,H k ) 2 , that is, calculate the cluster C k All data points in F i To the nearest selected cluster center H k The sum of the squares of the distances; then for each data point F i , calculate its distance to the nearest selected cluster center H k The square of the distance D(F i ,H k ) 2 , and divide it by the denominator to get the selection probability P(F i);

[0134] Next, generate a random number r uniformly distributed in the interval [0,1], and then select k Starting from the first data point in the , the selection probability P(F i ) until the cumulative sum is greater than or equal to the random number r, at which point the corresponding data point is selected as the next cluster center.

[0135] In a specific embodiment, if there are three data points F1, F2, and F3 in the data set F, their selection probabilities are calculated to be P(F1)=0.2, P(F2)=0.3, and P(F3)=0.5 respectively; at this time, a random number r=0.4 is randomly generated, and P(F1)=0.2 is first accumulated. At this time, the accumulated sum is less than r, and P(F2)=0.3 is continued to be accumulated. At this time, the accumulated sum is 0.2+0.3=0.5, which is greater than or equal to r, and F2 is selected as the next cluster center.

[0136] Step 34, repeat steps 42 and 43 until K clusters C1, C2, ..., C k ,…,C K and their corresponding cluster centers H1, H2,…, H k ,…,H K .

[0137] Step 4, clustering iteration: perform the clustering iteration process until the cluster center does not change or the change is less than the set threshold, including:

[0138] Step 41, cluster assignment: assign each data point F i Assign to the cluster C that is closest to it k :

[0139]

[0140] Step 42, cluster center update: update the cluster center according to the weighted average of the data points in the cluster:

[0141]

[0142] It should be noted that the present invention simultaneously considers the density weight of data points through the updating method, which can make the cluster center more inclined to the data points in the high-density area, thereby enhancing the stability and accuracy of clustering.

[0143] Step 43, convergence judgment: if the change of the cluster center is less than the set threshold δ, the algorithm converges and stops iterating; otherwise, continue to iterate from step 41 to step 43;

[0144] In the present invention, the goal of K-means clustering is to minimize the following objective function:

[0145]

[0146] in:

[0147] r ik : represents data point F i Whether it is assigned to cluster C k , 1 if yes, 0 otherwise.

[0148] ||F i -H k ||: represents data point F i With cluster center H k The Euclidean distance between .

[0149] It can be understood that those skilled in the art can set the threshold δ through various methods, and the present invention does not make any specific limitation here.

[0150] Step 5: Output the clustering results. When the clustering stops iterating, output the cluster C to which each data point belongs. k And the cluster C k Center H k .

[0151] In the whole clustering process, C k represents the i-th data point F i The cluster to which it is assigned. When the clustering stops iterating, output each data point F i The final cluster C k And the cluster centers H of all clusters k (k=1,…,K, K is the number of clusters). Among them, the cluster center H k To some extent, it represents cluster C k The typical characteristics of the data points in the cluster. In the early processing such as data standardization, the location of the cluster center will affect the distance calculation from the data point to the cluster center (such as when calculating the weighted distance d(F i ,H k ), and the distribution of data points (related to the cluster center) will affect the statistical characteristics of the feature quantity (such as variance, etc.), and then affect the density weight coefficient ω i and standardized scoring SC i,j In addition, the comprehensive response score SC i It is calculated based on the data processing results related to the cluster center, which reflects the cluster C k The overall response level, and the cluster center plays a fundamental role in this process, affecting the quantitative evaluation of the cluster characteristics by the entire scoring system.

[0152] In another embodiment, when outputting the cluster C to which each data point belongs, k and cluster center H k Previously, it also included dynamically adjusting the number of clusters K, including:

[0153] Step 51', cluster merging: If the distance between two cluster centers is ‖H k′ -H k″ ‖ is smaller than the set merging threshold d merge , then merge the two clusters.

[0154] In a specific embodiment, how to calculate the merging threshold d merge Specifically, if there are N = 1000 data points consisting of the maximum sound pressure, effective range of the sound field and total power density obtained by the data extraction module in the brain stimulation effect simulation data, then the Euclidean distance between the 1000 data points is calculated to obtain a distance matrix. If the mean value μ of all distances in the distance matrix is ​​calculated d =5.0, standard deviation σ d =2.0, then first set d merge =μ d =5.0.

[0155] At this time, the first clustering algorithm calculation is performed to obtain the clustering result C1, and the number of clusters K1 and the distribution of data points in the clusters (such as the mean and variance of the data points in each cluster in terms of maximum sound pressure, effective range of the sound field, and total power density characteristics) are recorded; then change d merge The value of merge =3.0 to perform the second clustering algorithm calculation, obtain the clustering result C2, record the corresponding cluster number K2 and the distribution statistics of the data points in the cluster; similarly, take d merge =7.0, perform the third clustering algorithm calculation, obtain the clustering result C3 and the number of clusters K3, and record the relevant information; compare the changes of K1, K2, and K3. If K2 changes greatly compared with K1 and K3, it means that d merge =3.0 is not appropriate, resulting in unstable clustering results. If K1 and K3 are relatively close, further compare the changes in the distribution statistics of the data points in the clusters C1 and C3. If the change is small, it means that the clustering results are unstable. merge =5.0 and d merge =7.0, the clustering results are relatively stable, and we can continue to conduct more detailed value experiments within this stable range, such as taking d merge =6.0, etc., until a d is found that makes the clustering result most stable in terms of the number of clusters and the distribution of data points. merge value.

[0156] In addition, in another embodiment, how to determine the merging threshold d merge , for cluster C k Each data point F i (Here F i Corresponding to the characteristic vector F composed of the maximum sound pressure, effective range of the sound field and total power density i ), calculate all data points F i The Euclidean distance between each pair is used to obtain the distance matrix. If the mean μ of all distances in the distance matrix is ​​calculated, d , standard deviation σ d , then set the merge threshold d merge The range is μ d -σ d <d merge <μ d +σ d .

[0157] Step 52', cluster splitting: If the local density variance of the data points within a cluster is Exceeding the set split threshold The cluster is split and a new cluster center is added.

[0158] In a specific embodiment, how to determine the splitting threshold For cluster C k Each data point F i (Here F i Corresponding to the characteristic vector F composed of the maximum sound pressure, effective range of the sound field and total power density i ), calculate all data points F i The Euclidean distance between each pair is used to obtain the distance matrix. If the mean μ of all distances in the distance matrix is ​​calculated, d , standard deviation σ d , then set the split threshold Range is

[0159] In one embodiment, the present invention further proposes a multi-dimensional scoring method for weighting and summarizing the scores of each key feature to obtain a comprehensive response score for each data point, thereby achieving the division of response levels of the simulation results, which specifically includes:

[0160] Step 1, calculate the data point F i The weight coefficients of each set of feature quantities in include:

[0161]

[0162] Among them, w i is the weight coefficient of the i-th group of features, is the variance of the feature quantity of the i-th group.

[0163] In addition, in another embodiment, the weight coefficients of the characteristic quantities may also be determined by pre-customization. For example, in one embodiment, the weight coefficient of the maximum sound pressure value is 0.25, the weight coefficient of the average sound pressure value is 0.2, the weight coefficient of the sound field volume is 0.15, the weight coefficient of the sound field shape parameter is 0.1, the weight coefficient of the power density peak is 0.2, and the weight coefficient of the temperature rise estimation is 0.1. For example, in another embodiment, the weight coefficient of the maximum sound pressure value is 0.3, the weight coefficient of the average sound pressure value is 0.15, the weight coefficient of the sound field volume is 0.1, the weight coefficient of the sound field shape parameter is 0.05, the weight coefficient of the power density peak is 0.25, and the weight coefficient of the temperature rise estimation is 0.15.

[0164] Step 2: Calculate F for each data point i In the cluster C k The composite response score SC i , specifically including:

[0165]

[0166] Among them, SC i is the data point F i In the cluster C k The overall response score in SC i,j is the data point F i In the cluster C k The normalized score of the j-th feature in .

[0167] In one embodiment, the standardized score SC can be calculated by the following method: i,j :First calculate the cluster C k All data points in F i The mean of the jth feature (where n k is cluster C k The number of data points in the inner part), and then according to the data point F i The jth feature x i,j With the above mean The standardized score SC is determined by the proportional relationship i,j For example, if the data point F i The jth feature x i,j Greater than the feature mean 75% of SC i,j A higher score, such as 4 or 5, indicates that the cluster has a higher degree of response on this feature; if the data point F i The jth feature x i,jLess than the feature mean 25%, then SC i,j Assigned a lower score, such as 1 or 2, indicating that the cluster has a lower response degree on this feature; when the data point F i The jth feature x i,j Less than the feature mean 75% and greater than the feature mean When the SC i,j A medium score, such as 3, indicates that the cluster has a medium degree of response on this feature.

[0168] Step 3: Obtain a corresponding response level according to the comprehensive response score.

[0169] In one embodiment, based on the biological effect threshold and safety standards, the comprehensive score is 0-20 for level one, 21-40 for level two, and so on, the response level can be divided into the following five levels:

[0170] Level 1 (safety low response): All key characteristic values ​​are within the safe low range and no significant biological effects are expected. The overall score is between 0 and 20 points.

[0171] Level 2 (mild response): Some characteristic values ​​are slightly increased, which may produce slight biological effects, but it is still within the safe range. The overall score is 20-40 points.

[0172] Level 3 (moderate response): The key characteristic values ​​are at an intermediate level, with obvious biological effects and high overall safety, with a comprehensive score of 40-60 points.

[0173] Level 4 (significant response): The characteristic value is high and the biological effect is significant, but attention should be paid to potential safety risks. The overall score is between 60 and 80 points.

[0174] Level 5 (dangerous response): The key characteristic values ​​far exceed the safety standards, and there is a risk of excessive stimulation. It needs to be used with caution or prohibited. The overall score is 80-100 points.

[0175] In one embodiment, each cluster is assigned to a corresponding response level according to its comprehensive response score. For example, cluster C1 has a comprehensive score of 10 points, which is assigned to level 1; cluster C3 has a comprehensive score of 35 points, which is assigned to level 3.

[0176] The present invention verifies the brain stimulation effect simulation system through a water tank experiment. The human ex vivo skull is placed in a water tank. The experimental equipment includes an ATA-8202 power amplifier, an EDU33212A signal source, etc. The ultrasonic parameters are set according to the simulation results. In the spatial peak time average sound intensity test, the parameters are kept stable during the stimulation depth test. According to the peak-to-peak value of the specific excitation signal voltage and the power amplifier gain, the average sound power before and after the skull attenuation is measured with a radiation force balance, and the response level of the water tank experiment result is calculated and compared with the predicted result of the brain stimulation effect simulation system of the present invention.

[0177] Table 1

[0178]

[0179] As shown in Table 1, it is a comparison between the response level corresponding to the water tank experiment results and the response level corresponding to the simulation results proposed by the present invention by integrating multi-dimensional features such as sound pressure, power density, sound field volume, spatial position and biological effect. It can be seen that the prediction algorithm proposed by the present invention is completely consistent with the water tank experiment results, indicating that the brain stimulation effect prediction system proposed by the present invention can achieve accurate quantitative and graded prediction of brain response.

[0180] Table 2

[0181]

[0182] As shown in Table 2, it is a comparison between the response level corresponding to the water tank experiment results and the response level corresponding to the simulation results of only selecting the three characteristic quantities of maximum sound pressure, effective range of sound field and total power density proposed in the present invention. It can be seen that the simulation results of only selecting the three characteristic quantities of maximum sound pressure, effective range of sound field and total power density are also completely consistent with the water tank experiment results, indicating that the brain stimulation effect prediction system proposed in the present invention can still achieve accurate quantitative graded prediction of brain response when selecting fewer characteristic quantities, thereby improving the computational efficiency.

[0183] like Figure 2 As shown, in another embodiment of the present disclosure, a multi-parameter combined brain stimulation effect prediction method based on transcranial ultrasound simulation is provided, the method comprising: an image acquisition step, acquiring a high-resolution head MRI or CT image from a medical database;

[0184] A model generation step, generating a three-dimensional model including a skull and brain structure based on the head MRI or CT image;

[0185] Simulate the calculation steps and perform multiple preset parameter combinations A i Sound field simulation, calculate each parameter combination A iThe corresponding sound field simulation result data S i , wherein each set of parameter combinations A i Including power intensity, pulse repetition frequency, duty cycle and focal depth, sound field simulation result data S i Including each parameter combination A i The sound pressure distribution and power density obtained by simulation calculation below;

[0186] Data extraction step, from each parameter combination A i The corresponding sound field simulation result data S i Extract each corresponding set of feature quantities, each set of feature quantities includes each set of parameter combination A i The sound field simulation result data S corresponding to the following i Multiple features x extracted from ij Each set of feature quantities includes at least three features of the maximum sound pressure, average sound pressure value, effective range of the sound field, sound field volume, sound field shape parameters, power density peak, temperature rise estimation, and energy absorption rate; each set of feature quantities constitutes each data point F i =(x i1 , x i2 , x ij ,…,x iM ), all data points F i Form a data set F, that is, F={F1,F2,…,F i …,F N}, where N is the number of parameter combinations, M is the number of extracted features, i is 1 to N, and j is 1 to M;

[0187] A data processing step, after clustering the data set F, divides the clustering results into response levels and outputs the division results; wherein the clustering of the data set F in the data processing step specifically includes:

[0188] Step 1: Density calculation. For each data point F i , calculate its local density ρ i :

[0189]

[0190] in, is the kernel function, which adopts Gaussian kernel function, and ε is the preset distance threshold;

[0191] Step 2, calculate the density weight:

[0192]

[0193] Among them, ω i is the density weight, reflecting the data point Fi The density of the area;

[0194] Step 3, initialize cluster centers: select K initial cluster centers, denoted as H1, H2, …, H k ,…,H K , where k is 1 to K;

[0195] Step 4, clustering iteration: perform the clustering iteration process until the cluster center does not change or the change is less than the set threshold, including:

[0196] Step 41, cluster assignment: assign each data point F i Assign to the cluster C that is closest to it k :

[0197]

[0198] Step 42, cluster center update: update the cluster center according to the weighted average of the data points in the cluster:

[0199]

[0200] Step 43, convergence judgment: if the change of the cluster center is less than the set threshold δ, the algorithm converges and stops iterating; otherwise, continue to iterate from step 41 to step 43;

[0201] Step 5: Output the clustering results. When the clustering stops iterating, output the cluster C to which each data point belongs. k And the cluster C k Center H k .

[0202] In summary, the present invention proposes a brain stimulation effect simulation system and method based on transcranial focused ultrasound technology, and uses the TUSX toolbox and k-Wave acoustic simulation package to realize individualized sound field simulation based on high-resolution medical images for the first time, providing more accurate prediction capabilities for transcranial ultrasound stimulation. In addition, the present invention also proposes an adaptive clustering algorithm based on weighted density and spatial features, which realizes the automated analysis of complex sound field data and the classification of neural responses, greatly improving the research efficiency and the repeatability of the results.

[0203] Those skilled in the art will appreciate that embodiments of the present invention may be provided as methods, systems, or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0204] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes 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 a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0205] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0206] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0207] The embodiments of the present invention are described above in conjunction with the accompanying drawings, but the present invention is not limited to the above-mentioned specific implementation methods. The above-mentioned specific implementation methods are merely illustrative and not restrictive. Under the enlightenment of the present invention, ordinary technicians in this field can also make many forms without departing from the scope of protection of the purpose of the present invention and the claims, which all fall within the protection of the present invention.

[0208] It should be noted that, in this article, 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 such actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device. In the absence of further restrictions, the elements defined by the sentence "comprise a ..." do not exclude the existence of other identical elements in the process, method, article or device including the elements.

Claims

1. A multi-parameter brain stimulation effect simulation system based on transcranial ultrasound simulation, characterized in that: include: An image acquisition module obtains high-resolution head MRI or CT images from a medical database; A model generation module generates a three-dimensional model including skull and brain structures based on the head MRI or CT image; Simulation calculation module, to perform multiple preset parameter combinations A i The sound field simulation under the condition of i The corresponding sound field simulation result data S i , wherein each set of parameter combinations A i Including power intensity, pulse repetition frequency, duty cycle and focal depth, sound field simulation result data S i Including each parameter combination A i The sound pressure distribution and power density obtained by simulation calculation below; The data extraction module extracts the i The corresponding sound field simulation result data S i Extract each corresponding set of feature quantities, each set of feature quantities includes each set of parameter combination A i The sound field simulation result data S corresponding to the following i Multiple features x extracted from ij Each set of feature quantities includes at least three features of the maximum sound pressure, average sound pressure value, effective range of the sound field, sound field volume, sound field shape parameters, power density peak, temperature rise estimation, and energy absorption rate; each set of feature quantities constitutes each data point F i =(x i1 , x i2 , x ij ,…,x iM ), all data points F i Form a data set F, that is, F={F1,F2,…,F i …,F N }, where N is the number of parameter combinations, M is the number of extracted features, i is 1 to N, and j is 1 to M; A data processing module, which performs clustering processing on the data set F, divides the clustering results into response levels, and outputs the division results; The clustering process of the data set F in the data processing module specifically includes: Step 1: Density calculation. For each data point F i , calculate its local density ρ i : in, is the kernel function, which adopts Gaussian kernel function, and ε is the preset distance threshold; Step 2, calculate the density weight: Among them, ω i is the density weight, reflecting the data point F i The density of the area; Step 3, initialize cluster centers: select K initial cluster centers, denoted as H1, H2, …, H k ,…,H K , where k is 1 to K; Step 4, clustering iteration: perform the clustering iteration process until the cluster center does not change or the change is less than the set threshold, including: Step 41, cluster assignment: assign each data point F i Assigned to the cluster C that is closest to it k : Step 42, cluster center update: update the cluster center according to the weighted average of the data points in the cluster: Step 43, convergence judgment: if the change of the cluster center is less than the set threshold δ, the algorithm converges and stops iterating; otherwise, continue to iterate from step 41 to step 43; Step 5: Output the clustering results. When the clustering stops iterating, output the cluster C to which each data point belongs. k And the cluster C k Center H k .

2. The multi-parameter brain stimulation effect simulation system according to claim 1, characterized in that: Each group of characteristic quantities includes maximum sound pressure, effective range of sound field and total power density.

3. The multi-parameter brain stimulation effect simulation system according to claim 2, characterized in that: It also includes a data standardization module, which performs data standardization processing on each feature extracted by the data extraction module, including: For each extracted feature x ij , calculate its mean μ ij and standard deviation σ ij ; Normalization is performed using the following formula:

4. The multi-parameter brain stimulation effect simulation system according to claim 3, characterized in that: The selecting of K initial cluster centers specifically includes: Step 31, randomly select the first cluster center: randomly select a data point from the data set as the first cluster center H1; Step 32, calculate the distance: for each data point F in the data set i , calculate its distance D(F i ,H k ): Among them, H k is the center of the kth cluster, D(F i ,H k ) is the weighted distance; Step 33, select the next cluster center: according to the probability distribution The probability of selecting the next cluster center is proportional to the square of the distance; Step 34, repeat steps 42 and 43 until K clusters C1, C2, ..., C k ,…,C K and their corresponding cluster centers H1, H2,…, H k ,…,H K .

5. The multi-parameter brain stimulation effect simulation system according to claim 4, characterized in that: The classification of the clustering results into response levels specifically includes: Step 1, calculate the data point F i The weight coefficients of each set of feature quantities in include: Among them, w i is the weight coefficient of the i-th group of features, is the variance of the feature quantity of the i-th group; Step 2: Calculate F for each data point i In the cluster C k The composite response score SC i , specifically including: Among them, SC i is the data point F i In the cluster C k The overall response score in SC i,j is the data point F i In the cluster C k The standardized score of the jth feature in ; Step 3: Based on the comprehensive response score SC i Get each data point F i The corresponding response level.

6. The multi-parameter brain stimulation effect simulation system according to claim 5, characterized in that: Output the cluster C to which each data point belongs k and cluster center H k Previously, it also included dynamically adjusting the number of clusters K, including: Step 51', cluster merging: If the distance between two cluster centers is ‖H k′ -H k″ ‖ is smaller than the set merging threshold d merge , then merge these two clusters; Step 52', cluster splitting: If the local density variance of the data points within a cluster is Exceeding the set split threshold The cluster is split and a new cluster center is added.

7. A multi-parameter brain stimulation effect simulation method based on transcranial ultrasound simulation, the method comprising: An image acquisition step, obtaining high-resolution head MRI or CT images from a medical database; A model generation step, generating a three-dimensional model including a skull and brain structure based on the head MRI or CT image; Simulate the calculation steps and perform multiple preset parameter combinations A i The sound field simulation under the condition of i The corresponding sound field simulation result data S i , wherein each set of parameter combinations A i Including power intensity, pulse repetition frequency, duty cycle and focal depth, sound field simulation result data S i Including each parameter combination A i The sound pressure distribution and power density obtained by simulation calculation below; Data extraction step, from each parameter combination A i The corresponding sound field simulation result data S i Extract each corresponding set of feature quantities, each set of feature quantities includes each set of parameter combination A i The sound field simulation result data S corresponding to the following i Multiple features x extracted from ij Each set of feature quantities includes at least three features of the maximum sound pressure, average sound pressure value, effective range of the sound field, sound field volume, sound field shape parameters, power density peak, temperature rise estimation, and energy absorption rate; each set of feature quantities constitutes each data point F i =(x i1 , x i2 , x ij ,…,x iM ), all data points F i Form a data set F, that is, F={F1,F2,…,F i …,F N }, where N is the number of parameter combinations, M is the number of extracted features, i is 1 to N, and j is 1 to M; A data processing step, after clustering the data set F, divides the clustering results into response levels and outputs the division results; wherein the clustering of the data set F in the data processing step specifically includes: Step 1: Density calculation. For each data point F i , calculate its local density ρ i : in, is the kernel function, which adopts Gaussian kernel function, and ε is the preset distance threshold; Step 2, calculate the density weight: Among them, ω i is the density weight, reflecting the data point F i The density of the area; Step 3, initialize cluster centers: select K initial cluster centers, denoted as H1, H2, …, H k ,…,H K , where k is 1 to K; Step 4, clustering iteration: perform the clustering iteration process until the cluster center does not change or the change is less than the set threshold, including: Step 41, cluster assignment: assign each data point F i Assign to the cluster C that is closest to it k : Step 42, cluster center update: update the cluster center according to the weighted average of the data points in the cluster: Step 43, convergence judgment: if the change of the cluster center is less than the set threshold δ, the algorithm converges and stops iterating; otherwise, continue to iterate from step 41 to step 43; Step 5: Output the clustering results. When the clustering stops iterating, output the cluster C to which each data point belongs. k And the cluster C k Center H k .

Images