Rehabilitation regulation method based on multi-modal brain network

Abstract

The application discloses a rehabilitation regulation method based on a multi-modal brain network, and relates to the technical field of neural network regulation and intelligent medical auxiliary technology. By collecting multi-modal physiological signals, K-SVD dictionary learning and variational modal decomposition are used to separate artifacts. Based on extended Kalman filtering, a neural activation index is estimated. Through conditional denoising diffusion probability model, an individualized whole brain functional connectivity matrix is generated, and precise target optimization is realized through a graph attention network. A physical information neural network is used to solve the Helmholtz equation electromagnetic positive problem and solve the optimal complex weight vector. Through triple causal decoupling, effective responses are purified. The regulation problem is modeled as a constrained Markov decision process, a Lyapunov safety function is introduced, and a safety closed-loop regulation is realized by using a constrained strategy optimization algorithm. A multi-modal large language model is called to generate a clinically interpretable rehabilitation report. The individualized neural rehabilitation regulation is precise, safe and interpretable.

Description

Technical Field

[0001] This invention relates to the field of neural network regulation and intelligent medical assistance technology, specifically to a rehabilitation regulation method based on multimodal brain networks. Background Technology

[0002] Existing rehabilitation systems that use EEG signal feedback primarily utilize EEG feedback sensors to simultaneously collect mixed signals from target neural networks, non-target brain regions, and environmental electromagnetic interference. Then, through specific neural networks and target points, specific frequencies, wavelengths, and repetitive magnetic waves are used to influence the neurotransmitter levels in the brain, regulate the release and absorption of neurotransmitters, and control different states of the user.

[0003] However, the amplitude of the neural response in the target area is inherently weak (µV level), generally below the noise floor. Furthermore, the quality of a single EEG signal fluctuates significantly due to scalp impedance, motion artifacts, and electromagnetic crosstalk from the stimulus itself, leading to feedback distortion in the control terminal and making it impossible to accurately determine the true response state of the target area. Simultaneously, the beamwidth of existing application units (coils / electrodes / transducers) is too large, and secondary diffusion after transmission caused by skull scattering allows physical field energy to enter non-target brain regions. Excitation interference in these non-target brain regions weakens the specificity of the control. Under safety constraints (SAR / thermal effects), this also limits the space for inputting higher effective energy into the target area, resulting in poor control effects.

[0004] Meanwhile, existing systems may misinterpret stimulus artifacts on EEG or respiratory sinus arrhythmia in HRV as patient sensitivity to that intensity. They may assume that all amplitude changes are valid "responses" of the brain to the stimulus, generating incorrect feedback that drives the controller to increase the stimulus intensity or change the phase, causing abnormal discharges or severe palpitations in non-target areas. In other words, the lack of dynamic perception of real-time signal-to-noise ratio (SNR) will affect the control accuracy of the control terminal. Summary of the Invention

[0005] The purpose of this invention is to provide a rehabilitation regulation method based on multimodal brain networks to solve the problems mentioned in the background art.

[0006] The specific technical solution provided by this invention is as follows: a rehabilitation regulation method based on multimodal brain networks, comprising the following operational steps: Multimodal physiological signals were acquired, and structured artifact separation was performed on the multimodal physiological signals to obtain multi-channel clean EEG signals, de-breathing HRV feature vectors, and fNIRS concentration signals. The real-time signal-to-noise ratio of each acquisition channel is calculated based on the clean EEG signals from multiple channels. The neural drive intensity is used as an augmented state variable based on the extended Kalman filter to fuse the multimodal physiological signals and estimate the neural activation index and its posterior estimation error variance. A conditional denoising diffusion probability model is introduced to generate an individualized whole-brain functional connectivity matrix. The individualized whole-brain functional connectivity matrix is ​​then fused with the structural connectivity matrix into a multimodal adjacency matrix. This matrix is ​​input into a graph attention network, and the regulatory influence score of each brain region is calculated. The brain region cluster with the highest score is selected as the individualized stimulation target and the effective connectivity pathway is output. A heterogeneous electromagnetic model is constructed based on individual magnetic resonance imaging. The electromagnetic forward problem of the heterogeneous electromagnetic model is solved using a physical information neural network to predict the electric field distribution of the whole head. The optimal complex weight vector is solved by combining time-reversal focusing and SAR-constrained convex optimization. A triple identification model is constructed, which includes a structural causal model, a convergent cross-mapping model, and a time series prediction model, to suppress false responses and extract and output valid responses. Finally, the regulation problem is modeled as a constrained Markov decision process, and a candidate Lyapunov function is introduced to use the effective response as feedback input for closed-loop regulation. The training graph convolutional autoencoder maps the sequence of functional connectivity matrices to a low-dimensional latent space to generate a rehabilitation progress trajectory. It then calls a multimodal large language model and combines interpretable feature contribution values ​​to generate a clinically interpretable rehabilitation report.

[0007] Furthermore, the multimodal physiological signals include electroencephalogram (EEG) signals, functional near-infrared spectroscopy (FIR) signals, electrocardiogram (ECG) signals, respiratory signals, and acceleration signals. Structured artifact separation of the multimodal physiological signals involves: training a dictionary learning algorithm on artifact waveform segments extracted during and after a preset time period following stimulation to construct a complete artifact dictionary; performing variational mode decomposition on the time series of the original EEG signal to obtain multiple intrinsic mode functions (EMFs); calculating the coherence between the EMFs and the stimulation rhythm; and selecting EMFs containing artifact components. Modal functions; solve the sparse coding problem of intrinsic modal functions containing artifact components on a complete artifact dictionary, reconstruct and subtract artifact components to obtain EEG signals without stimulation artifacts; using the acceleration signal as a reference, use a recursive least squares adaptive filter to perform online filtering on the EEG signals without stimulation artifacts to remove motion artifacts; using the respiratory signal and its orthogonal components after Hilbert transform as a reference, use a recursive least squares adaptive filter to perform adaptive noise cancellation on the HRV time series, and output the HRV index without respiratory modulation.

[0008] Furthermore, the calculation of the real-time signal-to-noise ratio (SNR) for each acquisition channel includes: within the time-locked window triggered by the stimulus, a preset duration after the stimulus is taken as the signal window, and a preset duration before the stimulus is taken as the baseline noise window; the EEG signals within the signal windows under multiple identical stimulus conditions are averaged using time-locked methods, and the energy within the time-locked average signal window is calculated as the signal energy; the mean of the signal variance within the baseline noise window of each trial is calculated as the noise energy; the single-channel SNR is calculated based on the signal energy and noise energy; when the single-channel SNR is lower than a set threshold, the single channel is marked as a low-confidence channel; and Bayesian change point detection is used simultaneously to monitor SNR mutations in real time.

[0009] Furthermore, using neural drive intensity as an augmented state variable based on extended Kalman filtering includes: defining an augmented state vector containing oxyhemoglobin concentration, the rate of change of oxyhemoglobin concentration, and neural drive intensity, wherein the neural drive intensity satisfies the random walk hypothesis, and the initial augmented state vector is set by resting-state data; using EEG signals... The band average power is used as the observation variable to establish the state evolution equation and the observation equation. An extended Kalman filter recursion is performed, including a prediction step and an update step. The neural drive intensity component in the posterior state estimate is output as the neural activation index, and the variance of the posterior estimation error is used as a confidence quantification of the neural activation index.

[0010] Furthermore, the introduction of a conditional denoising diffusion probability model to generate an individualized whole-brain functional connectivity matrix includes: performing fiber tract tracking on individual diffusion tensor imaging to obtain a normalized structural connectivity matrix; calculating the resting-state average functional connectivity matrix and extracting typical frequency band power spectra as spectral feature vectors; flattening the structural connectivity matrix and concatenating it with the spectral feature vectors, then compressing it into a low-dimensional conditional embedding via an encoder, which serves as the generation condition for the conditional denoising diffusion probability model; training the denoising network with the low-dimensional conditional embedding and the number of diffusion steps as input to predict the noise added in each diffusion step; after training, iteratively executing the denoising process starting from Gaussian noise to generate an individualized whole-brain functional connectivity matrix.

[0011] Furthermore, the calculation of the regulatory influence score for each brain region includes: using the sum of the product of the attention coefficient and the change in the weighted phase lag index between nodes after the pre-application of stimulus on the target network node as the regulatory influence score for each brain region.

[0012] Furthermore, the electromagnetic forward problem of the heterogeneous electromagnetic model is solved using a physical information neural network, which includes: segmenting six tissues—skull, cerebrospinal fluid, gray matter, white matter, skin, and air boundary—from T1-weighted images of individual magnetic resonance imaging; tracing the direction of white matter fiber bundles using diffusion tensor imaging, and assigning frequency-related conductivity and relative permittivity to each voxel; constructing a physical information neural network with spatial coordinates as input and electric field vector as output, the loss function including internal physical residual terms and boundary condition terms; and training the physical information neural network to obtain a whole-head electric field distribution prediction model.

[0013] Furthermore, the combination of time-reversal focusing and SAR-constrained convex optimization includes: placing a virtual dipole source at the individualized stimulus target location; calculating the impulse response generated by the virtual dipole source at each unit of the array using the physical information neural network; performing time reversal on the impulse response; convolving the time-reversed signal with complex weighted coefficients to obtain the excitation signal for each unit; and solving for the optimal complex weight vector with the objective function of maximizing the ratio of the quadratic form of the target region focusing matrix to the quadratic form of the specific absorption rate penalty matrix, and with the constraint that the quadratic form of the specific absorption rate penalty matrix does not exceed the maximum specific absorption rate limit.

[0014] Furthermore, a triple identification model is constructed, comprising a structural causal model, a convergent cross-mapping model, and a time-series prediction model. This includes: constructing a structural causal model to decouple multi-channel clean EEG signals and de-breathing HRV feature vectors from multimodal causality; using a convergent cross-mapping model to remove the unidirectional driving force of respiration on heart rate variability; and using a time-series prediction model to predict the EEG signal of the current stimulation segment and identify unmodeled artifacts based on the prediction residuals.

[0015] Furthermore, the candidate Lyapunov function includes a cortical excitability risk term and a specific absorptivity risk term. The sum of the squares of the amount by which the cortical excitability risk index exceeds the safety threshold and the amount by which the specific absorptivity exceeds the safety limit is used as the candidate Lyapunov function value. The constrained policy optimization algorithm is then used to update the stochastic policy, maximizing the reward within the trust domain while forcibly satisfying the safety constraints.

[0016] Compared with the prior art, the beneficial effects achieved by the present invention are: (1) This invention constructs a complete artifact dictionary and solves sparse coding through a structured artifact separation method, removes motion artifacts and respiratory sinus arrhythmia, and classifies transcranial magnetic stimulation discharge artifacts, motion artifacts and respiratory sinus arrhythmia into different dictionary atoms to achieve structured separation, thus overcoming the defect of existing single-mode denoising methods that cannot distinguish mixed interference with highly overlapping time-frequency features.

[0017] (2) This invention constructs a physical information neural network and embeds the dielectric properties of individualized tissues into the constraints of the Helmholtz equation, thereby realizing the real-time solution of individualized electromagnetic positive problems. By combining time-reversal focusing and specific absorptivity constraint convex optimization, the specific absorptivity in the non-target area is constrained within a safe limit while maximizing the focused energy in the target area. This solves the problem that existing methods lack the means to constrain the specific absorptivity at the level of solving electromagnetic positive problems and have the risk of energy leakage in the non-target area from the perspective of solving electromagnetic fields.

[0018] (3) This invention solves the problem that the existing closed-loop system is prone to being misjudged as a valid neural response by the controller due to residual stimulus artifacts and respiratory confusion signals by triple identification markers, thus ensuring that all signals entering the closed-loop controller are real brain state responses that have been causally purified. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of the steps of the rehabilitation and regulation method provided in the embodiment of the present invention; Figure 2 This is a schematic diagram of the multimodal signal synchronous acquisition and structured artifact separation structure provided in an embodiment of the present invention; Figure 3 This is a schematic diagram illustrating the principle of generating individualized functional connectivity matrices using the Conditional Denoising Diffusion Probability Model (DDPM) provided in this embodiment of the invention. Figure 4 This is a schematic diagram of the target optimization structure of the graph attention network (GAT) provided in an embodiment of the present invention. Detailed Implementation

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

[0021] Example 1: like Figure 1 As shown in the figure, the rehabilitation modulation method based on multimodal brain networks described in this embodiment includes the following steps: Step S1: Construct a multimodal signal synchronous acquisition system for EEG-fNIRS-ECG-respiration-accelerometer, and adaptive structured artifact separation.

[0022] In this embodiment, to address the problem that single-modal denoising in existing systems cannot distinguish between mixed interference from stimulus artifacts, motion artifacts, and respiratory sinus arrhythmia (RSA), this invention uses K-SVD dictionary learning in conjunction with variational mode decomposition (VMD) to classify various artifacts into different dictionary atoms, accurately removing the confusion between TMS artifacts, motion artifacts, and respiratory sinus arrhythmia, and achieving low-latency structured separation.

[0023] For example, specific implementations of multimodal signal synchronous acquisition and adaptive structured artifact separation include: Multimodal synchronous acquisition and preprocessing: In this embodiment, combined with Figure 2 As shown, this invention simultaneously records 32-channel EEG (covering bilateral frontal, central, parietal, and temporal lobes, including target areas such as the left DLPFC), 18-channel functional near-infrared spectroscopy (fNIRS, light source-detector spacing 3cm, covering the frontal and motor cortex), single-lead electrocardiogram (ECG), respiratory-induced plethysmography (RIP), and a triaxial accelerometer (placed in the frontal or temporal region) at a 4kHz sampling rate. All devices are synchronized via hardware TTL pulses. The EEG channel outputs scalp voltage (μV) signals, and the fNIRS outputs changes in oxygenated hemoglobin concentration. ECG provides the RR interval, and RIP outputs respiratory phase signals. Accelerometer output triaxial accelerometer signal This ensures that the reference signal is strictly aligned with the EEG sample during subsequent artifact removal.

[0024] Constructing a dictionary of stimulus artifacts: This invention extracts artifact waveform segments from all EEG channels within 1 ms after the firing of stimulators such as TMS / tACS (before the neuron generates an evoked response). (Length L is approximately 4 sampling points, i.e., 1ms × 4kHz). A large number of pseudo-trace segments are then used to form a training matrix. A complete dictionary of artifacts is learned using the K-SVD algorithm. ( For example, K=256, meaning the total number of dictionary atoms K in the complete pseudo-dictionary is 256, and each dictionary atom Describe a typical pseudotrace waveform basis, such that any pseudotrace segment can be sparsely represented as ,in A sparse coefficient vector (number of non-zero elements) ).

[0025] Variational Mode Decomposition (VMD) and Artifact Sparse Projection: Further analysis of the original EEG time series for each channel Through VMD adaptive decomposition into One intrinsic mode function (IMF) :

[0026] The center frequencies of each intrinsic mode function are arranged from low to high. Stimulus artifacts, due to their fixed pulse repetition frequency (e.g., 10Hz), tend to concentrate in a specific intrinsic mode function. By calculating the coherence between each intrinsic mode function and the stimulus rhythm, modes containing artifact components are automatically selected. The selected modes are then discretized into mode vectors containing stimulus artifacts. In a complete dictionary of forgeries Solve the sparse coding problem and calculate the optimal solution. :

[0027] in, The regularization parameter is used to reconstruct the artifact after obtaining the optimal solution, thus discretizing the mode into a vector. Subtracting the reconstructed artifacts The corrected intrinsic mode functions are then superimposed with other intrinsic mode functions to reconstruct the EEG signal with destimulation artifacts removed. .

[0028] Eliminating motion artifacts and respiratory sinus arrhythmia (RSA): using triaxial accelerometer signals For reference, a recursive least squares (RLS) adaptive filter was used to process the EEG signal after destimulation artifacts. Online filtering is performed, and the RLS filter weight vector is updated recursively to minimize the power of the filtered output, thereby eliminating artifacts related to head movement. Similarly, respiratory phase signals are used... Using the orthogonal components derived from the ECG and its Hilbert transform as reference input, another RLS filter is used to adaptively cancel noise in the ECG-derived HRV time series (such as instantaneous heart rate sequences). The filter output is the net autonomic HRV index after removing respiratory modulation, such as the RMSSD (root mean square of the difference between adjacent RR intervals) and the LF / HF ratio after RSA removal. The fNIRS concentration signal is also processed using the same reference signal to remove low-frequency oscillations caused by respiration. Finally, the output is a multi-channel clean EEG signal, a respiratory-free heart rate variability (HRV) eigenvector, and a functional near-infrared spectroscopy (fNIRS) concentration signal.

[0029] Step S2: Real-time sensing of signal-to-noise ratio and simultaneous extraction of multimodal spatiotemporal fusion features.

[0030] In this embodiment, the present invention utilizes SNR dynamic sensing to prevent low-quality segments from entering the control loop, and at the same time utilizes EKF to fuse EEG / fNIRS to generate robust neural drive estimation, providing high-quality input features with uncertainty for brain network target identification.

[0031] For example, the specific implementation process includes: Estimate the real-time signal-to-noise ratio (SNR) and label low confidence levels: Within the time-locked window triggered by the stimulus, define the time of each stimulus pulse as... Take it after For the "signal window", take the pre-stimulation value - This is referred to as the "baseline noise window". For a single channel, time-locked averaging of EEG signals within multiple signal windows under the same stimulus conditions is performed to obtain the event-related potential (ERP), and the energy within the time-locked average signal window is calculated. :

[0032] in, This represents the total number of sampling points within the signal window. The average ERP during the lockout period is the first time within the window. The amplitude of each sampling point, for each trial Calculate the variance of the signal within its baseline window. (The average of the squared variances after removing the mean) is used to calculate the baseline noise energy. :

[0033] Further calculations yielded the single-channel signal-to-noise ratio. :

[0034] If the single-channel signal-to-noise ratio If the signal energy does not significantly exceed the noise floor, the channel is marked as "low confidence" and its weight is reduced or it is directly excluded during subsequent fusion. At the same time, Bayesian change point detection is used to monitor the sudden changes in signal-to-noise ratio in real time to cope with changes in electrode contact state.

[0035] Multimodal neural drive fusion based on extended Kalman filter (EKF): First, a state-space model is established, fusing the neural-vascular coupling dynamics of fNIRS and EEG, and defining the state vector. for Using the latent variable "neural drive" to be estimated as the control input, the oxygenated hemoglobin concentration and its rate of change at time t is used to obtain the following state evolution equation:

[0036] in, Here is the state transition matrix. for The state vector at time -1 For the input matrix, for The implied neural drive strength to be estimated at time step. This is the process noise vector.

[0037] Based on the observation equation, EEG Band power Related to blood oxygenation status:

[0038] in, For a moment The observed variables, i.e. EEG at all times Average power in the band (30-80Hz), For the observation matrix, To observe noise (scalar), due to the implicit neural drive strength Since it cannot be directly measured, this invention further extends it to include time. The concentration of oxyhemoglobin and its rate of change and the state vector of the control input. In the middle, it is assumed that it follows a random walk. The new equation of state is:

[0039] in, for The augmented state vector at time step 1. for The augmented state vector at time -1 Random walk noise on neural drives The observation equation is further expressed as:

[0040] Real-time recursive estimation is performed using extended Kalman filtering: The prediction step is as follows:

[0041]

[0042] The update steps are as follows:

[0043]

[0044]

[0045] in, For prior state estimation, i.e. based on the first Information at time -1 is relevant to the... The predicted value of the state at time step. For the augmented state transition matrix, , The posterior state estimate of the previous time step, i.e., the th... The optimal state estimate at time -1 incorporates observation information, with the initial value... Set by resting state data. To estimate the error covariance matrix a priori, Let $\mathbf{a}$ be the posterior estimation error covariance matrix of the previous time step. Let be the process noise covariance matrix.

[0046] in, Here is the Kalman gain matrix. To augment the observation matrix, , To observe the noise variance, For posterior state estimation, For prior state estimation, Represents the observation residual term. For the posterior estimation of the error covariance matrix, It is the identity matrix. This represents a factor that reduces the uncertainty introduced by the observation.

[0047] This invention, through each time step : Prediction: Using a blood oxygen-neurodynamic model ( And the random walk assumption, based on the optimal state at the previous time step. Derive the prior guess at the current moment. At the same time, it conveys uncertainty. .

[0048] Update: When the new EEG is obtained Power observation Then, the difference between the model's predicted values ​​and the actual values ​​is calculated using Kalman gain. To determine the degree to which the model believes predictions or observations, we obtain the optimal posterior state. At the same time, tighten the covariance .

[0049] The final output is a neural activation index that integrates EEG rapid kinetics and fNIRS blood metabolic constraints. That is, the posterior state The third component is used, and the variance of the posterior estimation error is used as the confidence level of the neural activation index. The larger the variance, the less reliable the current estimate is.

[0050] Step S3: Introduce the Conditional Denoising Diffusion Probability Model (Conditional DDPM) to generate an individualized whole-brain high-resolution functional connectivity matrix. Select the optimal stimulus target and effective connectivity pathway in multimodal connectivity through a graph attention network.

[0051] In this embodiment, the present invention uses individual structural connectivity and resting-state spectrum as conditions, generates a high-resolution individualized functional connectivity matrix through a Conditional Denoising Diffusion Probability Model (Conditional DDPM), and then uses a Graph Attention Network (GAT) to select the optimal regulatory target and effective connectivity pathway in multimodal connectivity.

[0052] For example, the specific implementation steps include: Constructing an individual multimodal connectome: This invention uses diffusion tensor imaging (DTI) to trace fiber tracts in an individual, obtaining the number / probability of fibers in brain regions, and then normalizing them to form a structural connectivity matrix. N represents the number of brain regions. The resting-state average functional connectivity matrix is ​​then calculated based on the dynamic adjacency matrix. And extract the power spectrum of typical frequency bands as the spectral feature vector. Finally, define the condition vector, which is the structure connection matrix. Flattening and spectral eigenvectors The data is spliced ​​together and compressed by an encoder into a low-dimensional conditional embedding. This serves as a condition for generating the diffusion model.

[0053] Individualized Full-Time (FC) generated based on conditional denoising diffusion probability model: combined with Figure 3 As shown, the diffusion model is used to learn the inverse denoising process, from pure noise... Gradually generate individualized functional connectivity matrices The conditions are And the actual functional connection matrix. The process of gradually adding Gaussian noise to achieve forward diffusion is repeated while simultaneously training a conditional... A denoising network with step number t as input Used to predict added noise After training, iterative execution yields the individualized functional connectivity matrix. .

[0054] Target optimization for graph attention networks: combining Figure 4 As shown, the generated individualized functional connection matrix and structural connection matrix Merged into a multimodal adjacency matrix ,in As a weighted fusion coefficient, the input is a graph attention network (GAT). Node importance is assessed based on the efficiency of excitation / inhibition balance regulation in the target region, and the regulatory influence score of each brain region is calculated. :

[0055] in, The attention coefficient dynamically reflects the strength of effective connections and the priority of regulation. For nodes after pre-application of stimulation With downstream nodes The change in the weighted phase lag index wPLI between (which can be obtained from preliminary data or digital twin simulation in this embodiment) reflects the change from arrive Effective connectivity plasticity, and the score of regulatory influence of various brain regions. A higher score indicates greater efficiency in stimulating the target network in that region. The cluster with the highest score is designated as the individualized target point, and the output of this target point's... The strongest effective connection path revealed.

[0056] Step S4: Solve the electromagnetic positive of the individualized head model based on the physical information neural network, combined with time inversion (TR), holographic wavefront synthesis and SAR-constrained convex optimization.

[0057] In this embodiment, to address the issues of energy dissipation and non-target leakage, the present invention utilizes a Physical Information Neural Network (PINN) to solve the electromagnetic positive problem of the individualized head model, and combines time reversal (TR), holographic wavefront synthesis, and SAR constraint optimization.

[0058] For example, the specific implementation process includes: Constructing an individual heterogeneous electromagnetic model: T1-weighted and DTI images were extracted from individual magnetic resonance imaging (MRI). T1-weighted images were used to accurately segment head tissue types—distinguishing between skull, cerebrospinal fluid, gray matter, white matter, skin, and air boundaries. DTI images were used to trace the direction of white matter fiber bundles, constructing a structural connectivity matrix. Frequency-dependent dielectric properties (conductivity) were assigned to each voxel. Relative permittivity Establish the electromagnetic solution domain for the head. (Internal voxel set), the electric field intensity vector E in the brain satisfies the Helmholtz equation:

[0059]

[0060] in, The vector Laplace operator describes the curvature and diffusion of the electric field in space and is used to calculate the difference between the electric field at a point in space and the average electric field around it. Let be the square of the complex wave number, and define the propagation characteristics of electromagnetic waves in brain tissue, a lossy medium. Angular frequency, For vacuum permeability and dielectric constant, The imaginary unit, Constitutes the real part, It constitutes the imaginary part.

[0061] Conductivity for each individual and each brain region and relative permittivity These are different problems. This invention addresses them by spatially diffusing the electric field intensity vector E. ) and propagation attenuation ( ) are connected, and And through and Individualized brain region anatomy information (from T1 / DTI) and stimulation frequency They are merged into a single unified variable.

[0062] Introducing PINN to solve electromagnetic forward problems: Constructing a neural network Input spatial coordinates Predict and output the electric field vector at that point. By automatically differentiating the spatial derivative and embedding physical constraints, the loss function is defined as:

[0063] in, Let be the total loss function of the physical information neural network. This represents the number of internal sampling points. The Laplace operator represents the electric field. The number of boundary sampling points, For the boundary point The reference electric field value is given by measurement or simulation.

[0064] The first term (physical residual): forces the neural network to satisfy the Helmholtz equation within the solution domain.

[0065] The second condition (boundary condition) is to make the network approximate the tangential electric field given by the array elements through actual measurement or simulation at the scalp / air boundary. .

[0066] After training, PINN can instantly predict the full-head electric field distribution under arbitrary array excitation, replacing the time-consuming iteration of traditional mesh simulation.

[0067] Time-reversal focusing and SAR-constrained optimization: A virtual dipole source is placed at the target location, and the impulse response generated by the virtual dipole source at each element of the array is calculated using PINN. The response is obtained by time reversal. The excitation signal for each unit is ,in The weighting coefficients are complex, and wavefront shaping is achieved by adjusting these coefficients. To maximize the target region energy and minimize the non-target region specific absorptivity (SAR), the following constrained optimization problem is solved:

[0068]

[0069] in, Complex weight vectors are used to control the amplitude and phase of each unit. The conjugate transpose of the complex weight vector is used together as the optimization variable. For target area focusing matrix ( (a positive semi-definite matrix), the electric field vector at the target point predicted by PINN is backpropagated to the array end, obtaining the contribution vector of each array element to the field at that point, making the measurement of the focal area received power... It is proportional to the focal energy density. SAR penalty matrix ( (A positive semi-definite matrix), jointly constructed from the conductivity of each tissue and the predicted array response field, enabling effective SAR indices. Represents peak spatial SAR or global SAR; in the constraints, To define the maximum specific absorption rate, this invention maximizes the objective function to determine a set of weights. This maximizes the energy generated by each unit of SAR at the target point, achieving the most efficient and safest focusing.

[0070] Step S5: Decouple the multimodal causality of stimulus-response and suppress spurious responses.

[0071] In this embodiment, the problem of spurious responses such as residual stimulus artifacts and respiratory sinus arrhythmia (RSA) being misidentified as valid responses by the controller is addressed. This invention employs a triple identification method: constructing a structural causal model (SCM), converging cross-mapping (CCM), and using a large temporal prediction model (PatchTST) to predict residuals. This decouples genuine neural responses, residual artifacts, and respiratory confusion.

[0072] For example, the specific implementation process includes: Constructing a structural causal model (SCM): Defining key variables and their causal relationships: S: Stimulation parameters (intensity, frequency, etc.); N: Real neural response (target area excitation / inhibition); A: Residual irritant artifacts; R: Respiratory phase; H: HRV indicator (such as RMSSD); Causality edge: (Stimulus produces artifacts). (HRV is modulated by both respiratory and neural pathways), and the SCM diagram explicitly encodes obfuscated paths to provide prior information for subsequent decoupling.

[0073] Converging Cross-Mapping (CCM) Causal Stripping RSA: The Converging Cross-Mapping (CCM) algorithm is used to delay embedding and reconstruct attractors to detect causal drivers between variables. For respiratory signals and HRV time series, the convergent cross-mapping is calculated. If the calculated CCM correlation coefficient... and If significantly lower, it confirms unidirectional respiratory drive HRV. Then, using counterfactual reasoning, the respiratory phase is fixed at the end of expiration, and the expected HRV under this intervention is calculated using a regression model as a net autonomic indicator for "de-RSA" to eliminate the spurious "response" caused by respiratory fluctuations.

[0074] Residual identification of the large-scale temporal prediction model: The PatchTST large-scale temporal prediction model, pre-trained on general neural signals, is used to predict the EEG of the current stimulus segment. The model takes the historical window and current stimulus parameters as input, and outputs a one-step prediction with RSA-degraded features as covariates. If the fluctuation range of the prediction residual exceeds a set multiple of the baseline noise standard deviation, the fluctuation is marked as an unmodeled artifact or a precursor to abnormal discharge, and is removed from the effective response and not included in the closed-loop control.

[0075] Extract and output valid responses: Observed power changes or potential amplitudes are marked as valid responses only if all three of the following conditions are met simultaneously. : Condition 1: EEG after removing artifacts and removing causal factors related to breathing There is a significant time lock between band power and stimulus; Condition 2: The neural response → blood flow response in the CCM is significant in one direction, which meets the necessary causal condition; Condition 3: If the predicted residual does not exceed the threshold, abnormal artifacts are excluded.

[0076] Therefore, the closed-loop controller only receives the real brain state responses that have been purified of causality and verified for reliability.

[0077] Step S6: Closed-loop adaptive control based on constrained Markov decision and Lyapunov safety.

[0078] In this embodiment, under the premise of ensuring absolute safety (no abnormal discharge, no exceeding SAR, no severe discomfort), the present invention adaptively adjusts the stimulation parameters (intensity, frequency, pulse width, phase, etc.) to make the neural response in the target area tend towards the optimal state. The modulation problem is modeled as a reinforcement learning problem with safety constraints, and Lyapunov stability theory is used to ensure that the safe state converges exponentially.

[0079] For example, specific implementations include: Constrained Markov Decision Process (CMDP) Modeling: Define state At time t, the state consists of multimodal feature vectors, including neural drive estimation, current signal-to-noise ratio (mean or minimum of each channel), target region functional connectivity strength (e.g., mean wPLI with the target network), cumulative SAR value (average SAR within the past time window), cortical excitability index, etc. Define Action : Stimulus parameter vector, which may include magnetic field strength (% of maximum output), frequency (Hz), pulse width ( The phase (in radians) is output by the policy network after being discretized or made continuous. The action space is limited by the physical constraints of the device.

[0080] Define rewards Defined as the degree of matching between the neural response intensity of the target area and the target response curve, such as the expected... The power increase is reduced by the penalty for power leakage in the non-target area.

[0081] Security Cost The safety cost of each step consists of two parts: cortical excitability risk and SAR risk, namely the degree to which cortical excitability risk exceeds the set safety threshold and the degree to which instantaneous SAR exceeds the limit.

[0082] Constructing Lyapunov security constraints: To enforce security, this invention introduces a candidate Lyapunov function. The "energy" that measures the distance between the current state and a dangerous state:

[0083] in, It uses the ReLU function, which only calculates the cost when a threshold is exceeded. The current risk index for cortical excitability is obtained by detecting a sudden drop in deoxyhemoglobin using the EEG high-frequency oscillation (>80Hz) envelope and fNIRS. These are set clinical / physiological safety thresholds; exceeding them indicates a risk of abnormal discharge. The peak space at the current moment is predicted in real time by PINN. To set safety standard limits, the safety constraint requires that after each policy update, the expected Lyapunov function value decays exponentially. This constraint ensures that even if the state temporarily deteriorates, it will quickly return to the safe domain.

[0084] Constrained Policy Optimization (CPO) and Large Model Pre-training: A Constrained Policy Optimization algorithm is used to update the stochastic policy. CPO maximizes rewards within the trust domain while enforcing safety constraints, obtaining safe and monotonically improving policy parameters through second-order approximation and Lagrange duality. To avoid unsafe exploration in the early stages of online learning, a pre-trained Transformer model is used to predict a safe and effective range of initial stimulus parameters (e.g., an intensity upper limit of 90% of the resting motion threshold) based on individual resting-state multimodal characteristics (such as power spectrum and connectivity matrix). This model is trained on a large amount of offline simulation data, enabling online CPO to start learning from a safe starting point and reducing exploration risks.

[0085] Step S7: Multi-scale neural mechanism co-evaluation and interpretability report generation In this embodiment, changes in cortical plasticity and neurotransmitter levels are quantified through spatiotemporal decomposition of multimodal networks, and an interpretable rehabilitation report that conforms to clinical thinking is automatically generated by calling a multimodal large language model (MLLM).

[0086] For example, specific implementations include: Multiscale mechanism quantification: Extracting excitation / inhibition balance (non-periodic component of EEG1 / f slope), cortical plasticity (long-term changes in N100 / P300), indirect neurotransmitter indicators (decoupling of fNIRS suggests dopamine / norepinephrine regulation), cerebrovascular reactivity (breath-holding index correlated with CO2), and HRV mitochondrial efficiency index (ratio of total HRV power to oxygen consumption) from fusion features.

[0087] Brain network reorganization assessment: Calculate the graph theory differences (clustering coefficients, global efficiency) of individualized connection matrices before and after modulation, then train a graph convolutional autoencoder to encode the high-dimensional connection matrix into a two-dimensional latent space. Input the connection matrix sequence during the modulation process into the encoder to obtain a rehabilitation progress trajectory, visualizing the brain network's shift towards a healthy state.

[0088] Interpretable generation of multimodal large language models (MLLM): Input serialization formats the following information into a text prompt: Patient's basic information and changes in anxiety scale scores (such as GAD-7); Control parameter sequence (intensity and frequency each time); Time series of quantitative mechanism indicators (E / I index, cerebrovascular reactivity, total HRV power, etc.); Brain network graph theory index changes and trajectory description; Significant functional connectivity changes (e.g., a 40% decrease in DLPFC-sgACC connectivity).

[0089] MLLM processing: Using a finely tuned visual-language large model (such as LLaVA-Med or GPT-4V), the above data tables and graphs are used as input. The model locates key changes through a causal attention mechanism and generates a natural language report.

[0090] SHAP Interpretability: To explain which features contribute most to the prediction of mood improvement, the model output is analyzed using SHAP values. For example, the Shapley value for the probability of anxiety relief is calculated for each input feature, the top three features contributing the most (such as "left DLPFC1 / f slope decreased by 0.15", "HRV total power increased by 30%", etc.), and the weight explanation is given in the report.

[0091] By breaking down the above mechanisms and combining them with the natural language interpretation generated by MLLM, clinicians or patients can clearly understand how each modulation affects the "neurotransmitter-network-emotion" pathway and optimize subsequent plans based on the reports.

[0092] It should be noted that, in this invention, relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus.

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

Claims

1. A rehabilitation modulation method based on multimodal brain networks, characterized in that, The following steps are included: Multimodal physiological signals were acquired, and structured artifact separation was performed on the multimodal physiological signals to obtain multi-channel clean EEG signals, de-breathing HRV feature vectors, and fNIRS concentration signals. The real-time signal-to-noise ratio of each acquisition channel is calculated based on the clean EEG signals from multiple channels. The neural drive intensity is used as an augmented state variable based on the extended Kalman filter to fuse the multimodal physiological signals and estimate the neural activation index and its posterior estimation error variance. A conditional denoising diffusion probability model is introduced to generate an individualized whole-brain functional connectivity matrix. The individualized whole-brain functional connectivity matrix is ​​then fused with the structural connectivity matrix into a multimodal adjacency matrix. This matrix is ​​input into a graph attention network, and the regulatory influence score of each brain region is calculated. The brain region cluster with the highest score is selected as the individualized stimulation target and the effective connectivity pathway is output. A heterogeneous electromagnetic model is constructed based on individual magnetic resonance imaging. The electromagnetic forward problem of the heterogeneous electromagnetic model is solved using a physical information neural network to predict the electric field distribution of the whole head. The optimal complex weight vector is solved by combining time-reversal focusing and SAR-constrained convex optimization. A triple identification model is constructed, which includes a structural causal model, a convergent cross-mapping model, and a time series prediction model, to suppress false responses and extract and output valid responses. Finally, the regulation problem is modeled as a constrained Markov decision process, and a candidate Lyapunov function is introduced to use the effective response as feedback input for closed-loop regulation. The training graph convolutional autoencoder maps the sequence of functional connectivity matrices to a low-dimensional latent space to generate a rehabilitation progress trajectory. It then calls a multimodal large language model and combines interpretable feature contribution values ​​to generate a clinically interpretable rehabilitation report.

2. The rehabilitation regulation method based on multimodal brain networks according to claim 1, characterized in that, The multimodal physiological signals include: electroencephalogram (EEG) signals, functional near-infrared spectroscopy (FIR) signals, electrocardiogram (ECG) signals, respiratory signals, and acceleration signals; The structured artifact separation of multimodal physiological signals includes: A dictionary learning algorithm is used to train the artifact waveform segments captured during the stimulation discharge and within a preset time after the stimulation discharge begins, thereby constructing a complete artifact dictionary. Variational mode decomposition was performed on the time series of the original EEG signal to obtain multiple intrinsic mode functions. The coherence between the intrinsic mode functions and the stimulus rhythm was calculated, and the intrinsic mode functions containing artifact components were selected. The sparse coding problem is solved on the eigenmode function containing artifact components on a complete artifact dictionary. The artifact components are then reconstructed and subtracted to obtain the EEG signal with destimulated artifacts. Using acceleration signals as a reference, a recursive least squares adaptive filter is used to filter EEG signals for destimulation artifacts online to remove motion artifacts. Using the respiratory signal and its orthogonal components after Hilbert transform as references, a recursive least squares adaptive filter is used to adaptively eliminate noise in the HRV time series, and the output HRV index is de-modulated.

3. The rehabilitation regulation method based on multimodal brain networks according to claim 2, characterized in that, Calculating the real-time signal-to-noise ratio of each acquisition channel includes: Within the time-locked window triggered by the stimulus, the preset duration after the stimulus is taken is the signal window, and the preset duration before the stimulus is taken is the baseline noise window. Time-locked averaging was performed on EEG signals within signal windows under multiple identical stimulus conditions, and the energy within the time-locked average signal window was calculated as the signal energy. The mean of the signal variance within the baseline noise window for each trial is calculated as the noise energy. Calculate the signal-to-noise ratio of a single channel based on signal energy and noise energy; When the signal-to-noise ratio of a single channel is lower than a set threshold, the single channel is marked as a low-confidence channel. Simultaneously utilize Bayesian change point detection to monitor signal-to-noise ratio mutations in real time.

4. The rehabilitation regulation method based on multimodal brain networks according to claim 3, characterized in that, Based on the extended Kalman filter, neural drive strength is used as an augmented state variable, including: Define an augmented state vector that includes oxyhemoglobin concentration, oxyhemoglobin concentration change rate, and neural drive strength, wherein the neural drive strength satisfies the random walk hypothesis, and the initial augmented state vector is set by resting state data; With EEG signals The average power of the band is used as the observation variable to establish the state evolution equation and the observation equation; The extended Kalman filter recursion is performed, including a prediction step and an update step. The neural drive strength component in the posterior state estimate is output as the neural activation index, and the variance of the posterior estimation error is used as a confidence quantification of the neural activation index.

5. The rehabilitation regulation method based on multimodal brain networks according to claim 4, characterized in that, The introduction of a conditional denoising diffusion probability model to generate a personalized whole-brain functional connectivity matrix includes: Fiber tract tracing was performed on individual diffusion tensor imaging to obtain the normalized structural connectivity matrix; Calculate the resting-state average functional connectivity matrix and extract the power spectrum of typical frequency bands as spectral feature vectors; The structural connection matrix is ​​flattened and concatenated with the spectral feature vector, then compressed by the encoder into a low-dimensional conditional embedding, which serves as the generation condition for the conditional denoising diffusion probability model. The denoising network is trained with low-dimensional conditional embedding and diffusion steps as inputs to predict the noise added at each diffusion step. After training, the denoising process is iteratively executed starting from Gaussian noise to generate an individualized whole-brain functional connectivity matrix.

6. The rehabilitation regulation method based on multimodal brain networks according to claim 5, characterized in that, The calculation of the regulatory influence score for each brain region includes: using the sum of the product of the attention coefficient and the change in the weighted phase lag index between nodes after the pre-application of stimulation on the target network node as the regulatory influence score for each brain region.

7. The rehabilitation regulation method based on multimodal brain networks according to claim 6, characterized in that, Solving electromagnetic forward problems of heterogeneous electromagnetic models using physical information neural networks includes: Segmenting six tissue types—skull, cerebrospinal fluid, gray matter, white matter, skin, and air boundary—from individual T1-weighted magnetic resonance imaging images; Diffusion tensor imaging was used to track the direction of white matter fiber bundles, and each voxel was assigned frequency-dependent conductivity and relative permittivity. A physical information neural network is constructed, with spatial coordinates as input and electric field vector as output. The loss function includes internal physical residual terms and boundary condition terms. Train a physical information neural network to obtain a prediction model of the electric field distribution across the entire head.

8. The rehabilitation regulation method based on multimodal brain networks according to claim 7, characterized in that, Combining time-reversal focusing and SAR-constrained convex optimization includes: A virtual dipole source is placed at a personalized stimulation target location, and the impulse response generated by the virtual dipole source at each unit of the array is calculated by the physical information neural network. The impulse response is time-reversed, and the time-reversed signal is convolved with the complex weighted coefficients to obtain the excitation signal of each unit. The objective function is to maximize the ratio of the quadratic form of the target region focusing matrix to the quadratic form of the specific absorption rate penalty matrix, with the constraint that the quadratic form of the specific absorption rate penalty matrix does not exceed the maximum specific absorption rate limit. The optimal complex weight vector is then solved.

9. The rehabilitation regulation method based on multimodal brain networks according to claim 8, characterized in that, The construction of a triple identification model, comprising a structural causal model, a convergent cross-mapping model, and a time series prediction model, includes: A structural causal model was constructed to decouple multimodal causality between multichannel clean EEG signals and debreathing HRV feature vectors; A convergent cross-mapping model was used to isolate the unidirectional driving force of respiration on heart rate variability; A time-series prediction model is used to predict the EEG signal of the current stimulus segment, and unmodeled artifacts are identified based on the prediction residuals.

10. The rehabilitation regulation method based on multimodal brain networks according to claim 9, characterized in that, The candidate Lyapunov function includes: cortical excitability risk term and specific absorption rate risk term, and the sum of the squares of the amount by which the cortical excitability risk index exceeds the safety threshold and the amount by which the specific absorption rate exceeds the safety limit is used as the candidate Lyapunov function value; Further, a constraint policy optimization algorithm is used to update the random policy, which maximizes the reward within the trust domain while forcibly satisfying the safety constraints.

Images