Tms motor hotspot automatic search method and system fusing multi-modal priori and bayesian optimization

Abstract

This application discloses an automatic motion hotspot search method and system for TMS that integrates multimodal priors and Bayesian optimization, relating to the fields of medical robotics, neuromodulation, and medical image processing. The method includes: acquiring medical image data of the subject, extracting prior center coordinates, delineating a search space based on these coordinates, and generating a set of candidate stimulus poses; iteratively performing the search operation until a termination condition is met. Each iteration includes: based on the current Gaussian process model and the candidate set, using Bayesian optimization to determine the next target pose and calculating a comprehensive physiological efficacy score; updating the posterior distribution of the Gaussian process model using this pose and score. When the iteration terminates, the pose with the highest score in historical observations is output as the motion hotspot; if multiple poses have the same highest score, the one closest to the prior center coordinates is selected. This application achieves fully automatic and precise localization from individualized prior guidance to data-driven closed-loop optimization.

Description

Technical Field

[0001] This application relates to the fields of medical robot technology, neuromodulation technology and medical image processing technology, and in particular to a TMS motion hotspot automatic search method and system that integrates multimodal priors and Bayesian optimization. Background Technology

[0002] Transcranial magnetic stimulation (TMS), a non-invasive neuromodulation technique, works by inducing an electric field in specific cortical regions using a time-varying magnetic field, thereby modulating neuronal activity. This technique is widely used in neuroscience research (such as brain functional connectivity analysis) and clinical treatment (such as depression and post-stroke rehabilitation), and its efficacy largely depends on the precise localization of the stimulation target. Among these, the localization of motor hotspots is a crucial aspect of TMS application. A motor hotspot typically refers to a cortical location that can elicit the maximum amplitude motor evoked potential (MEP) in the target muscle with the lowest stimulation intensity. Precise localization of this hotspot is fundamental to determining an individual's resting motor threshold (RMT), which in turn is the core physiological basis for formulating all subsequent treatment stimulation intensity protocols. Due to significant differences in individual brain anatomy, hotspots are not always located in the standard "hand knot" anatomical region, thus requiring individualized searches.

[0003] Currently, the positioning of sports hotspots mainly relies on the following technical solutions: 1. Manual Search: This is the most commonly used method in current clinical practice. With the assistance of a neuronavigation system, the operator holds the TMS coil and first locates a pre-defined anatomical landmark (e.g., a hand knot) as a starting point. Then, based on experience, the operator manually moves and adjusts the coil angle near the target area on the scalp, performing trial stimulation. The operator observes the MEP waveform in real time, subjectively judges and compares the response intensity at different sites, and finally marks the location that evokes the maximum average MEP amplitude as a hotspot. This method is highly dependent on the operator's experience and feel, and is inherently subjective.

[0004] 2. Semi-automatic grid search: This method represents a higher level of standardization and is typically performed using a robot or fixed scaffold. Its core is to predefine a grid covering the representative area of ​​the target muscle in a standard brain space (such as the Montreal Neurological Institute space, or MNI space). A typical workflow employs a two-stage strategy: first, a preliminary scan is performed on a sparse grid (e.g., 10mm spacing) to identify the approximate regions with strong responses; then, using this region as the center, a finer search is performed within a denser grid (e.g., 5-7mm spacing) to exhaustively find the optimal point. While this method offers a high degree of standardization, it is essentially a blind, traversal search.

[0005] 3. Traditional Automated Search: Some cutting-edge research attempts to introduce more intelligent closed-loop optimization algorithms to reduce unnecessary sampling. These methods typically use random locations or locations based on population-averaged brain maps as the starting point for the search, and then use standard optimized acquisition functions to sample within the entire pre-defined search space, attempting to fit a model of the relationship between MEP response and stimulus location. None of these methods effectively incorporate multimodal imaging information such as functional magnetic resonance imaging (fMRI) and diffusion tensor imaging (DTI) of individual subjects as prior guidance.

[0006] In summary, existing technologies have the following inherent defects in terms of positioning mechanisms, search strategies, and evaluation criteria: 1. Lack of personalized prior knowledge and poor adaptability to pathological brains: Existing automatic or semi-automatic methods typically use random initialization or initialization based on the average position of standard healthy human brain maps. However, for subjects with stroke, brain tumors, etc., the brain functional areas often undergo reorganization or displacement, and the anatomical structure and functional location no longer correspond. This causes the algorithm to easily waste a lot of time in the early stages of the search at incorrect "standard" locations, making it difficult to quickly capture the true functional areas that have shifted, resulting in low search efficiency and even getting stuck in local optima. This fails to meet the clinical needs of special subjects, which is the "cold start" problem.

[0007] 2. High subjectivity and low repeatability: Manual search relies entirely on the operator's personal experience, hand-eye coordination, and instantaneous judgment. Because MEP (Medium Emission Point) itself exhibits significant physiological fluctuations between trials, and different operators have different habits in adjusting coil angles and positions, search results are easily misled by accidental high-amplitude MEP values ​​or affected by inter-operator differences. Studies have shown that the average deviation of hotspot locations measured by different operators or at different times can reach more than 10 millimeters, seriously affecting the consistency of treatment protocols and the repeatability of scientific research experiments.

[0008] 3. Blind and inefficient search strategy: While semi-automatic grid search has a standardized process, it is essentially an exhaustive method based on a fixed grid. To improve localization accuracy, the grid must be densified, leading to a quadratic increase in the number of stimulus points to be evaluated. The system is forced to perform ineffective tests on a large number of low-probability or unresponsive areas, failing to intelligently focus on high-potential areas. This makes a single search time lengthy (usually exceeding 20 minutes), not only reducing equipment efficiency but also easily causing subject fatigue. Fatigue itself alters cortical excitability, further interfering with measurement accuracy.

[0009] 4. The evaluation criteria are too simplistic, leading to false positives: Current technologies generally use the maximum MEP amplitude as the sole criterion for locating hotspots. However, high-amplitude MEPs may sometimes be accompanied by long latency periods (suggesting the signal transmission path may not be the most direct) or high background noise (such as artifacts generated by muscle pretension). The corresponding cortical sites may not be the physiologically optimal sites with the highest activation efficiency and the lowest required RMT. Using such false positive hotspots as the basis for developing treatment parameters may result in an overestimation of the actual required stimulation intensity, increasing the risk of discomfort and side effects in subjects. Summary of the Invention

[0010] In view of the above-mentioned defects and shortcomings of the existing technology, it is desirable to provide a TMS motion hotspot automatic search method and system that integrates multimodal prior and Bayesian optimization, which can achieve fast, accurate and adaptive fully automatic hotspot localization.

[0011] To achieve the above objectives, this application provides the following solution: Firstly, this application provides an automatic TMS motion hotspot search method that integrates multimodal priors and Bayesian optimization, including: S1: Acquire the subject's medical imaging data, extract the prior center coordinates related to the target muscle's motor function from the medical imaging data, define the search space based on the prior center coordinates, and generate a set of candidate stimulus poses based on the search space; The search iteration operation is performed repeatedly until the preset termination condition is met; Each search iteration includes: S2: Based on the current Gaussian process model and the set of candidate stimulus poses, the Bayesian optimization algorithm is used to determine the target pose of the next stimulus. Wherein, if the current search iteration operation is the first search iteration operation, the mean function of the current Gaussian process model is obtained by initialization based on the prior center coordinates; if the current search iteration operation is not the first search iteration operation, the current Gaussian process model is the Gaussian process model that updated the posterior distribution in the previous search iteration operation. S3: Control the actuator to move to the target pose, apply TMS, and simultaneously collect electromyographic signals of the target muscle; S4: Extract the multi-dimensional physiological features of the electromyographic signal, and calculate the comprehensive physiological efficacy score based on the multi-dimensional physiological features; use the target pose and the corresponding comprehensive physiological efficacy score as new observation data. S5: Use the target pose and its corresponding comprehensive physiological efficacy score as observation data to update the posterior distribution of the current Gaussian process model; When the preset termination condition is met, the search iteration operation is stopped, and the target pose with the highest comprehensive physiological efficacy score among all historical observation data is taken as the TMS motion hotspot output; if there are multiple target poses with the same highest comprehensive physiological efficacy score, the one whose spatial position is closest to the prior center coordinates is taken as the final TMS motion hotspot output.

[0012] Secondly, this application provides a TMS motion hotspot automatic search system that integrates multimodal priors and Bayesian optimization, including: a data acquisition and processing module, a model building and decision-making module, a stimulus execution and signal acquisition module, a signal processing and scoring module, a model update and iteration control module, and a hotspot output module; The aforementioned data acquisition and processing module is used to acquire the subject's medical imaging data, extract the prior center coordinates related to the target muscle's motor function from the aforementioned medical imaging data, delineate the search space based on the aforementioned prior center coordinates, and generate a set of candidate stimulus poses based on the aforementioned search space. The aforementioned model building and decision-making module is used to build a Gaussian process model and initialize the mean function of the Gaussian process model based on the aforementioned prior center coordinates during the first search iteration operation; and in each search iteration operation, based on the current Gaussian process model and the aforementioned candidate stimulus pose set, the Bayesian optimization algorithm is used to determine the target pose of the next stimulus. The aforementioned stimulus execution and signal acquisition module is used to control the actuator to move to the target position, apply TMS, and simultaneously acquire the electromyographic signals of the target muscle. The aforementioned signal processing and scoring module is used to extract the multi-dimensional physiological features of the aforementioned electromyographic signals and calculate a comprehensive physiological efficacy score based on the aforementioned multi-dimensional physiological features. The aforementioned model update and iteration control module is used to update the posterior distribution of the current Gaussian process model by taking the target pose and its corresponding comprehensive physiological efficacy score as observation data; and to control the repeated execution of the search iteration operation consisting of the aforementioned model construction and decision module, the aforementioned stimulus execution and signal acquisition module, the aforementioned signal processing and scoring module and model update function until the preset termination condition is met. The aforementioned hotspot output module is used to stop the search iteration operation when the above termination condition is met, and to output the target pose with the highest comprehensive physiological efficacy score among all historical observation data as the TMS motion hotspot; if there are multiple target poses with the same highest comprehensive physiological efficacy score, the one whose spatial position is closest to the above prior center coordinates is taken as the final TMS motion hotspot.

[0013] According to the specific embodiments provided in this application, the following technical effects are disclosed: This application provides an automatic TMS motion hotspot search method and system that integrates multimodal priors and Bayesian optimization. Through step S1 (acquiring data, extracting prior center coordinates, and generating a set of candidate stimulus poses), it solves the problems of traditional methods relying on standard brain atlases, ignoring individual anatomical variations and pathological brain remodeling, and the lack of individualized physical boundaries in robot search. It achieves a fundamental shift in the localization benchmark from a group standard to a subject-specific one, laying a unique anatomical and functional starting point for subsequent precise searches. Through iterative search operations (especially the decision-making step S2 based on the current Gaussian process model and the model update step S5), it solves the problems of low efficiency due to random initialization in traditional searches ("cold start") and the lack of quantification and active learning capabilities for uncertainty in the search process. In the first iteration, the mean function is initialized based on the prior center coordinates to achieve a "warm start" guidance; in subsequent iterations, the model is continuously updated online based on new observation data, forming a data-driven closed-loop adaptive optimization. This process can intelligently approach the global optimal target with the fewest stimulations and achieve efficient and intelligent search convergence through explicit termination conditions. By employing steps S3 (precise execution and synchronous data acquisition by the actuator) and S4 (multi-dimensional physiological feature extraction and comprehensive physiological efficacy scoring) in the search and iterative operation, the system addresses the problems of poor repeatability in manual operations, asynchronous stimulation and signal acquisition, and the susceptibility to interference and one-sidedness of traditional assessments relying on a single amplitude indicator. This achieves precise stimulus delivery, synchronous high-fidelity data acquisition, and comprehensive and robust physiological efficacy assessment based on multi-dimensional physiological features. The complete closed-loop process of "prior guidance - iterative search (perception - decision-making - execution - learning) - intelligent output" constituted by the above steps solves the systemic problems of fragmented processes, high dependence on operator experience, and poor repeatability in traditional TMS hotspot localization. Ultimately, the system outputs the optimal hotspot based on preset rules (such as the highest score and closest to the prior center), achieving personalized, automated, and standardized intelligent localization throughout the entire process, significantly improving the accuracy, efficiency, and reliability of TMS therapy in research and clinical practice. Attached Figure Description

[0014] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0015] Figure 1 This is a flowchart illustrating an automatic TMS motion hotspot search method that integrates multimodal priors and Bayesian optimization in one embodiment of this application. Figure 2 This is a schematic diagram of the discrete candidate stimulus pose set on the individualized three-dimensional head model in the embodiments of this application; Figure 3 This is a schematic diagram of multi-physiological parameter feature extraction in an embodiment of this application. Detailed Implementation

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

[0017] To make the technical solution of this application clearer and easier to understand and consult, the main technical terms involved in this application are hereby uniformly explained and described. The definitions of these terms are all based on the technical field of this application (medical execution institutions and neuromodulation technology) and the specific context.

[0018] 1. Terminology related to this application Unless the context otherwise requires, this method / system refers specifically to the "Automatic Search Method and System for Motion Hotspots in TMS Integrating Multimodal Priors and Bayesian Optimization" claimed in this application.

[0019] Transcranial magnetic stimulation (TMS): a neuromodulation technique that uses a time-varying magnetic field to non-invasively penetrate the skull and induce an electric field in a specific area of ​​the cerebral cortex, thereby regulating neuronal activity.

[0020] Motor focus: In this application, it specifically refers to the position in the primary motor cortex of the brain that can produce the "optimal comprehensive physiological efficacy score" for the target muscle under a fixed stimulation intensity. This is the basis for determining the individualized treatment target and resting motor threshold.

[0021] Motor evoked potentials: These are complex muscle action potentials recorded by electromyography on the target muscle on the contralateral side of the limb after TMS pulse stimulation of the motor cortex. They are a core physiological indicator for assessing the excitability and conduction integrity of the corticospinal tract.

[0022] Resting Motion Threshold (RMT): An important individualized physiological parameter. It is generally defined as the minimum TMS stimulation intensity (as a percentage of the device's maximum output intensity) that can elicit a peak-to-peak MEP greater than 50 μV in at least 5 out of 10 consecutive stimuli, with the target muscle completely relaxed (at rest). RMT is the benchmark for determining the intensity of subsequent treatments or experiments.

[0023] Hand knot: A typical anatomical landmark area in the primary motor cortex of the precentral gyrus of the brain that controls hand movements. Morphologically, it often presents as an "Ω" or "ε" shaped curve. In standard healthy human brain atlases, it is often used as an initial reference location for searching motor hotspots.

[0024] Actuator: refers to an automated device capable of controlled movement and precise positioning of the TMS coil. In specific embodiments, the actuator may be an industrial robot, a robotic arm, or a dedicated TMS navigation device.

[0025] Pose: A complete spatial state description including position and orientation. In this application, it specifically refers to the pose of the TMS coil stimulation focus relative to the subject's head coordinate system, typically represented by a 4×4 homogeneous transformation matrix.

[0026] 2. Medical Imaging and Navigation Related Terminology Multimodal imaging prior: refers to different types of medical imaging data from the same subject. This application includes at least T1-weighted structural magnetic resonance imaging (hereinafter referred to as T1 structural imaging), functional magnetic resonance imaging, and diffusion tensor imaging data.

[0027] Functional magnetic resonance imaging (fMRI) is an imaging technique that non-invasively maps brain function by detecting changes in blood oxygen level-dependent signals (BOLD signals) associated with neural activity. In this application, based on fMRI data, the coordinates of brain regions that are significantly activated during the performance of a target muscle movement task are extracted.

[0028] Diffusion tensor imaging (DTI) is a magnetic resonance imaging technique that non-invasively displays and tracks the trajectory of nerve fiber bundles by measuring the directionality and anisotropy of water molecule diffusion in the brain's white matter. In this application, based on DTI data, the coordinates of the corticospinal tracts that innervate the target muscles are extracted from the cortical surface where the connections are most dense.

[0029] Corticospinal tract: The main descending neural pathway in the brain that controls voluntary body movement. It originates from areas such as the primary motor cortex, descends through structures such as the internal capsule and cerebral peduncle to the spinal cord, and its integrity is the anatomical basis for the generation of MEP.

[0030] 3. Algorithm and Model Related Prior center coordinates: Three-dimensional spatial coordinates extracted from the subject's multimodal imaging data, representing brain functional areas or structural connectivity nodes that are directly or indirectly associated with the voluntary motor control or innervation of the target muscle. The brain regions corresponding to these coordinates show significant blood oxygen level-dependent (BOLD) signal activation during the execution of the target muscle movement task, or have strong anatomical connections with the descending motor fiber bundles innervating the target muscle. Peak activation points or points with the highest fiber projection density can be located from fMRI or DTI data through statistical analysis methods (such as general linear models) or fiber tracing algorithms, and mapped onto the individualized three-dimensional head model space.

[0031] Gaussian process model: The probabilistic model used in this application as a Bayesian optimization surrogate model is used to model and predict the unknown functional relationship of "stimulus pose (x) → physiological score (y)".

[0032] Comprehensive physiological efficacy score: A scalar evaluation index obtained by extracting multi-dimensional physiological features from electromyographic signals induced by a single stimulus and then weighting and fusing them. It serves as the observation value y in the Bayesian optimization algorithm process. In the score calculation part, it can also be represented as a function S(x), i.e., y=S(x).

[0033] Acquisition function: In Bayesian optimization, the criterion function used to balance "exploration" and "exploitation". This application involves two strategies: Upper Confidence Bound (UCB) and Expected Improvement (EI).

[0034] Heatmap: An auxiliary data structure corresponding to the search space, whose values ​​are updated by spreading to its spatial neighborhood based on historical observation scores, used to encode the spatial continuity prior of physiological responses.

[0035] 4. Search process related Global exploration phase: The initial stage of the search focuses on discovering potential high-response areas and employs an exploration-oriented acquisition function (such as UCB).

[0036] Local fine localization stage: After a valid response is detected, the search space is narrowed to accurately determine the optimal pose. This stage uses a local optimization-biased acquisition function (such as EI).

[0037] Candidate stimulus pose set: A finite set of all possible stimulus poses generated by discretization within a defined scalp search space.

[0038] 5. Related to physiological characteristics Amplitude: Peak-to-peak voltage of a single electromyographic signal, measured in microvolts (μV).

[0039] Latency: The time interval between the TMS pulse trigger moment and the moment when the electromyographic signal first deviates significantly from the baseline, measured in milliseconds (ms).

[0040] Waveform similarity: A quantitative index of the degree of shape similarity between a single MEP waveform and a standard MEP template waveform, calculated through methods such as normalized cross-correlation.

[0041] Signal stability: A measure of the consistency of MEP amplitude obtained from multiple stimulations of the same stimulus pose, usually calculated based on the coefficient of variation.

[0042] To make the above-mentioned objectives, features and advantages of this application more apparent and understandable, the application will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0043] This application provides an automatic TMS motion hotspot search method that integrates multimodal prior and Bayesian optimization. This method is executed by a computer device, specifically a terminal or server, or both. In this application embodiment, for example... Figure 1 As shown, the method is applied to a terminal as an example, including the following steps 1 to 6: S1: Acquire the subject's medical imaging data, extract the prior center coordinates related to the target muscle's motor function from the medical imaging data, define the search space based on the prior center coordinates, and generate a set of candidate stimulus poses based on the search space; S2: Construct a Gaussian process model, wherein the mean function of the Gaussian process model is initialized based on the prior center coordinates; based on the initialized Gaussian process model and the candidate stimulus pose set, the target pose of the next stimulus is determined using a Bayesian optimization algorithm. S3: Control the actuator to move to the target pose, apply TMS, and simultaneously collect electromyographic signals of the target muscle; S4: Extract the multi-dimensional physiological features of the electromyographic signal and calculate the comprehensive physiological efficacy score based on the multi-dimensional physiological features; S5: Use the target pose and its corresponding comprehensive physiological efficacy score as observation data to update the posterior distribution of the Gaussian process model; S6: Repeat steps S2 to S5 until the preset termination condition is met, and output the historical best pose as a TMS motion hotspot.

[0044] By implementing steps 1 to 6 above, and acquiring the subject's medical imaging data, the subject's unique multimodal images (e.g., fMRI / DTI) are fused with a high-precision head model (i.e., an individualized 3D head model). This anchors the search starting point and scope to the subject's specific brain anatomy and functional state, completely resolving the "targeting inaccuracy" problem caused by reliance on standard atlases. Furthermore, by extracting prior central coordinates related to the target muscle's motor function from the medical imaging data, a Bayesian optimization algorithm based on "hot-start" prior central coordinates is introduced to achieve dynamic search and multi-target balance, minimizing stimulation. The method efficiently approaches the global optimum, overcoming the inefficiency and blindness of traditional manual or random searches. By processing electromyographic signals and extracting multi-dimensional physiological features, a comprehensive physiological efficacy score is obtained, providing an interference-resistant and verifiable "multi-dimensional physiological fingerprint" for motion hotspots, surpassing the one-sidedness and instability of single amplitude interpretation. The method integrates image analysis, intelligent decision-making, actuator execution, signal processing, and updates of the posterior distribution of the Gaussian process model into a fully automatic "perception-decision-execution-learning" closed loop, achieving end-to-end automation from data to pose, and significantly improving operational accuracy and efficiency.

[0045] Furthermore, the multimodal image adaptive prior extraction strategy designed in this application enables the scheme to capture both functional reorganization areas and structural preservation areas, thereby covering the entire spectrum of subjects from functionally intact to severely impaired, greatly expanding the clinical applicability of precision treatment.

[0046] In addition, the aforementioned terminals can be, but are not limited to, various desktop computers, laptops, smartphones, tablets, IoT devices, and portable wearable devices. IoT devices can include smart speakers, smart TVs, smart air conditioners, smart in-vehicle devices, etc. Portable wearable devices can include smartwatches, smart bracelets, head-mounted devices, etc. Servers can be implemented using independent servers or server clusters composed of multiple servers, or they can be cloud servers.

[0047] In an exemplary embodiment, in order to construct a high-precision, individualized 3D head model and provide individualized, high-precision initial guidance and physical boundaries for subsequent automatic search, step 1: acquire the subject's medical imaging data, extract the prior center coordinates related to the target muscle motor function from the medical imaging data, delineate the search space based on the prior center coordinates, and generate a set of candidate stimulus poses based on the search space. Specifically, this includes detailed steps 11 to 13, as follows: S11: Acquire medical imaging data and reconstruct a personalized 3D head model; Step S11 specifically includes the following operations: S111, Obtain data and check: T1-weighted images, fMRI, and DTI data of the subjects were acquired, and preliminary quality checks were performed (such as checking for obvious artifacts, signal dropout, etc.). S112, Preprocessing T1 structural images and segmenting tissues: The T1 structure image is subjected to deviation field correction to eliminate grayscale deviation caused by magnetic field inhomogeneity; The T1 structural image after bias field correction was segmented into five tissues: scalp, skull, cerebrospinal fluid, gray matter, and white matter. S113, Reconstructing an individualized 3D head model: Based on the segmentation results, the T1 structural images were reconstructed from the scalp and gray / white matter respectively using the traveling cube algorithm to reconstruct a three-dimensional triangular mesh model of the subject's scalp surface, which served as an individualized three-dimensional head model. This model was then used as the visual reference and coordinate mapping carrier for subsequent neural navigation. In this embodiment, prior center coordinates related to the target muscle motor function are extracted from the medical imaging data. Its core is to integrate fMRI data with DTI data to locate the most promising physiological response points.

[0048] Activation of center coordinates based on fMRI data extraction function This method is suitable for subjects who can cooperate in performing motor tasks. During data acquisition, subjects are required to perform standardized motor tasks related to the target muscles using an MRI scanner. Task-state BOLD signals are analyzed using a general linear model to identify statistically significant activation clusters within the hand representation area or functional reorganization area of ​​the primary motor cortex (M1 area) on the affected side. The coordinates of local peak points of statistics (such as t-values) within these significant activation clusters are extracted and defined as the coordinates of the functional activation center. This coordinate system reflects the functionally active areas resulting from brain plasticity and reorganization.

[0049] Based on DTI extraction of structural connection center coordinates This method is suitable for subjects with severely impaired motor function who are unable to cooperate in completing fMRI tasks. Fiber tracing analysis is performed on DTI data. Starting from deep key nodes of the corticospinal tract on the affected side (such as the posterior limb of the internal capsule), retrograde fiber tracing is performed to the surface of the cerebral cortex. The spatial distribution density of all traced fiber ends on the cortical surface is calculated, and the geometric center coordinates of the region with the highest density are defined as the structural connection center coordinates. This coordinate represents the cortical region that maintains the strongest anatomical connection with the descending motor pathway that innervates the target muscle.

[0050] The multimodal prior fusion strategy is to obtain the coordinates of the functional activation center. Center coordinates of the structure Then, the final prior center coordinates are determined using the following strategy. If the spatial distance between the two points is less than a preset threshold, they are considered to indicate the same area. The geometric midpoint between the two points can be directly taken, or the point with the higher signal-to-noise ratio can be selected (e.g., the center coordinates of the comparison function are activated). The t-value and the coordinates of the structural connection center (Fiber connection probability value). If the spatial distance between the two is greater than or equal to the preset threshold, the selection should be made according to the clinical focus: if the focus is on functional recovery, the coordinates of the functional activation center should be selected first. If the focus is on structural connections, then the coordinates of the structural connection center should be selected first. .

[0051] S12: Extract the prior center coordinates of the multimodal image; Step S12 specifically includes the following operations: S121, Activate center coordinates based on fMRI data extraction function (Applicable to subjects who can cooperate with the task, which requires subjects to perform movements relative to the TMS target muscle at a constant frequency during fMRI acquisition (i.e., task design)), specifically including the following steps: A1) After head motion correction and slice time correction, the fMRI data were spatially smoothed using a Gaussian kernel (e.g., 6mm FWHM), and the preprocessed fMRI data were registered with high-resolution T1 structural images. A2) For the time series of each voxel in the preprocessed fMRI data, a general linear model is constructed using the following formula to predict the activation intensity of each voxel. And calculate its corresponding Statistics, generating statistical graphs: in, The observed BOLD signal matrix; An orthogonal design matrix containing task design regressors and covariates such as head movement parameters; This is the error term; A3) On the subject's T1 structural image, manually delineate the hand representative area of ​​the primary motor cortex (M1 area) on the affected side (e.g., manually delineate the primary motor cortex area along the central sulcus and anterior central sulcus, and limit it to the area near the hand knuckles facing the hand), generate the region of interest (i.e. M1-ROI) of the primary motor cortex, and convert it to fMRI space. A4) Within the M1-ROI transformed to fMRI space, perform cluster-level multiple comparison correction on the statistical map (e.g., first set a voxel-level threshold p < 0.001, then perform cluster-level FWE correction p < 0.05). In the significantly activated clusters obtained after correction, locate the local peak points of the t-statistic using the following formula, and define the MNI coordinates obtained by registration and normalization transformation as the functional activation center coordinates. This function activates the center coordinates It reflects the functionally active areas after brain plasticity reorganization.

[0052] in, Activation intensity of independent variables in the regression model The corresponding regression coefficients; Regression coefficients The standard error of the estimate; S122, Based on DTI data extraction, connect the center coordinates of the structure. ( Specifically, it includes the following steps: B1) Perform head motion and eddy current correction preprocessing on the DTI data, and then register the preprocessed DTI data (such as the b0 image) with the high-resolution T1 structural image (rigid or affine). B2) For each voxel in the preprocessed DTI data, calculate its diffusion tensor using the following formula. And the anisotropy fraction FA: in, , , For diffusion tensor eigenvalues, ≥ ≥ ; B3) Preset the turning angle threshold (<45°) and the termination threshold of the anisotropy fraction FA (>0.2); use a fiber tracing algorithm (such as a probabilistic tracing algorithm) to select seed points from the predetermined deep structures of the corticospinal tract on the affected side (such as the posterior limb of the internal capsule) and perform reverse fiber tracing to the cortical surface; project the traced fiber ends to the cortical surface, calculate the projection area of ​​all traced fiber ends on the cortical surface, and then use the geometric center of the projection area or the weighted center with the highest fiber density as the coordinate of the structural connection center. This structure connects the central coordinates. This represents the cortical region that maintains the strongest anatomical connection with the descending motor pathway.

[0053] S123, Calculate the center coordinates. Center coordinates of the structure In this context, the coordinate point with the highest potential for physiological response is defined as the prior center coordinate. : If the function is activated, the center coordinates Center coordinates of the structure If the Euclidean distance is less than 10mm, the midpoint of the line connecting them can be directly taken, or a higher signal-to-noise ratio can be selected (e.g., the center coordinate of the comparison function is activated). The t-value and the coordinates of the structural connection center The center of the larger of the fiber density and connection probability values ​​is used as the final prior center coordinate. ; If the function is activated, the center coordinates Center coordinates of the structure If the Euclidean distance is greater than or equal to 10 mm, the selection should be made based on clinical assumptions. For example, if the focus is on functional recovery, the coordinates of the functional activation center should be selected first. As the final a priori central coordinates If the focus is on structural connections, then the coordinates of the structural connection center should be selected first. As the final a priori central coordinates ; S13: Discretizing the search space and constructing the candidate stimulus pose set; Step S13 specifically includes the following operations: S131, Calculate the prior center coordinates obtained in step S12. The shortest Euclidean distance to the individualized 3D head model reconstructed in step S11 is used to determine the coordinates of the individualized 3D head model from the prior center. The coordinates of the nearest vertex are defined as the coordinates of the scalp mapping point; S132, construct a tangent plane centered on the scalp mapping point coordinates; on the tangent plane, based on a preset range ( A two-dimensional grid is generated within a rectangular area; the two-dimensional grid is projected and mapped onto the surface of the individualized three-dimensional head model, and used as a global exploration space. The global exploration space is used to define the range of the scalp that the actuator end may contact (i.e., the activity boundary). S133, the global exploration space is discretized into a grid with a preset spacing (3mm) to obtain N discrete scalp surface locations. ,in , , Let N be the three-dimensional coordinates (i.e., translation vector) of the point in space, N≥1, i∈N; for each discrete point Calculate the external normal vector of its scalp surface. The normal vector here It can be obtained directly from the geometric properties of the triangular mesh (such as vertex normal or face normal interpolation). S134, for each discrete point location Construct its pose rotation matrix : C1) Determine the Z-axis: (Pointing to the inside of the brain); C2) Construct the initial X-axis: Take the X-axis vector of the world coordinate system Project it onto a plane perpendicular to On the tangent plane, obtain the intermediate quantity = Then normalize it to get = ; C3) Calculate the Y-axis: = ; C4) Adjust the handle direction: Calculation requires bypassing The angle of rotation θ (θ=45°) will and Around Rotate by an angle θ, so that the rotated... The projection of the (handle direction) onto the cutting plane is aligned with the preset anatomical direction (e.g., [-1, -1, 0] for the left hemisphere) to obtain the final result. and ; C5) Rotation Matrix Each column represents the unit direction vector of that axis in the world coordinate system; S135, the position of each discrete point Its rotation matrix By combining the following formulas, a complete TMS coil pose matrix is ​​generated. ; in, ~ These are the rotation / scaling / shearing submatrices of the voxel in the initial tangent plane coordinate system, respectively. The pose matrix generated by combining all discrete scalp surface position coordinates with their rotation matrix. The set of these can be used as the candidate stimulus pose set to be stimulated. See Figure 2 ; In an exemplary embodiment, step 2: constructing a Gaussian process model, wherein the mean function of the Gaussian process model is initialized based on the prior center coordinates; based on the initialized Gaussian process model and the candidate stimulus pose set, the target pose of the next stimulus is determined using a Bayesian optimization algorithm, specifically including a refinement of steps 21 to 23, as follows: S21: Define the Gaussian process model and initialize the mean function; Step S21 specifically includes the following operations: S211, the input variable x of the Gaussian process model is defined as the spatial coordinates of the cortical observation point to be stimulated, and the output variable y is the comprehensive physiological efficacy score obtained after stimulation of that point, denoted as... ; Assume the observed score y is composed of the true physiological response f(x) at that point superimposed with Gaussian observation noise ε, i.e., the observation model is: in, To observe the noise variance; Furthermore, assume that the actual physiological response f(x) to be determined follows a Gaussian process prior: in, The point to be observed The mean function at; The point to be observed and its neighboring points The covariance kernel function is used to describe the observation points. and its neighboring points The similarity between physiological responses; in this embodiment, the covariance kernel function Select the radial basis kernel function: in, It is a natural exponential function; The signal variance controls the fluctuation range of the response amplitude. The point to be observed and its neighboring points The Euclidean distance between them; The length scale controls the smoothness of the function's variation, and its physical meaning corresponds to the approximate spatial range of "hot spots" on the cortex. S212, Initialize the mean function of the Gaussian process model based on the prior center coordinates; As a key improvement of this application, the mean function Instead of a traditional zero-mean function, it is constructed as a function with the aforementioned prior center coordinates. The core is a Gaussian bias function: in, The preset confidence level range; The radius of the prior confidence range; The mean function Before the iteration begins (when there is no observation data), at the prior center coordinates This creates a high expectation value, guiding the Bayesian optimization algorithm to directly target the high-potential region with the first stimulus. This mechanism effectively utilizes individualized multimodal imaging prior information (i.e., fMRI and DTI data) to achieve precise "warm start" for pathological brains or individual differences, overcoming the inefficiency of "cold start" caused by random or standard atlas initialization in traditional methods.

[0054] S22: Calculate the posterior distribution; Step S22 specifically includes the following operations: S221, Collect historical observation data Assuming the current task is completed The second stimulus collects its historical observation dataset, denoted as... ; in, For the first The pose coordinates of each stimulus point; for position Corresponding rating; Let the observation input matrix Observation vector ; S222, Calculate the posterior distribution For the candidate stimulus pose set generated in step S13 In the dataset, for any unobserved point x, in the known historical observation dataset Under the condition that, its posterior distribution It is still a Gaussian distribution: = in, The predicted mean represents the expectation of obtaining a high score at the observed point x; The prediction variance represents the uncertainty in predicting the observation point x. Based on the covariance kernel function defined in step S211 and the mean function defined in step S212 The predicted mean is calculated using the following formula. and prediction variance : in: = A vector containing the mean function values ​​at known historical observation points; = Let x be the covariance vector, representing the similarity between the observed point x and each historical observation point; For each historical observation point and the point to be observed Similarity between them; Let be a t×t covariance matrix, representing the pairwise similarity between historical observation points, whose elements are... The matrix actually used in the posterior calculation is + Where I is the identity matrix, The observation noise variance is defined in S211.

[0055] To observe the residual vector; S23, Bayesian decision-making and sampling; Step S23 specifically includes the following operations: S231, Phase Monitoring and Switching Real-time monitoring of two state variables: the number of times a better score was not found in consecutive iterations. and the current best overall physiological efficacy score ; If both conditions are met and Then, the process switches from the global exploration phase to the local fine-tuning phase; Otherwise, remain in the global exploration phase; in, This is the global stage stagnation threshold; This is the physiological activation threshold; S232, Update heatmap matrix Each time at the stimulation point Get a rating Then, based on the prior knowledge of spatial continuity, the stimulus point is... Centered on, with a preset radius (For example Any grid point within ) The Gaussian weighted update is performed using the following formula: in, The learning rate; The length scale of the heatmap is (e.g., 3 mm). Stimulus point and its neighboring points The Euclidean distance between them; This operation creates "hot zones" around high-response points, thereby guiding the algorithm to prioritize exploring these potentially high-probability areas in subsequent decisions.

[0056] S233, Space Adaptive If we are currently in the local fine-tuning stage, then we use the historical best point. Define the shrunken search radius centered on the search center. : in, Minimum precision radius; This is the original search radius during the global exploration phase; Defined as the current historical best point Centered on, with The spherical sub-region is designated as the new region of interest (ROI); the candidate stimulus pose set generated in step S13 is... The cropped set of sub-candidate stimulus poses is obtained by cropping the poses to include only those located within the new region of interest (ROI) using the following formula. : in, Here is the pose matrix; It is a distance function; For mapping functions; If we are currently in the global exploration phase, we can directly use the candidate stimulus pose set. ; S234, Calculate the acquisition function Extracting the candidate stimulus pose set (the candidate stimulus pose set used in the global phase) Local phases using sub-candidate stimulus pose sets Generate a candidate coordinate set from all candidate coordinates in the given set; calculate the acquisition function for each spatial coordinate point x in the candidate coordinate set based on the current search stage. : In the global exploration phase, a confidence upper bound strategy based on fused heatmaps is adopted to calculate the candidate stimulus pose set generated in step S13. In the equation, the expected utility of each candidate spatial coordinate point x is... : in, To explore the weighting coefficients; These are the weighting coefficients for the heatmap; for Standard deviation at point; In the local fine localization stage, an expectation boosting strategy is adopted, and step S233 is used to obtain the set of sub-candidate stimulus poses. In the equation, the expected utility of each candidate spatial coordinate point x is... : in, The cumulative distribution function of the standard normal distribution; The probability density function is the standard normal distribution. For fine-tuning parameters; This strategy focuses on finding the pose with the greatest potential to surpass the current best point within its neighborhood.

[0057] S235, Calculate the acquisition function for all coordinate points in the candidate coordinate set used in the current stage. Value, select to enable the acquisition function The coordinates of the point with the largest value will be used as the target coordinates for the next stimulus. Output: in, For candidate stimulus pose set or sub-candidate stimulus pose set The set of spatial coordinates corresponding to all poses; Current expected utility or expected utility ; Based on the selected coordinates The corresponding preset pose matrix is ​​retrieved from the pre-generated candidate stimulus pose set and used as the next target pose to be stimulated to drive the actuator action. Output.

[0058] In another exemplary embodiment of this application, step 3: control the actuator to move to the target pose, apply TMS, and simultaneously acquire electromyographic signals of the target muscle. The actuator can be an industrial robot, a robotic arm, or a dedicated TMS positioning device; in this embodiment, a robot is used. Step 3 specifically includes detailed steps 31 to 33, as follows: S31, coordinate system transformation and target pose calculation; step S31 specifically includes the following operations: S311, Receive target pose Receive the target pose defined in the subject's head image coordinate system as output by the decision in step S235. To drive the robot to the target pose defined by the above image coordinate system. It needs to be converted to the robot's base coordinate system; S312, call the transformation matrix Call the pre-calibrated hand-eye calibration matrix and tool center point matrix ; Among them, hand-eye calibration matrix This describes the pose of the origin of the image coordinate system in the robot's base coordinate system; the tool center point matrix. , used to describe the pose of the TMS coil stimulation focus (usually the coil center) in the robot end flange base coordinate system; S313, Calculation Target pose: Through chain-link coordinate transformation, the calculation is performed to ensure that the stimulation focus of the TMS coil reaches the target pose. The corresponding target pose that the robot end effector flange should achieve in the base coordinate system. : in, To transform the target pose of the TMS coil from the image coordinate system to the robot base coordinate system; Tool center point matrix The inverse matrix; However, what is directly controlled by the robot joint is its end effector flange, not the coil itself. The robot end effector describes the pose of the stimulation focus of the TMS coil relative to the robot end effector flange. Therefore, to obtain the target pose of the flange in the base coordinate system, the influence of the TMS coil must be "subtracted," i.e., multiplied by the tool center point matrix. inverse matrix .

[0059] Ultimately, the target pose This refers to the target pose of the end effector flange, described in the base coordinate system, required to drive the robot's movement. When the robot flange precisely reaches the target pose... This ensures that the stimulation focus of the TMS coil is strictly located in the image space specified by the time. Position, and maintain the optimal stimulus pose required by the algorithm.

[0060] This formula ensures that when the robot flange reaches the target pose... At that time, the stimulation focus of the TMS coil will be precisely located at the target pose specified in the image coordinate system. Position and maintain the specified pose.

[0061] S32, Inverse kinematics solution and trajectory motion; Step S32 specifically includes the following operations: S321, the target pose calculated in S313 The robot is input into an algorithm, and it uses inverse kinematics to calculate the target pose. The required vectors of joint angle values ​​are used as the target configuration q. S322, Based on the current joint configuration of the robot and the target joint configuration q, plan the robot's motion trajectory; control the robot to move along the planned trajectory to the target pose. This ensures that the TMS coil maintains a preset tangential contact with the subject's scalp; S33, Synchronous stimulus triggering and multiple trial signal acquisition; Step S33 specifically includes the following operations: S331, TMS parameter settings and triggering The TMS stimulator is controlled to apply multiple TMS pulses (usually 3-5 times) at random intervals within a first preset time window (e.g., 3-5 seconds) with a preset stimulation intensity (usually based on an initial estimate or adaptive adjustment) to obtain a statistically significant response. S332, Synchronous Acquisition and Segmentation Each time a TMS pulse is emitted, a high-precision TTL synchronous trigger signal is simultaneously generated and sent to the electromyography amplifier. After receiving the TTL synchronous trigger signal, the electromyography amplifier automatically extracts a segment of the original electromyography signal centered on the stimulation trigger moment (e.g., 50ms before trigger to 100ms after trigger) as an independent single segment of electromyography signal. S333, Data Encapsulation After completing a preset number of stimulation-acquisition cycles (e.g., 5 times), all generated single-segment electromyographic signals are packaged into a multi-trial data package. In another exemplary embodiment of this application, step S4: extracting multi-dimensional physiological features of the electromyographic signal and calculating a comprehensive physiological efficacy score based on the multi-dimensional physiological features, specifically includes detailed steps 41 to 43, as follows: S41: Single-trial signal feature extraction; Step S41 specifically includes the following operations: To avoid phase cancellation caused by direct waveform averaging, a strategy of "feature extraction first, followed by statistical analysis" is adopted. This involves collecting single-segment electromyographic signals from each stimulus in a multi-trial data packet. The following processes are performed respectively: S411, Amplitude Extraction: In a single segment of electromyography signal Within the primary response time window (e.g., 20ms to 50ms), calculate its peak-to-peak value (P-value). (), used to reflect the intensity of the stimulus response in that trial; S412, Extracting Latency: Selecting a Single Segment of Electromyography Signal The time difference between the moment when more than G consecutive sampling points deviate from the resting baseline (e.g., the amplitude consistently exceeds ±2 standard deviations of the resting baseline) and the moment of stimulus triggering is taken as the latency of that trial. (), used to reflect the nerve conduction velocity of that trial; S413, Calculate waveform similarity: Call the pre-stored standard MEP waveform template T, and calculate the single-segment electromyography signal of the current trial using the following formula. The normalized cross-correlation coefficient between the waveform and the template T is used as the waveform similarity ( ): Where Cov is the covariance; The standard deviation of template T; Waveform of a single-segment electromyographic signal Standard deviation; S414, Validity Determination Set waveform similarity threshold (like =0.6), determine the single-segment electromyographic signal of the current trial. Is it effective? like( )≥ If so, the trial is deemed valid; Otherwise, the trial is deemed invalid and removed or its weights are reset to zero in subsequent analyses; S42, Trial Statistical Analysis and Stability Assessment; Step S42 specifically includes the following operations: S421, if all trials in the multi-trial data packet are deemed invalid, the current stimulus pose is determined to be invalid, and its comprehensive physiological efficacy score S(x) can be set to 0 or a very low value, and directly jump to S44 feedback. Otherwise, calculate the average amplitude of all trials. and average incubation period : S422, Calculate the coefficient of variation for all valid trials: Wherein, CV is the coefficient of variation, used to evaluate the stability of the signal; The standard deviation of the amplitude for all valid trials; S423, Define the signal stability score based on the coefficient of variation (CV). ; Specifically, the signal stability score The coefficient of variation (CV) is mapped to a stability measure between 0 and 1. A score of 1 is awarded when the response is perfectly consistent (CV=0), and a score of 0 is awarded when the fluctuation exceeds the mean (CV≥1). It is used to reward high repeatability, and its calculation formula is as follows: When the coefficient of variation (CV) is smaller (the response is more consistent), the stability score (STA) is closer to 1 (most stable); when CV ≥ 1, STA = 0, which means the response is extremely unstable.

[0062] S43, Multi-parameter normalization and construction of comprehensive physiological efficacy score; Step S43 specifically includes the following operations: Features with different physical units and dimensions are integrated into a single comprehensive physiological efficacy score to reflect the intensity and reliability of hotspots. S431, using the mean and standard deviation of the average amplitude / latency corresponding to all attempted stimulus poses up to the current moment, calculate the average amplitude for the current pose. and average incubation period Dynamic Z-score standardization is performed to eliminate dimensions, giving it zero mean and unit variance, thus obtaining the standardized amplitude. and incubation period ;amplitude Used to reward high reaction intensity; latency period Used to punish those with long incubation periods; S432 constructs a comprehensive physiological efficacy score using a weighted summation model. : in, The average waveform similarity of valid trials is used to reward waveform standards and eliminate artifacts. w1, w2, w3, and w4 are preset weighting coefficients for each physiological dimension, used to balance the influence of each physiological dimension on the final score. Specifically, the weighting coefficients can be set according to clinical experience (e.g., w1=0.4, w2=-0.3, w3=0.2, w4=0.1, where w2 is a negative value to reflect the penalty for long latency). This model rewards high amplitude, short latency, high waveform fidelity, and high stability.

[0063] Figure 3 This diagram illustrates the extraction of multiple physiological parameters, showing the differences and key feature parameters between the "current test waveform" and the "standard template waveform." The solid line represents the actual electrical signal waveform being tested, while the dashed line represents the preset "standard template waveform" (for comparison). The vertical axis represents voltage (in μV), reflecting the signal amplitude; the horizontal axis represents time (in ms), showing the signal's change over time. Figure 3 As can be seen, during the latency period (Vlat) of approximately 0 to 20 ms, the signal amplitude remains relatively stable near 0 with minimal fluctuations, representing the "stable period before signal triggering." The feature extraction window (20 to 50 ms) is the region where the signal begins to fluctuate dramatically after 20 ms, and is the core area for analyzing signal characteristics (such as amplitude and waveform shape). The peak voltage (Vamp) is the point of maximum signal amplitude within the feature window, reflecting the strength of the current test signal (the peak values ​​of the current waveform and the template waveform differ in the figure). The waveform similarity (Correlation) occurs after 50 ms, when the signal tends to stabilize. The current test waveform has a slightly lower peak amplitude than the standard template and a slightly delayed latency. The waveform similarity ρ value is 0.87, indicating that the two are highly similar in morphology but have a slight phase shift. Combining the standardized amplitude and latency indices, as well as the stability score STA=0.91 (CV=0.09), the comprehensive physiological efficacy score S(x)=0.83 is calculated, suggesting that the stimulation site has strong physiological response efficacy and good repeatability. This score can be dynamically updated and mapped to the individualized three-dimensional head model to form a real-time hotspot efficacy distribution map, which facilitates the precise location of the optimal stimulation target during surgery.

[0064] (From left to right along the timeline) S44, Closed-loop feedback execution; Step S44 specifically includes the following operations: S441 will integrate physiological efficacy scores Define the observed value y corresponding to the current stimulus pose x; S442, binds the data pair (x,y) to the new observation data; S443, output the new observation data (x,y) to step S5 (update the posterior distribution of the Gaussian process model) to complete the learning and update of this iteration; In an exemplary embodiment, in order to integrate new data into the historical observation set and efficiently update the Gaussian process model's predictions (mean and variance) for all candidate poses, step S5: uses the target pose and its corresponding comprehensive physiological efficacy score as observation data to update the posterior distribution of the Gaussian process model, specifically including refining steps 51 to 55, as follows: S51: Receiving and binding new observation data Receive output (new observation data (x,y)), and then the points of the new stimulus. and its corresponding comprehensive physiological efficacy score and bind it as a new data pair , ); S52, Expand the historical observation dataset; Step S52 specifically includes the following operations: S521, will put new data into... , Added to historical observation dataset In this process, an updated dataset is formed. : ∪( , ) S522, Update the observation input matrix and observation vector Obtain the new input matrix and new observation vector ; S53, Update model parameters; Step S53 specifically includes the following operations: To avoid repeatedly calculating the inverse of a large matrix during each prediction, this step is based on the new input matrix. and new observation vector Calculate and cache the following intermediate values ​​that depend only on historical observation point X: S531, according to the radial basis kernel function defined in step S211 The covariance matrix K among all historical observations (including new data pairs) is calculated using the following formula: in, ∈ , ∈ ; S532, according to the observation noise variance defined in step S211 Construct matrix A for posterior calculation: in, (t+1)× The identity matrix of t+1); S533, calculate the lower triangular matrix L of matrix A by Cholesky decomposition such that A = LL. T And solve for vector α using the following formula: Lα= in, For the mean function in A vector formed by the values ​​of each point; The triangular matrix L and the vector α are cached, which contain all the information needed to compute the posterior prediction for any point.

[0065] S54, Define the updated posterior prediction function; Step S54 specifically includes the following operations: Based on the triangular matrix L and vector α cached in step S53, for any point Define its updated posterior prediction function: S541, the covariance vector is calculated using the following formula: = S542, the predicted mean is calculated using the following formula: S543, according to get The prediction variance is calculated using the following formula: This prediction variance reflects the variance after adding new data for any point. Reduced uncertainty in forecasting.

[0066] S55, triggering the next decision cycle. The updated posterior prediction function (i.e. and Set it to the available state and start the next round of acquisition function calculation and pose decision.

[0067] In an exemplary embodiment, to intelligently decide whether the search process should continue or terminate after each iteration (completing one cycle of steps S2 to S5), and to validate the final result to ensure that the output is a reliable target with clinical significance, step S6: Repeat steps S2 to S5 until a preset termination condition is met, and output the historical optimal pose as a TMS motion hotspot; during the cyclic execution of S2 to S5, monitor the optimization status in real time, and determine whether the search has converged based on preset multi-dimensional criteria. Once the termination condition is met, stop the search, and output the optimal stimulation target found to date with the strongest physiological evidence. Step S6 specifically includes detailed steps 61 to 63, as follows: S61, Preset multi-dimensional termination conditions Before the search begins, the following two types of judgment conditions are preset: Limit termination (hard safety constraint): If the total number of iterations t reaches the preset maximum limit of iterations T max If the search is terminated immediately to prevent an infinite loop and ensure that the operation time is controllable, this is the highest priority termination condition. Process convergence and termination (efficiency and performance assessment): Global exploration phase: If continuous The second iteration did not improve the historical best overall physiological efficacy score. ,and ≥ (If a valid point has been found), it will trigger a switch from the global exploration phase to the local fine-tuning phase. This process itself does not directly terminate the search. Local fine-tuning stage: If during this stage, continuous Second-rate( < Iterations have not yet improved the historical best overall physiological efficacy score. If the search is successful, it is determined that the optimization has been fully completed in the local area, the performance is saturated, and convergence termination is triggered.

[0068] The key parameters in the above termination conditions can be determined based on clinical experience and experimental data. Typical values ​​or determination methods are as follows: Maximum number of iterations T max The number of searches is usually set between 20 and 50 to ensure that the total time for a single search does not exceed 15-30 minutes. Global stage stagnation threshold Typically 3 to 6 times; Local stage stagnation threshold : Usually 2 to 4 times, and satisfy < ; Final physiological effectiveness threshold Its calibration method is the same as the physiological activation threshold. However, the values ​​are higher, typically ranging from 0.6 to 0.8 after normalization of the overall score.

[0069] Specifically, physiological activation threshold This refers to a pre-defined comprehensive physiological efficacy score scalar value used during the search process to determine whether a target muscle motor cortex activation response with clear physiological significance has been found, thereby triggering a switch in the search strategy from the global exploration phase to the local fine localization phase. When the comprehensive physiological efficacy score S(x) of a certain stimulus pose first reaches or exceeds the physiological activation threshold... When the pose has elicited a sufficiently strong, stable, and waveform-consistent motor evoked potential (MEP), it indicates that the pose is within the effective neighborhood of the target muscle's motor functional area, and subsequent searches should focus on this region for fine-tuning. Typically, the distribution of composite scores obtained from known motor hotspots in historical healthy subjects or subject groups is statistically analyzed (e.g., taking the 20th to 40th percentile); or through preliminary experiments, stimulation is applied to representative cortical areas of the target muscle, and the lowest composite score corresponding to a MEP with a stable amplitude >50µV, latency within the normal physiological range, and typical waveform morphology is defined as the physiological activation threshold. With the comprehensive physiological efficacy score S(x) normalized to the [0,1] interval, the physiological activation threshold is... It is usually set between 0.4 and 0.6.

[0070] Final physiological effectiveness threshold It is the final output threshold, used to determine whether the search has successfully found a "hotspot" that meets the requirements of clinical treatment.

[0071] S62, Iterative Loops and Real-time Monitoring After each step S5, which updates the posterior distribution of the Gaussian process model, the following judgment is immediately executed: S621, update iteration t=t+1; S622, Assessment Termination Conditions: Determine whether the current state satisfies any of the preset termination conditions in step S61: If none of the termination conditions are met, the process returns to step S2, and the next round of the "decision-stimulus-evaluation" cycle begins based on the latest updated Gaussian process model. If any termination condition is met, exit the loop and execute step S63; S63, generates the final hotspot When the algorithm terminates, a structured TMS motion hotspot localization report is automatically generated and output: S631, Read the final state: Obtain the historical best pose and its historical best comprehensive physiological efficacy score ; S632, Validity Verification: like ≥ When the search is successful, the current optimal pose is determined. As the final TMS motion hotspot output; Otherwise, the search will be considered a failure, and a message will be displayed. S633, output the multi-dimensional decomposition features of this pose (mean amplitude, mean latency, mean waveform similarity and stability score corresponding to this pose). S64, Process Termination and Data Archiving Send a command to move the robot end effector away from the subject's head to a preset safe position, and save all raw data, intermediate results and the "TMS Motion Hotspot Location Report" generated by S63 for this search session in a standardized format (such as BIDS) to ensure that the process is traceable and the results are reproducible.

[0072] This application addresses the "anatomical-functional mismatch" problem caused by traditional methods relying on standard atlases or random initialization by introducing multimodal imaging priors (fMRI / DTI) to extract subject-specific functional areas as the search starting point, significantly improving adaptability to pathological brains (such as post-stroke remodeling). Based on this, the system employs a two-stage adaptive search strategy of "global exploration (UCB + heatmap) and local fine localization (EI + spatial contraction)," greatly improving localization efficiency while balancing search breadth and accuracy, overcoming the time-consuming nature of grid search and the algorithm's susceptibility to local optima. To ensure the physiological authenticity and clinical reliability of the located points, this invention constructs a multidimensional physiological efficacy evaluation model integrating "amplitude, latency, waveform similarity, and signal stability," replacing the traditional single amplitude standard and effectively eliminating false positives. Ultimately, the entire process of "image registration - intelligent decision-making - robot execution - signal evaluation - model learning" has achieved full closed-loop automation. Driven by algorithms and relying on the high-precision execution of robots, it completely eliminates the subjective errors and inconsistencies of human operation, providing a highly objective, standardized and repeatable precise positioning method for clinical research and treatment.

[0073] This application also provides an application scenario in which the above-mentioned automatic TMS motion hotspot search method integrating multimodal priors and Bayesian optimization is applied. Specifically, the automatic TMS motion hotspot search method integrating multimodal priors and Bayesian optimization provided in this embodiment can be applied in the rehabilitation treatment of stroke subjects. A subject with left middle cerebral artery infarction and right-sided hemiplegia requires TMS treatment to promote motor function recovery, but traditional manual search methods are difficult to accurately locate motion hotspots after functional reorganization. The subject's T1 structural image, fMRI data, and DTI data are imported into the system after quality control, and then registered, segmented, and reconstructed to generate an individualized three-dimensional head model that can be used for navigation.

[0074] Step 1: Extraction of prior center coordinates and construction of search space - The system automatically analyzes fMRI data, identifies significant activation clusters in the left primary motor cortex hand region, and extracts the coordinates of functional activation centers. (Coordinates: MNI[-36,-22,58]); Simultaneously, based on DTI data, reverse fiber tracing was performed to extract the cortical points with the densest connections to the corticospinal tract as the coordinates of the structural connection center. (Coordinates: MNI[-38,-20,60]); The distance between the two is approximately 3mm, and the system automatically takes the midpoint as the prior center coordinates. (MNI[-37,-21,59]); based on prior center coordinates The point mapped onto the scalp surface is used as the geometric center, and a preset rectangular search space is defined as the region of interest. The region of interest is discretized on the scalp surface into 196 candidate stimulus poses (grid spacing 3mm), and each pose includes three-dimensional coordinates and coil normals.

[0075] Step 2: Hot-start Bayesian optimization algorithm search Initialize the Gaussian process model, with the mean function using the prior center coordinates. A Gaussian bias is set at the center to achieve a "hot start"; the system enters the global exploration phase and uses the UCB acquisition function with fused heatmap to select the first stimulus point.

[0076] Step 3: Automatic Stimulation and Signal Acquisition by the Robot The robot automatically moves to the target pose, with the coil fitting against the scalp; the TMS applies 5 pulses at 70% of its maximum output intensity (with an interval of 4±0.5 seconds); the electromyography system simultaneously acquires the electromyographic signals of the right first interosseous dorsal muscle, and extracts data from -50ms to +100ms after each stimulation.

[0077] Step 4: Multidimensional physiological efficacy scoring Five electromyographic (EMG) signals were analyzed to extract the following: amplitude (peak-to-peak value), latency, and waveform similarity to a standard template. The average amplitude was calculated to be 1.2 mV, the average latency to be 24.3 ms, the average waveform similarity to be 0.82, and the signal stability score to be 0.75. A comprehensive physiological efficacy score was also calculated. x) = 0.68; Step 5: Model Update and Iterative Decision Making The pose and score are added to the observation dataset, and the posterior distribution of the Gaussian process is updated. The model prediction shows that the high-potential region is concentrated in the front left of the prior center coordinate. The system continues to select the next stimulus pose and repeats steps 3 to 5.

[0078] Step 6: Stage Switching and Convergence After the 8th iteration, the optimal comprehensive physiological efficacy score failed to improve by 3 for three consecutive iterations, and the optimal comprehensive physiological efficacy score... =0.71> =0.5, the system switches to the local fine positioning stage; the search space shrinks to a 9mm radius area around the current optimal point, and the acquisition function is changed to Expected Improvement (EI); after 4 more iterations, the score fails to improve for 2 consecutive times in the local stage, and the termination condition is met.

[0079] Output results: Final hotspot coordinates: MNI[-35,-23,57] (approximately 3mm from the initial prior center coordinates), optimal comprehensive physiological efficacy score 0.73, corresponding physiological characteristics: average amplitude: 1.35mV, average latency: 23.8ms, waveform similarity: 0.85, stability score: 0.78; Report: Total iterations: 12; Total time: approximately 5 minutes (including robot movement and stimulus intervals); Motion hotspots meeting the final physiological effectiveness threshold were successfully located and are recommended as target points for subsequent resting-state motion threshold measurement.

[0080] This application example demonstrates the complete workflow of the TMS motion hotspot automatic search system, which integrates multimodal priors and Bayesian optimization, in the rehabilitation treatment of stroke subjects. By integrating multimodal image priors and adaptive Bayesian optimization algorithms, the system achieves rapid, accurate, and automatic localization of individualized and reconstructed motor functional areas, providing reliable technical support for personalized and precise TMS treatment.

[0081] Based on the same inventive concept, this application also provides a system for implementing the above-mentioned automatic TMS motion hotspot search method that integrates multimodal priors and Bayesian optimization. The solution provided by this system is similar to the implementation described in the above method. Therefore, the specific limitations of one or more embodiments of the automatic TMS motion hotspot search system integrating multimodal priors and Bayesian optimization provided below can be found in the limitations of the automatic TMS motion hotspot search method integrating multimodal priors and Bayesian optimization described above, and will not be repeated here.

[0082] In an exemplary embodiment, a TMS motion hotspot automatic search system integrating multimodal priors and Bayesian optimization is provided, including: a data acquisition and processing module, a model building and decision-making module, a stimulus execution and signal acquisition module, a signal processing and scoring module, a model update and iteration control module, and a hotspot output module; The aforementioned data acquisition and processing module is used to acquire the subject's medical imaging data, extract the prior center coordinates related to the target muscle's motor function from the aforementioned medical imaging data, delineate the search space based on the aforementioned prior center coordinates, and generate a set of candidate stimulus poses based on the aforementioned search space. The aforementioned model building and decision-making module is used to build a Gaussian process model and initialize the mean function of the Gaussian process model based on the aforementioned prior center coordinates during the first search iteration operation; and in each search iteration operation, based on the current Gaussian process model and the aforementioned candidate stimulus pose set, the Bayesian optimization algorithm is used to determine the target pose of the next stimulus. The aforementioned stimulus execution and signal acquisition module is used to control the actuator to move to the target position, apply TMS, and simultaneously acquire the electromyographic signals of the target muscle. The aforementioned signal processing and scoring module is used to extract the multi-dimensional physiological features of the aforementioned electromyographic signals and calculate a comprehensive physiological efficacy score based on the aforementioned multi-dimensional physiological features. The aforementioned model update and iteration control module is used to update the posterior distribution of the current Gaussian process model by taking the target pose and its corresponding comprehensive physiological efficacy score as observation data; and to control the repeated execution of the search iteration operation consisting of the aforementioned model construction and decision module, the aforementioned stimulus execution and signal acquisition module, the aforementioned signal processing and scoring module and model update function until the preset termination condition is met. The aforementioned hotspot output module is used to stop the search iteration operation when the above termination condition is met, and to output the target pose with the highest comprehensive physiological efficacy score among all historical observation data as the TMS motion hotspot; if there are multiple target poses with the same highest comprehensive physiological efficacy score, the one whose spatial position is closest to the above prior center coordinates is taken as the final TMS motion hotspot.

[0083] As an optional implementation, in the data acquisition and processing module, acquiring the subject's medical imaging data specifically involves acquiring the subject's T1 structural image, fMRI data, and DTI data, and processing the T1 structural image to reconstruct an individualized three-dimensional head model.

[0084] As an optional implementation, the data acquisition and processing module extracts prior center coordinates related to the target muscle motor function from the medical imaging data, including: a functional activation center coordinate extraction unit, used to locate brain functional areas activated by the task based on the fMRI data using statistical analysis methods to obtain functional activation center coordinates; and / or, a structural connectivity center coordinate extraction unit, used to locate cortical regions with the strongest anatomical connections to the corticospinal tract based on the DTI data using fiber tracing methods to obtain structural connectivity center coordinates; and a fusion unit, used to fuse or select the functional activation center coordinates and the structural connectivity center coordinates to determine the final prior center coordinates.

[0085] As an optional implementation, in the data acquisition and processing module, a search space is defined based on the prior center coordinates, and a candidate stimulus pose set is generated according to the search space. Specifically, this is used to: calculate the mapping point from the prior center coordinates to the individualized three-dimensional head model; define a spherical search region centered on the mapping point; discretize the surface where the spherical search region intersects with the scalp surface into a mesh to obtain multiple discrete scalp position points; and calculate the coil pose conforming to preset clinical standards based on the scalp surface normal of each discrete position point to generate the candidate stimulus pose set.

[0086] As an optional implementation, the signal processing and scoring module is specifically used to: extract amplitude features, latency features, and waveform similarity features from single-segment electromyographic signals acquired for each stimulus; calculate the average value and coefficient of variation of amplitude features for multiple stimulus signals of the same target pose, and convert the coefficient of variation into a signal stability score; standardize the average value of amplitude features and the average value of latency features; and perform weighted summation of the standardized amplitude value, standardized latency value, average waveform similarity, and signal stability score to obtain the comprehensive physiological efficacy score.

[0087] As an optional implementation, the model update unit in the model update and iteration control module is specifically used to efficiently update the parameters of the posterior distribution of the Gaussian process model through an incremental update algorithm.

[0088] As an optional implementation, the preset termination conditions include at least one of the following: the total number of iterations reaches a preset maximum value; the historical best comprehensive physiological efficacy score is not improved for a preset number of consecutive iterations in the local fine localization stage; the historical best comprehensive physiological efficacy score exceeds a preset effectiveness threshold.

[0089] In an exemplary embodiment, a computer device is provided, which may be a server or a terminal. The computer device includes a processor, memory, input / output (I / O) interfaces, and a communication interface. The processor, memory, and I / O interfaces are connected via a system bus, and the communication interface is connected to the system bus via the I / O interfaces. The processor of the computer device provides computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and internal memory. The non-volatile storage medium stores an operating system, computer programs, and a database. The internal memory provides an environment for the operation of the operating system and computer programs in the non-volatile storage medium. The database of the computer device stores video tag processing data. The I / O interfaces of the computer device are used for exchanging information between the processor and external devices. The communication interface of the computer device is used for communicating with external terminals via a network connection. When the computer program is executed by the processor, it implements a TMS motion hotspot automatic search method that integrates multimodal priors and Bayesian optimization.

[0090] In one exemplary embodiment, a computer device is also provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the above-described method embodiments.

[0091] In one exemplary embodiment, a computer-readable storage medium is provided storing a computer program that, when executed by a processor, implements the steps in the above-described method embodiments.

[0092] In one exemplary embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps in the above-described method embodiments.

[0093] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, data stored, data displayed, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties, and the collection, use and processing of the relevant data must comply with relevant regulations.

[0094] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium. When executed, the computer program can include the processes of the embodiments described above. Any references to memory, databases, or other media used in the embodiments provided in this application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM can take many forms, such as Static Random Access Memory (SRAM) or Dynamic Random Access Memory (DRAM).

[0095] The databases involved in the embodiments provided in this application may include at least one type of relational database and non-relational database. Non-relational databases may include, but are not limited to, blockchain-based distributed databases. The processors involved in the embodiments provided in this application may be general-purpose processors, central processing units, graphics processing units, digital signal processors, programmable logic devices, quantum computing-based data processing logic devices, etc., and are not limited to these.

[0096] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0097] This document uses specific examples to illustrate the principles and implementation methods of this application. The descriptions of the above embodiments are only for the purpose of helping to understand the methods and core ideas of this application. Furthermore, those skilled in the art will recognize that, based on the ideas of this application, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of this application.

Claims

1. A method for automatic search of motion hotspots in TMS that integrates multimodal priors and Bayesian optimization, characterized in that, Includes the following steps: S1: Acquire the subject's medical imaging data, extract the prior center coordinates related to the target muscle's motor function from the medical imaging data, define the search space based on the prior center coordinates, and generate a set of candidate stimulus poses based on the search space; The search iteration operation is performed repeatedly until the preset termination condition is met; Each search iteration includes: S2: Based on the current Gaussian process model and the set of candidate stimulus poses, the Bayesian optimization algorithm is used to determine the target pose of the next stimulus. Wherein, if the current search iteration operation is the first search iteration operation, the mean function of the current Gaussian process model is obtained by initialization based on the prior center coordinates; if the current search iteration operation is not the first search iteration operation, the current Gaussian process model is the Gaussian process model that updated the posterior distribution in the previous search iteration operation. S3: Control the actuator to move to the target pose, apply TMS, and simultaneously collect electromyographic signals of the target muscle; S4: Extract the multi-dimensional physiological features of the electromyographic signal, and calculate the comprehensive physiological efficacy score based on the multi-dimensional physiological features; use the target pose and the corresponding comprehensive physiological efficacy score as new observation data. S5: Use the target pose and its corresponding comprehensive physiological efficacy score as observation data to update the posterior distribution of the current Gaussian process model; When the preset termination condition is met, the search iteration operation is stopped, and the target pose with the highest comprehensive physiological efficacy score among all historical observation data is output as the TMS motion hotspot. If there are multiple target poses with the same highest comprehensive physiological efficacy score, the one whose spatial position is closest to the prior center coordinates is taken as the final TMS motion hotspot.

2. The automatic TMS motion hotspot search method integrating multimodal prior and Bayesian optimization as described in claim 1, characterized in that, Step S1 involves acquiring the subject's medical imaging data, extracting prior center coordinates related to the target muscle's motor function from the medical imaging data, defining a search space based on the prior center coordinates, and generating a set of candidate stimulus poses based on the search space; specifically including: S11, acquire the subject's T1 structural image, fMRI and DTI data, and perform quality checks; preprocess and segment the T1 structural image, and reconstruct an individualized three-dimensional head model; S12, Based on the fMRI data and / or DTI data, extract the prior center coordinates; S13, Calculate the shortest Euclidean distance from the prior center coordinates to the individualized 3D head model. Determine the coordinates of the scalp mapping point; construct a tangent plane centered on the coordinates of the scalp mapping point; generate two-dimensional grid points on the tangent plane based on a preset range; project the two-dimensional grid points onto the surface of the individualized three-dimensional head model to obtain multiple discrete scalp surface position coordinates; for each discrete scalp surface position coordinate, construct the corresponding stimulus pose based on its surface normal vector, and all stimulus poses constitute the candidate stimulus pose set.

3. The automatic TMS motion hotspot search method integrating multimodal prior and Bayesian optimization as described in claim 2, characterized in that, In step S12, extracting the prior center coordinates based on the fMRI data specifically includes: The fMRI data were preprocessed and registered with the T1 structural images; The activation statistics of each voxel in the fMRI data are calculated based on a general linear model; the hand representative area of ​​the primary motor cortex on the affected side is delineated as the region of interest on the individual T1 structural image; multiple comparison corrections are performed on the statistics within the region of interest, and the local peak points of the statistics are located in the significant activation clusters, and their coordinates are defined as the coordinates of the functional activation center.

4. The automatic TMS motion hotspot search method integrating multimodal prior and Bayesian optimization as described in claim 2, characterized in that, In step S12, extracting the prior center coordinates based on the DTI data specifically includes: The DTI data is preprocessed and registered with the T1 structural image; The anisotropy fraction FA of the DTI data is calculated, and a fiber tracing algorithm is used to trace fibers backward from the predetermined deep structures of the corticospinal tract on the affected side to the cortical surface. The traced fiber ends are projected onto the cortical surface, and the density of all fiber ends or their projection areas on the cortical surface are calculated. The point of highest density or the geometric center of the projection area is defined as the coordinate of the structural connection center.

5. The automatic TMS motion hotspot search method integrating multimodal prior and Bayesian optimization according to claim 2, characterized in that, In step S12, the extraction of prior center coordinates specifically involves: When the functional activation center coordinates of the fMRI data and the structural connectivity center coordinates of the DTI data are obtained simultaneously, if the spatial distance between the two is less than a threshold, the midpoint or the one with the higher signal-to-noise ratio is taken as the final prior center coordinates; otherwise, one of them is selected as the final prior center coordinates based on the clinical hypothesis that emphasizes functional recovery or structural connectivity.

6. The automatic TMS motion hotspot search method integrating multimodal prior and Bayesian optimization according to claim 1, characterized in that, The mean function mentioned in step S2 is specifically as follows: The mean function Using the aforementioned prior center coordinates The core is a Gaussian bias function: in, The preset confidence level range; Prior center coordinates; The radius of the prior confidence range; x is the spatial coordinate; It is a natural exponential function.

7. The automatic TMS motion hotspot search method integrating multimodal prior and Bayesian optimization as described in claim 1 or 6, characterized in that, In step S2, the determination of the next target pose to be stimulated using the Bayesian optimization algorithm also includes a stage switching strategy: The number of times a better score was not found in the monitoring. and the current best overall physiological efficacy score ; If satisfied Then switch to the exploration phase. Otherwise, remain in the current exploration phase; in The stagnation tolerance threshold, This is the physiological activation threshold.

8. The automatic TMS motion hotspot search method integrating multimodal prior and Bayesian optimization according to claim 1, characterized in that, In step S5, the step of using the target pose and its corresponding comprehensive physiological efficacy score as observation data to update the posterior distribution of the current Gaussian process model specifically involves: The newly obtained target pose and its comprehensive physiological efficacy score were added as data pairs to the historical observation dataset. Based on the updated historical observation dataset, the posterior distribution of the Gaussian process model is updated by Cholesky decomposition to obtain the updated prediction mean function and prediction variance function.

9. The automatic TMS motion hotspot search method integrating multimodal prior and Bayesian optimization according to claim 1, characterized in that, The preset termination condition includes at least one of the following: The total number of iterations has reached the preset maximum value; The best overall physiological efficacy score was not improved by the number of consecutive preset attempts during the local fine-tuning phase.

10. A TMS motion hotspot automatic search system integrating multimodal priors and Bayesian optimization, characterized in that, include: The system includes modules for data acquisition and processing, model building and decision making, stimulus execution and signal acquisition, signal processing and scoring, model updating and iteration control, and hotspot output. The data acquisition and processing module is used to acquire the subject's medical imaging data, extract the prior center coordinates related to the target muscle motor function from the medical imaging data, delineate the search space based on the prior center coordinates, and generate a set of candidate stimulus poses based on the search space. The model building and decision-making module is used to build a Gaussian process model and initialize the mean function of the Gaussian process model based on the prior center coordinates during the first search iteration operation; and in each search iteration operation, based on the current Gaussian process model and the candidate stimulus pose set, use the Bayesian optimization algorithm to determine the target pose of the next stimulus. The stimulation execution and signal acquisition module is used to control the actuator to move to the target pose, apply TMS, and simultaneously acquire the electromyographic signals of the target muscle. The signal processing and scoring module is used to extract multi-dimensional physiological features of the electromyographic signal and calculate a comprehensive physiological efficacy score based on the multi-dimensional physiological features. The model update and iteration control module is used to update the posterior distribution of the current Gaussian process model by using the target pose and its corresponding comprehensive physiological efficacy score as observation data; and to control the repeated execution of the search iteration operation consisting of the model construction and decision module, the stimulus execution and signal acquisition module, the signal processing and scoring module and the model update function until the preset termination condition is met. The hotspot output module is used to stop the search iteration operation when the termination condition is met, and to output the target pose with the highest comprehensive physiological efficacy score among all historical observation data as the TMS motion hotspot; if there are multiple target poses with the same highest comprehensive physiological efficacy score, the one whose spatial position is closest to the prior center coordinates is taken as the final TMS motion hotspot.

Images