A brain-computer interaction neural patient rehabilitation training method based on digital twinning

Abstract

The application discloses a brain-computer interaction neural patient rehabilitation training method based on digital twinning, and belongs to the technical field of intelligent rehabilitation and brain-computer interface, and comprises the following steps: taking real-time brain electrical signals of a patient as driving, performing online state-parameter joint estimation on a multi-node Jansen-Rit neural mass model through an unscented Kalman filtering algorithm, constructing an active digital twin, tracking a brain network state in real time, and quantifying and predicting a neural plasticity window period according to a neural parameter evolution sequence; driving the neural mass model simulation by applying virtual parameter disturbance to the active digital twin, generating individualized labeled synthetic electroencephalogram data, performing zero-sample virtual pre-training on a decoding network, feeding back a measured response to the active digital twin for correction, dynamically adjusting training intensity and training frequency according to the neural plasticity window period, and promoting neural pathway reconstruction of the patient; and the application realizes dynamic feedback of neural state individualized perception and training conditions, and improves the accuracy of individualized neural rehabilitation.

Description

Technical Field

[0001] This invention relates to the field of intelligent rehabilitation and brain-computer interface technology, specifically to a brain-computer interaction-based neurorehabilitation training method for patients based on digital twins. Background Technology

[0002] Motor dysfunction caused by neurological injury is one of the leading causes of disability worldwide. Modern neurorehabilitation theory indicates that the closed-loop integration of repeated, active motor intention activation and sensory feedback is the core mechanism for promoting neural plasticity remodeling and restoring motor function. Brain-computer interface (BCI) technology, by acquiring and decoding the patient's brain neural activity in real time, transforms motor intention into control commands for external devices and combines this with functional electrical stimulation (FES) to provide sensory feedback, has become one of the most promising intervention methods in the field of neurorehabilitation.

[0003] Existing brain-computer interface-based rehabilitation systems suffer from the following core bottlenecks:

[0004] First, the targeting of neural signal processing is insufficient. Existing systems generally employ broadband filtering strategies, which easily preserve baseline drift and electrooculography artifacts in the low-frequency band, while introducing electromyographic noise in the high-frequency band, severely interfering with the extraction quality of the core frequency bands of motor imagery (μ rhythm and β rhythm). Traditional methods rely on pre-stored fixed artifact templates for online removal, which completely fails when faced with dynamic random artifacts such as head micro-movements, swallowing, and limb movements commonly found in rehabilitation patients, lacking adaptive removal capabilities.

[0005] Secondly, individual neural models are lacking, and training data is scarce. Motor intention decoding models rely on a large amount of labeled data for training, but neurological patients, limited by the degree of their motor impairment, often find it difficult to complete standard calibration training paradigms, resulting in a severe shortage of labeled samples. Existing data augmentation schemes lack neurophysiological constraints, and the generated data deviates significantly from the actual EEG distribution of patients, making it difficult to effectively improve decoding performance.

[0006] Third, the training program lacks a quantitative understanding of the time window of neuroplasticity. The intensity and frequency of rehabilitation training are not dynamically adjusted according to the patient's neuroplasticity state, and the timing of training is not accurately grasped, affecting rehabilitation efficiency and safety.

[0007] Based on this, the present invention designs a brain-computer interface-based neurorehabilitation training method for patients with digital twins to solve the above problems. Summary of the Invention

[0008] In view of the above-mentioned shortcomings of the existing technology, the present invention provides a brain-computer interface-based neurorehabilitation training method for patients based on digital twins.

[0009] To achieve the above objectives, the present invention provides the following technical solution:

[0010] A brain-computer interface-based neurorehabilitation training method for patients based on digital twins includes the following steps:

[0011] Step S1: Collect multi-channel electroencephalogram (EEG) signals and functional near-infrared spectral signals from neurological patients, and form neural signal data frames after adaptive hierarchical preprocessing;

[0012] Step S2: Using the neural signal data frames as the driving force, construct an active digital twin of the patient's nervous system based on a neural quality network model with multiple cortical column nodes; the active digital twin includes a short-term neural state and parameter joint estimation module and a long-term neural plasticity evolution module.

[0013] Step S3: Based on the neural parameter evolution sequence output by the neural plasticity evolution module, quantify and predict the patient's neural plasticity window period; apply controlled virtual perturbation to the model parameter estimates output by the neural state and parameter joint estimation module within the neurophysiological constraint domain, and combine the spatiotemporal activation template of the motor imagination brain network to drive the positive simulation of the active digital twin to generate an labeled synthetic EEG dataset with individual neurophysiological characteristic constraints;

[0014] Step S4: The motor intention decoding network is virtually pre-trained using the synthetic EEG dataset to generate individual adaptation weights. Based on the adaptation weights, the neural signal data frames are decoded to generate motor intentions. The decoded motor intention categories are converted into control commands to drive the external rehabilitation device to perform corresponding actions.

[0015] Step S5: Based on the neural state estimate and model parameter estimate output by the neural state and parameter joint estimation module, and combined with the simulated and predicted response of the somatosensory cortex to different functional electrical stimulation parameter combinations in the active digital twin, the functional electrical stimulation parameter combination that makes the predicted somatosensory cortex activation response intensity tend towards the optimization target value is determined through optimization search. The parameter combination is sent to an external rehabilitation device to apply sensory feedback, and the measured somatosensory cortex response is fed back to the active digital twin for error correction. According to the neural plasticity window period, the training intensity and training frequency of subsequent training are dynamically adjusted.

[0016] Furthermore, cascaded bandpass filters of 8-13Hz and 14-30Hz were applied sequentially to the EEG signals to extract the μ-rhythm and β-rhythm neural oscillation components, respectively. An online adaptive sparse component analysis (SCA) algorithm was used to update the artifact dictionary in real time with a sliding time window to adaptively identify and remove dynamic random artifacts. After z-score normalization of each channel, sliding window framing was performed with a preset frame length and inter-frame overlap rate. The functional near-infrared spectral signal was calculated based on the modified Beer-Lambert law to calculate the changes in oxyhemoglobin and deoxyhemoglobin concentrations and then bandpass filtered. The processed EEG data frames and functional near-infrared spectral data frames were merged into neural signal data frames after time alignment.

[0017] Furthermore, based on the Jansen-Rit type neural quality model, an augmented state vector is jointly constructed by the state vectors of each node, the model parameters within the node, and the connection weights between nodes. Based on the unscented Kalman filter algorithm, the augmented state vector is subjected to online state-parameter joint estimation to obtain the posterior estimation sequence of short-term neural states and parameters.

[0018] The median and interquartile range statistics of intra-node synaptic gain parameters and inter-node connection weight parameters are extracted from the posterior estimation sequence for each training segment to obtain candidate long-term parameter evolution sample points. For each sample point, its Mahalanobis distance anomaly score is calculated. After removing abnormal samples, they are included in the long-term parameter evolution time series. Trend fitting and extrapolation prediction are performed on the parameter evolution time series to achieve weekly-monthly tracking and prediction of neural plasticity evolution.

[0019] Furthermore, using the unscented Kalman filter algorithm, the state variable vector and model parameter vector of the neural mass model are combined to form an augmented state vector. Based on the augmented state vector, a set of Sigma points is generated. The set of Sigma points is then predicted in one step using a discretized nonlinear state transition function and then weighted and summed to obtain the prior state estimate and prior covariance. Subsequently, the prior predicted Sigma points are mapped to the observation space through an observation function. The Kalman gain is calculated using the measured EEG signal at the current time as the observation, and the prior estimate is corrected posteriorly. This process is iterated to track the brain state.

[0020] Furthermore, synaptic connection strength parameters are extracted from the model parameter estimates for each training segment to form a long-term parameter evolution time series. The time series is then subjected to trend fitting and extrapolation prediction using an autoregressive moving average (ARMA) model. When the rate of change of the excitation-inhibition ratio, represented by the predicted excitatory synaptic gain parameters, inhibitory synaptic gain parameters, and corresponding connection strength coefficients, exceeds a preset threshold, the patient is determined to have entered the neuroplasticity window period, and training reinforcement suggestions are issued to the clinician.

[0021] Furthermore, Monte Carlo sampling is performed within the neighborhood of the posterior distribution of the parameter estimation mean to generate multiple sets of perturbation parameter vectors. Each set of perturbation vectors follows a Gaussian distribution with the zero vector as the mean and the perturbation variance adaptively determined according to the diagonal elements of the parameter posterior covariance matrix as the scale. The thalamic cortex input mean corresponding to the three types of motor category labels (left-hand motor imagery, right-hand motor imagery, and resting) is used as the external driver. Based on the neurophysiological laws of motor imagery, the inter-node connection weight modulation mode of enhanced contralateral sensorimotor event-related desynchronization (ERD) and enhanced ipsilateral event-related synchronization (ERS) is superimposed to drive the positive simulation of the active digital twin, generating a multi-channel synthetic EEG time series carrying motor category labels. After screening, a synthetic EEG dataset is formed.

[0022] The screening process includes: calculating the deviation between the synthetic EEG dataset and the current real EEG data frames of patients in terms of power spectral density in the μ / β band and Frechet feature distance; when the deviation exceeds a preset threshold, shrinking the perturbation variance scale and resampling until the dataset passes the distribution consistency check and is included in the dataset.

[0023] Furthermore, the confidence score of each decoding output is calculated, and only neural signal data frames with confidence scores exceeding a preset gating threshold and their corresponding decoding labels are included in the fine-tuning sample pool. The samples in the fine-tuning sample pool are used as supervision signals to perform small-batch incremental updates on the individual adaptation weights. At the same time, the individual adaptation weights generated by the virtual pre-training are used as regularization anchors to apply parameter drift penalties, thereby realizing the dynamic evolution of the individual adaptation weights with the patient's neural state.

[0024] Furthermore, during the patient's formal rehabilitation training, the confidence level of the motor intention decoding results is continuously assessed. Samples corresponding to decoding results where the component with the highest value in the classification probability vector exceeds a preset threshold are pseudo-labeled. Once the preset threshold is reached, an online fine-tuning of the individual-specific adaptation layer is triggered, enabling the decoding network to continuously converge towards the patient's current neural state.

[0025] Furthermore, the peak amplitude and latency of somatosensory evoked potentials in the corresponding leads of the somatosensory cortex after the application of functional electrical stimulation were extracted as measured response indicators. The measured response indicators were subtracted from the twin predicted response indicators to form a residual vector. Using the residual vector as the observation, a heavy-tailed noise robust Kalman filter mechanism was used to perform posterior correction on the synaptic gain parameters and cortical excitability parameters of the somatosensory cortex sub-model, and the corrected parameters were incorporated into the neural parameter evolution sequence.

[0026] Furthermore, based on the mean ratio of the excitatory to inhibitory synaptic gain estimates of each cortical node in the model parameter estimates, when the rate of change of the excitatory-inhibition ratio is lower than the plasticity activation threshold, the current training intensity and frequency, as well as the previous functional electrical stimulation parameters, remain unchanged; when the rate of change of the excitatory-inhibition ratio exceeds the plasticity activation threshold, the duration of a single training session is increased to a preset threshold and the training interval is shortened to a preset minimum safe interval, while triggering full parameter space re-optimization of the functional electrical stimulation parameters; when the mean decoding confidence or the rate of prolongation of somatosensory evoked potential latency exceeds a preset fatigue warning threshold, the training intensity is reduced and the rest interval is extended, while a fatigue warning is sent to the clinician.

[0027] Compared with the prior art, the beneficial effects of this invention are as follows:

[0028] 1. By constructing a digital twin of the patient's active nervous system, we can achieve real-time tracking of brain state and precise dynamic generation of sensory feedback parameters, providing a technical basis for activating bidirectional Hebbian plasticity and improving the efficiency of neural pathway reconstruction.

[0029] 2. By modeling the evolution of synaptic parameters over a long timescale and using autoregressive moving average prediction, we can quantitatively identify the window of neural plasticity, providing clinicians with an objective basis for deciding on the timing of training and avoiding a decline in rehabilitation effectiveness due to biased training timing.

[0030] 3. By generating virtual perturbation synthetic data driven by digital twins and pre-training with zero-sample individual adaptation, the mandatory calibration step before training of traditional brain-computer interface systems is eliminated, reducing the system activation time from tens of minutes to minutes, and significantly improving the system availability for patients for whom traditional decoding methods have failed. Attached Figure Description

[0031] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are merely some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without any creative effort.

[0032] Figure 1 This is a flowchart of the brain-computer interface-based neurorehabilitation training method for patients based on digital twins according to the present invention. Detailed Implementation

[0033] 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. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0034] Example 1: Figure 1 This invention provides a brain-computer interface-based neurorehabilitation training method for patients based on digital twins, comprising the following steps:

[0035] Step S1: Collect multi-channel electroencephalogram (EEG) signals and functional near-infrared spectral signals from neurological patients, and form neural signal data frames after adaptive hierarchical preprocessing.

[0036] Step S1 includes the following steps:

[0037] Step 1.1: Using standard electrode caps, multi-channel electroencephalography (EEG) signals were acquired from the neurological patient, focusing on leads corresponding to the bilateral sensorimotor cortex regions. Given that the core neural oscillation components of the motor imagery task are concentrated in the μ rhythm (8–13 Hz) and β rhythm (14–30 Hz), cascaded bandpass filters of 8–13 Hz and 14–30 Hz were sequentially applied to the raw EEG signals to extract the μ rhythm and β rhythm neural oscillation components, respectively.

[0038] Step 1.2: Regarding artifact removal, for dynamic random artifacts commonly found in rehabilitation patients, such as head micro-movements, swallowing, and voluntary limb movements, an online adaptive sparse component analysis (SCA) algorithm is employed. This algorithm updates the artifact dictionary in real-time using a sliding time window to adaptively identify and remove dynamic random artifacts. Compared to traditional independent component analysis methods that rely on pre-stored fixed templates, this approach continuously tracks the time-varying statistical characteristics of artifacts, ensuring the effectiveness of artifact removal under diverse movement states of rehabilitation patients. After artifact removal for each channel, z-score normalization is performed to eliminate baseline differences between leads.

[0039] The resting power of patients in a relaxed or closed-eye resting state is collected as the resting threshold, and a non-uniform dynamic sliding window mechanism driven by neural state perception is introduced: the event-related desynchronization index of μ / β rhythm in each channel is calculated in real time. ,when When the threshold is below the resting state, the patient is determined to be in a neural resting state. At this time, the window step size is automatically expanded to 100ms, and the inter-frame overlap is reduced to 60% to reduce redundant computation. If the threshold for resting motion is exceeded and continues to rise, the patient is determined to be in a state of active motor intention. The step size is adaptively compressed to 25ms, and the inter-frame overlap is increased to 90% to achieve precise capture of the moment of sudden change in motor intention. The frame length is kept constant at 250ms, and the frequency resolution of time-frequency analysis is 4Hz to distinguish the boundary between μ rhythm and β rhythm.

[0040] This strategy enables dynamic matching of data frame density and neural activation level, ensuring high temporal resolution during key periods of motion intent while effectively alleviating real-time processing pressure and avoiding high data redundancy and real-time pressure caused by excessive overlap.

[0041] Step 1.3: For the functional near-infrared spectral signal, calculate the changes in oxyhemoglobin and deoxyhemoglobin concentrations based on the modified Beer-Lambert law, and perform bandpass filtering to extract hemodynamic response components.

[0042] The modified Beer-Lambert law formula is as follows:

[0043]

[0044] in, This represents the change in hemoglobin concentration. The change in optical density It represents the molar extinction coefficient of oxyhemoglobin or deoxyhemoglobin at a specific wavelength. The source-detector distance between the emitter and detector electrodes. It is the differential path factor.

[0045] It is used to correct the deviation between the actual optical path and the source-detector distance caused by light scattering in the tissue, and its value is determined based on the tissue optical characteristics of the detection site.

[0046] Step 1.4: Merge the processed EEG data frames and functional near-infrared spectroscopy data frames after time alignment to form neural signal data frames. Time alignment is achieved by using a unified sampling timestamp to ensure strict correspondence between the two modalities in the time dimension.

[0047] Step S2: Using the neural signal data frames as the driving force, construct an active digital twin of the patient's nervous system based on a neural quality network model with multiple cortical pillar nodes; the active digital twin includes a short-term neural state and parameter joint estimation module and a long-term neural plasticity evolution module.

[0048] Step S2 includes the following steps:

[0049] Step 2.1: Simultaneously perform dynamic modeling and parameter estimation of the patient's nervous system state on two time dimensions: short time scale (milliseconds to seconds) and long time scale (weeks to months).

[0050] Each functional region, including the sensorimotor cortex, premotor cortex, supplementary motor area, and thalamus, corresponds to a cortical column node. Within each node, the Jansen-Rit neural quality model is used to describe the dynamic synaptic interactions between excitatory pyramidal cell groups, excitatory interneuron groups, and inhibitory interneuron groups. The nodes are coupled through a directed connection weight matrix to form a multi-node neural quality network model.

[0051] The nodal state equations are:

[0052]

[0053] in, It is a set of Jansen-Rit type nonlinear differential equations. To convert the average membrane potential into a Sigmoid-type nonlinear function of the average discharge rate, To act on the first Input signals from the outer cortex of the node. This represents the derivative of the state vector with respect to time. The model parameter vector within the node. For from the first Node to the The strength of incoming connections to a node. For the first Each cortical column node at time [time] The state vector.

[0054] Specifically, no. Each cortical column node at time [time] The state vector is This includes the average membrane potential and its first-order time derivative of each neuron group at the node; the intra-node model parameter vector is... Includes excitatory synaptic gain Inhibitory synaptic gain Excitatory synaptic time constant Inhibitory synaptic time constant The connection weight matrix between nodes is denoted as... Its elements Indicates from the first Node to the The strength of incoming connections to a node.

[0055] The theoretical observation output of the model is:

[0056]

[0057] in, This involves concatenating the state vectors of all nodes. To map the model state to the observation matrix in the EEG observation channel space, This is the observed noise vector.

[0058] To achieve online real-time tracking of individual patient neural states and model parameters, the state vectors of each node are... Model parameters within nodes and the connection weights between nodes Jointly construct augmented state vectors:

[0059]

[0060] in, This involves concatenating the state vectors of all nodes. This involves concatenating the model parameter vectors within all nodes. To represent the vectorized form of the connection weight matrix, the superscript... This indicates the transpose operation.

[0061] Based on the Unscented Kalman Filter (UKF) algorithm, the augmented state vector is... Perform online state-parameter joint estimation:

[0062]

[0063]

[0064]

[0065] in, To augment the dimension of the state vector, For composite scaling parameters, As auxiliary scaling parameters, The first square root of a matrix represents the square root of the matrix. List, For the prior mean of the augmented state vector, The prior covariance matrix, .

[0066] Each Sigma point is predicted in one step using a discretized nonlinear state transition function, and then weighted and summed to obtain the prior state estimate and prior covariance. The prior predicted Sigma points are then mapped to the observation space using an observation function, based on the measured EEG signal at the current time. Calculate Kalman gain for observations Posterior correction is applied to the prior estimate:

[0067]

[0068]

[0069] in, For prior prediction of observations, The new information covariance matrix, The mean of the posterior state estimate. Let be the posterior covariance matrix.

[0070] By iterating the above process, the patient's brain state is tracked in real time to obtain a posterior estimate sequence of short-term neural states and parameters.

[0071] On a per-training-segment basis, extract intra-node synaptic gain parameters from the posterior estimation sequence. , The median and interquartile range statistics of the connection weight parameters between nodes are used as candidate long-term parameter evolution sample points. For each sample point, its Mahalanobis distance anomaly score is calculated. After removing anomalies, the valid sample points are included in the long-term parameter evolution time series.

[0072] Mahalanobis distance anomaly score is defined as:

[0073]

[0074] in, For the first The parameter statistics vector extracted from each training segment This is the mean vector of the historical samples. The covariance matrix of the historical samples, For sample points Mahalanobis distance relative to the historical distribution center.

[0075] The autoregressive moving average (ARIMA) model is applied to the cleaned long-term parameter evolution time series for trend fitting and extrapolation prediction, enabling weekly-monthly tracking of neural plasticity evolution.

[0076] Step S3: Based on the neural parameter evolution sequence output by the neural plasticity evolution module, quantify and predict the patient's neural plasticity window period; apply controlled virtual perturbation to the model parameter estimates output by the neural state and parameter joint estimation module within the neurophysiological constraint domain, and combine the spatiotemporal activation template of the motor imagination brain network to drive the positive simulation of the active digital twin to generate an labeled synthetic EEG dataset with individual neurophysiological characteristic constraints.

[0077] Step S3 includes the following steps:

[0078] Step 3.1

[0079] At the end of each training segment, the excitatory synaptic gain parameters of each cortical node are extracted from the posterior estimation sequence. With inhibitory synaptic gain parameters Calculate the excitation-inhibition (E / I) ratio:

[0080]

[0081] in, For the first The number of cortical nodes in the first The excitation-inhibition ratio of each training segment, This is an estimate of the excitatory synaptic gain. This is an estimate of the inhibitory synaptic gain.

[0082] The global excitation-inhibition ratio is obtained by averaging over all nodes. Furthermore, a comprehensive E / I dynamic index sequence is constructed by combining the connection strength coefficients between each node, and incorporated into the long-term parameter evolution time series.

[0083] The ARIMA model was applied to the time series of the above E / I dynamic index to perform trend extrapolation and calculate the predicted rate of change of the E / I ratio:

[0084]

[0085] in, The E / I ratio at time step The rate of change within, It represents the difference in the E / I ratio between adjacent training segments.

[0086] When the predicted rate of change of the E / I ratio Exceeding the preset plasticity activation threshold When the system determines that the patient has entered the neuroplasticity window, it sends training reinforcement suggestions to the clinician and triggers dynamic adjustment of training parameters.

[0087] Before generating synthetic EEG data, a parameter constraint domain was first established based on the clinical neurophysiological value ranges of each parameter in the Jansen-Rit model. Specifically, the parameter vectors of each cortical node were... Define the lower bound vector for each component. With upper bound vector This ensures that the parameters generated by subsequent perturbations always satisfy:

[0088]

[0089] The introduction of this constraint domain ensures that the generated perturbation parameters do not exceed the physiologically reasonable range of the human nervous system, thus avoiding the generation of invalid pseudo-EEG data that violates the laws of neurophysiology.

[0090] The posterior covariance matrix of parameters obtained based on unscented Kalman filtering The diagonal elements are extracted as reference variances for each parameter dimension, thereby adaptively determining the scale of the perturbation variance. Monte Carlo sampling is performed within the neighborhood of the posterior distribution of the parameter estimation to generate multiple sets of perturbation parameter vectors:

[0091]

[0092] in, This represents the posterior mean of the parameters at the end of the current training segment. For the first The perturbation vectors are such that their components follow a Gaussian distribution with a mean of zero and a variance equal to an adaptive scaling factor of the posterior variance of the corresponding parameters.

[0093]

[0094] in, For the first The first group of perturbation vectors Quantity, The posterior covariance matrix of the parameters is the first... diagonal elements, This is the disturbance amplitude control coefficient. It follows a normal distribution.

[0095] The above-mentioned neurophysiological constraint domain pruning is applied to all sampling parameters to ensure the physiological rationality of the parameters.

[0096] For the three categories of motor imagery—left-hand motor imagery, right-hand motor imagery, and resting motor imagery—corresponding spatiotemporal activation templates were constructed based on the neurophysiological laws of motor imagery to ensure the physiological reliability of the synthesized spatiotemporal features of EEG.

[0097] Specifically: for left-hand motor imagery categories, the corresponding thalamic cortex input mean was used. As an external driver, an event-related desynchronization (ERD) enhancement mode is superimposed on the right (contralateral) sensoromotor area, which modulates the down-modulation of the connection weights between nodes in the right sensoromotor area, and an event-related synchronization (ERS) enhancement mode is superimposed on the left (same side) sensoromotor area.

[0098] Category of right-hand movement imagination, As an externally driven, the internode connection weight modulation mode is superimposed with the ERD enhancement of the left (contralateral) sensorimotor area and the ERS enhancement of the right (same side) sensorimotor area.

[0099] For the resting category, As an externally driven system, the connection weights between nodes maintain the basic estimated values ​​without applying additional modulation.

[0100] Based on each set of disturbance parameters Using spatiotemporal activation templates of the corresponding categories, the active digital twin is driven to perform forward simulation, generating a multi-channel synthetic EEG time series with motion category annotations.

[0101] Calculate the deviation of the synthetic EEG time series from the current patient's real EEG data frames in the μ / β band power spectral density, and the Frechet distance between their characteristic distributions:

[0102]

[0103] in, The Frechet distance is the difference between the synthetic EEG and the actual EEG feature distribution. This is the mean vector of the synthetic EEG feature distribution. For the covariance matrix, This is the mean vector of the feature distribution of real EEG data frames. It is the covariance matrix; This is the trace operation of a matrix.

[0104] When the power spectral density deviation or Frechet distance exceeds the preset threshold, the contraction perturbation amplitude control coefficient The samples were then resampled until they passed the distribution consistency check before being included in the dataset, thus forming a labeled synthetic EEG dataset with individual neurophysiological characteristics as constraints.

[0105] Step S4: The motor intention decoding network is virtually pre-trained using the synthetic EEG dataset to generate individual adaptation weights. Based on the adaptation weights, the neural signal data frames are decoded to generate motor intentions. The decoded motor intention categories are converted into control commands to drive external rehabilitation devices to perform corresponding actions.

[0106] Step S4 includes the following steps:

[0107] The motor intention decoding network was virtually pre-trained using the aforementioned synthetic EEG dataset to optimize network weights to adapt to the individual neurophysiological characteristics of the target patients, thereby generating individually adapted weights. This process eliminates the need for patients to undergo real-world annotation training, effectively addressing the problem of scarce labeled data for neurological patients.

[0108] During formal rehabilitation training, based on individual adaptation weights Decode the motion intent of the neural signal data frame and output a classification probability vector:

[0109]

[0110] in The total number of sports categories, For the first Predicted probability of class This is the classification probability vector.

[0111] Specifically, Take the maximum value This serves as the confidence score for this decoding process.

[0112] Given the generally low decoding confidence levels in the early stages of patient recovery, a confidence threshold is required. Adaptive setting strategy:

[0113]

[0114] in, For the first The current confidence gating threshold for the training segment. The preset minimum security threshold, This represents the mean of the maximum predicted probabilities within the recent sliding window. The corresponding standard deviation is... An adaptive coefficient for controlling the tightness of the gate.

[0115] Only those with confidence levels exceeding the current threshold will be considered. The neural signal data frames and their corresponding decoding labels are included in the fine-tuning sample pool to avoid low-quality samples contaminating the online update process.

[0116] Using samples from the fine-tuning pool as supervision signals, the individual fitness weights are updated incrementally in small batches. A parameter drift penalty term is introduced into the update objective function, using virtual pre-trained individual fitness weights. For regularization anchor points:

[0117]

[0118] in, For the total loss function, Based on the current weight The cross-entropy classification loss is calculated. The regularization coefficient is . The squared Frobenius norm of the current weights relative to the pre-trained anchor weights. As the current weight, Weights are assigned to individuals.

[0119] Once the accumulated pseudo-labeled samples reach the preset batch size, in addition to triggering the small-batch fine-tuning of the adaptation layer, the system simultaneously evaluates the overall network performance. When the mean decoding confidence of several consecutive training segments shows a continuous downward trend, it indicates that the patient's neurological state has undergone a long-term and significant change. At this point, periodic deep iterative optimization of the network backbone layer is triggered: the network backbone is re-fine-tuned with all parameters using the latest accumulated synthetic EEG dataset to prevent the continuous degradation of decoding performance caused by the long-term disconnect between the general features of the backbone and the patient's real features, thus achieving long-term continuous optimization of the decoding network.

[0120] Step S5: Based on the neural state estimate and model parameter estimate output by the neural state and parameter joint estimation module, and combined with the simulated and predicted response of the somatosensory cortex to different functional electrical stimulation parameter combinations in the active digital twin, the functional electrical stimulation parameter combination that makes the predicted somatosensory cortex activation response intensity tend towards the optimization target value is determined through optimization search. The parameter combination is sent to an external rehabilitation device to apply sensory feedback, and the measured somatosensory cortex response is fed back to the active digital twin for error correction. According to the neural plasticity window period, the training intensity and training frequency of subsequent training are dynamically adjusted.

[0121] Step S5 includes the following steps:

[0122] Before rehabilitation training begins, a standardized, stepwise, incremental stimulation protocol is used to assess the individual sensory thresholds of each patient. Calibration: Starting with a safe initial current intensity, gradually increase the stimulation current with a fixed step size, record the current intensity at which the patient first reports sensation, repeat three times, and take the average value as the initial sensory threshold. In subsequent training, by combining real-time estimation of the peak amplitude and latency of somatosensory evoked potentials (SEPs) and based on the prediction of the somatosensory cortex activatability by the active digital twin, the sensory threshold is periodically and dynamically updated to maintain the current effective sensory threshold. This ensures that the current intensity of functional electrical stimulation is always limited to a safe range, avoiding excessive stimulation that could cause stinging or muscle spasms, or ineffective stimulation due to threshold changes.

[0123] Based on the current neural state estimate and model parameter estimate output by the joint neural state and parameter estimation module, and combined with the simulation of the somatosensory cortex sub-model in an active digital twin, the system predicts the response of the somatosensory cortex to different combinations of functional electrical stimulation parameters (including stimulation frequency). Pulse width and current intensity ) activation response strength .

[0124] The optimization search was used to determine the target value that would make the predicted intensity of the somatosensory cortex activation response tend towards the optimization target value. Functional electrical stimulation parameter combinations:

[0125]

[0126]

[0127] in, As constraints, For safety margin coefficient, The current dynamic sensing threshold, For the optimal combination of stimulus parameters, For stimulation frequency, The pulse width. For current intensity, To optimize the target activation strength, This is a predicted value for the intensity of the somatosensory cortex activation response. The independent variable is the one that minimizes the value.

[0128] Following the application of functional electrical stimulation, somatosensory evoked potential (SEP) responses in the corresponding leads of the somatosensory cortex were extracted using an adaptive time window. To address individual differences in somatosensory cortex response latency among patients, the start time of the response acquisition window was based on the average historical SEP latency of the current patient. Adaptive settings, rather than using a fixed time window, are used to avoid missing the true response due to individual differences in response latency.

[0129] Extracting the peak amplitude of somatosensory evoked potentials and incubation period As a measured response index, it constitutes the measured response vector. .

[0130] The measured response vector Predicting response vectors with active digital twins The difference is calculated to form the residual vector. .

[0131] To distinguish between errors caused by signal noise and errors caused by actual changes in neural state, and to avoid wasting computing power due to frequent and meaningless updates, a residual significance test is introduced into the error correction trigger judgment: only when the weighted norm of the residual vector exceeds a preset noise tolerance threshold... At that time, twin parameter correction is triggered, that is:

[0132]

[0133] in, The noise covariance matrix obtained from the observation noise estimation. For The weighted norm of the weight matrix. To preset the noise tolerance threshold, This is the transpose of the residual vector. For the residual vector, This is the inverse matrix of the observation noise covariance.

[0134] When the triggering condition is met, the residual vector is used. For observation purposes, a heavy-tailed noise robust Kalman filter mechanism was used to perform posterior correction on the synaptic gain parameters and cortical excitability parameters of the somatosensory cortex sub-model. The corrected parameters were then incorporated into the neural parameter evolution sequence to continuously improve the modeling accuracy of the somatosensory cortex dynamics by the active digital twin.

[0135] Based on the rate of change of the global excitation-inhibition ratio output by the active digital twin Mean confidence score for decoding motion intent and the latency prolongation rate of somatosensory evoked potentials The subsequent training plan will be dynamically adjusted in three levels:

[0136] Stable maintenance state: when That is, the rate of change of the E / I ratio is lower than the plasticity activation threshold, and and All values ​​are within the normal range. The system determines that the patient's current neurological state is stable and maintains the current training intensity and frequency as well as the previous functional electrical stimulation parameters to consolidate the existing training effect.

[0137] Plasticity window enhancement state: when When the rate of change of the E / I ratio exceeds the plasticity activation threshold, the system determines that the patient has entered the neuroplasticity window period, increases the duration of a single training session to the preset reinforcement duration threshold, and shortens the training interval to the preset minimum safe interval. At the same time, it triggers full-parameter spatial re-optimization of functional electrical stimulation parameters to maximize the neuroremodeling benefits of rehabilitation training within the most active time window of neuroplasticity.

[0138] Fatigue warning status: when the mean confidence level of the decoded value is... The rate of decline to below the preset fatigue warning threshold, or the rate of increase in the latency of somatosensory evoked potentials. When the preset fatigue warning threshold is exceeded, the system determines that the patient is showing signs of fatigue, reduces the training intensity and extends the rest interval, and sends a fatigue warning to the clinician to ensure the patient's training safety.

[0139] The latency extension rate of somatosensory evoked potentials is defined as:

[0140]

[0141] in, The latency extension rate of somatosensory evoked potentials in the current training segment. This represents the measured latency of somatosensory evoked potentials during the current training segment. The baseline latency of somatosensory evoked potentials for each individual patient.

[0142] Example 2: Comparison of the effects of the present invention with existing methods.

[0143] Following the methodology of the recent BCI Rehabilitation Trial (Cheng et al., J NeuroEngineeringRehabil, 2025, 22:45), the Fugl-Meyer Upper Extremity Scale (FMA-UE) change over eight weeks was used as the endpoint. Based on the pre-trial effect size (Cohen's d≈1.4), α=0.05, and a statistical power of 80%, it was estimated that at least 22 cases were needed in each group. To support subgroup analysis and consider an approximately 15% dropout rate, 96 cases were actually enrolled. They were stratified and randomly assigned to an experimental group (n=48) and a control group (n=48) based on four factors: disease type (stroke / traumatic brain injury), duration (subacute 1-6 months / chronic 6-18 months), injury location (cortical / subcortical), and age (<55 / ≥55 years). The baselines of the two groups were balanced (mean FMA-UE 32.6 vs 33.1, p>0.1). Of the 72 patients, 24 had traumatic brain injury; 46 had subacute stroke and 50 had chronic stroke; 38 had cortical involvement and 58 had subcortical involvement; 42 were under 55 years old and 54 were 55 years old or older. Both groups received 8 weeks of training at the same frequency, 5 times a week for 45 minutes each time. The primary endpoint was the change in FMA-UE, and the secondary endpoints were the Wolf Motor Function Test (WMFT) time and the modified Barthel Index (MBI).

[0144] Existing systems such as Cerebus, recoveriX, and "Shen Gong-Ling Xi Zhi" only provide visual feedback, with stimulation parameters set to preset fixed values ​​and adjustment cycles exceeding several minutes. The feedback delay is far greater than the synaptic plasticity time window (approximately 20-100ms), making it difficult to activate Hebbian long-term enhancement. This invention utilizes digital twins to estimate neural states in real time and Bayesian optimization to dynamically generate somatosensory stimulation parameters within hundreds of milliseconds. This compresses the closed-loop delay of "intent decoding → motor execution → somatosensory feedback" to within 50ms (decoding ≤ 20ms, optimization ≤ 100ms, stimulus delivery ≤ 10ms). For the first time, the timing requirements of bidirectional Hebbian plasticity are achieved in engineering, with synchronous triggering of the proprioceptive and visual channels, significantly improving integrity.

[0145] Existing methods rely on bi-weekly scale assessments, which lag behind neurophysiological changes. This invention uses long-term UKF to continuously estimate excitation / inhibition gain and connectivity coefficients, forming an evolutionary sequence of the E / I ratio for successive training, and extrapolates and predicts using an ARMA model. When the rate of change in the E / I ratio exceeds a threshold, it is determined that the patient has entered a plasticity window. In 96 patients, this method identified the recovery acceleration point an average of 5.7 days earlier than scale assessment (95% CI 4.2–7.2 days), 5.3 days in the stroke subgroup, and 6.9 days in the traumatic brain injury subgroup. During the window, the average monthly improvement in FMA-UE in the intensive training subgroup was 3.0 times that of the subgroup outside the window (p<0.001), effectively avoiding plateaus, and this advantage was consistent across all strata.

[0146] Traditional BCI calibration takes 15-30 minutes (average measured 22.6 min), and fatigue leads to a decline in data quality, with approximately 15-30% of "BCI illiterate" users achieving decoding accuracy below 70%. This invention utilizes digital twin virtual pre-training: generating a synthetic EEG fine-tuning decoding network based on the mean UKF parameters, completing the entire process automatically in 2-5 minutes (average 3.8 min), requiring no patient intervention. In 96 cases, the zero-shot decoding accuracy reached 78.6%, significantly higher than the 73.8% of traditional CSP (p<0.001); among 19 illiterate users (16 with stroke, 3 with traumatic brain injury), 17 achieved a usable level of over 70% (89.5%), with consistent advantages across different disease, disease duration, location, and age subgroups. The effective duration of a single training session increased from 62% to 94%.

[0147] Primary endpoints for the full sample: The experimental group showed an improvement of 15.8 ± 3.2 points in FMA-UE over eight weeks, while the control group showed an improvement of 10.1 ± 3.5 points, with a difference of 5.7 points between the groups (p < 0.001, Cohen's d = 1.70); the improvement rate in WMFT was 52.7% vs 26.3%; and the improvement in MBI was 23.1 vs 13.8 points.

[0148] Disease stratification: Stroke experimental group (n=36) improved by 15.5 points vs. control group by 10.2 points (Δ5.3, p<0.001); Traumatic brain injury experimental group (n=12) improved by 16.7 points vs. control group by 10.0 points (Δ6.7, p<0.001). The interaction p=0.24, and there was no modifying effect on disease type.

[0149] Disease course stratification: The subacute experimental group (n=23) improved by 17.4 points vs 10.9 points (Δ6.5); the chronic experimental group (n=25) improved by 14.4 points vs 9.4 points (Δ5.0), both p<0.001, indicating that the chronic phase still showed significant benefits.

[0150] Injury site stratification: The cortical experimental group (n=19) improved by 14.9 points vs. 9.5 points (Δ5.4); the subcortical experimental group (n=29) improved by 16.4 points vs. 10.5 points (Δ5.9), both p<0.001.

[0151] Age stratification: <55 years old experimental group (n=21) improved by 16.3 points vs 10.6 points; ≥55 years old experimental group 27 cases improved by 15.4 points vs 9.7 points, the difference between the groups was 5.7 points, and the interaction p=0.95.

[0152] All subgroups showed gains exceeding the minimum clinically important difference in FMA-UE (4.25 points), with no significant interaction, indicating that the invention is robust and effective in a broad population with different disease stages, locations, and ages in stroke / traumatic brain injury.

[0153] Neuroelectrophysiology: In the experimental group, the amplitude of the affected side μ rhythm ERD increased from 18.2% to 31.5% (compared to 22.1% in the control group); the α-phase synchronization index between the affected and healthy sides M1 increased from 0.28 to 0.41 (compared to 0.27 to 0.31 in the control group); the amplitude of the subcortical TMS-MEP wave increased from 0.38mV to 1.12mV (compared to 0.36mV in the control group), and the elicitance rate increased from 52.2% to 82.6%, indicating the recovery of corticospinal tract conduction.

[0154] Resting-state fMRI (48 cases completed, 24 cases in each group): The experimental group showed significantly enhanced M1-SMA / PMC functional connectivity; the sensorimotor network clustering coefficient changed from 0.19 to 0.27 (p<0.001), the whole-brain global efficiency changed from 0.16 to 0.22 (p<0.001), and the small-world property σ changed from 1.31 to 1.52; the control group showed no significant changes in any imaging indicators.

[0155] DTI (30 cases completed, 15 cases in each group): In the experimental group, the FA of the corticospinal tract changed from 0.38 to 0.43 (p=0.003), and the MD decreased from 0.89×10⁻³ to 0.82×10⁻³ mm² / s (p=0.01). There was no significant change in the control group.

[0156] Cross-index validation: Changes in FMA-UE were moderately to highly positively correlated with changes in the E / I ratio, global efficiency, and FA (r=0.51-0.68, p<0.001), pointing to the recovery of neural plasticity from multiple levels of behavioral, electrophysiological, and imaging studies.

[0157] To distinguish the independent contributions of each technical feature, an additional 40 stroke patients were recruited and randomly divided into four groups (n=10 per group) for a 4-week ablation trial:

[0158] A Closed-loop feedback only 6.7±2.1 Δ4.2 (p=0.002) B Window-only identification 5.9±1.9 Δ5.0(p<0.001) C Zero-sample decoding only 6.4±2.0 Δ4.5(p<0.001) D Full module (this invention) 10.9±2.4 —

[0159] The overall module group gain was significantly higher than that of each individual module group (63%-85% higher), while the differences between groups A, B, and C were not significant. This confirms that there is a synergistic effect among closed-loop feedback, window period identification, and zero-shot decoding. The overall effect is not caused by a single factor, and the attribution is clear.

[0160] Summary table of effect comparison

[0161] Feeling feedback type Visual-oriented Visual + Proprioceptive Dual-Channel Closed Loop Added proprioceptive channel Feedback parameter update cycle Minutes ≤50ms Meets Hebbian timing requirements Window period identification method Clinical scales (intervals of two weeks or more) Twin parameter evolution + ARMA prediction (per training session) On average, identification is performed 5.7 days in advance. Calibration time 15-30 minutes 2-5 minutes (virtual pre-training) Patient participation time reduced to zero BCI Illiteracy Availability* Approximately 40%-70% 89.5% (17 out of 19) Coverage has increased significantly FMA-UE Improvement in Eight Weeks 10.1 points 15.8 points An increase of 56% (d=1.70) Improvement in the subgroup of traumatic brain injury 10.0 points (n=12) 16.7 points (n=12) Effective across etiologies Chronic phase subgroup enhancement 9.4 points (n=25) 14.4 points (n=25) Significant benefits still exist in the chronic phase E / I ratio recovery 0.41→0.55 0.43→0.77 Approaching a healthy level fMRI global efficiency No significant changes 0.16→0.22（p<0.001） Gold Standard Evidence for Brain Network Remodeling DTI-FA Improvement No significant changes 0.38→0.43（p=0.003） Improved white matter microstructure Module contribution Undistinguished Ablation confirms synergistic gain The attribution is clear and there is no exaggeration.

[0162] The above embodiments, through a stratified randomized controlled trial of 96 patients (72 with stroke and 24 with traumatic brain injury), combined with ablation control experiments and a multimodal evaluation system, systematically verified the substantial progress of the active digital twin-driven closed-loop rehabilitation method of this invention compared to existing technologies. Stratified analysis showed that this invention exhibited robust and consistent rehabilitation gains in subgroups including stroke and traumatic brain injury, subacute and chronic phases, cortical and subcortical injuries, and different age groups. The consistency results of multimodal evaluation indicators (behavioral scales + neurophysiology + resting-state fMRI + DTI) provided multi-level evidence for the mechanism of neural plasticity. The ablation experiment clearly distinguished the independent contributions and synergistic gains of each technical module.

[0163] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions will not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A brain-computer interface-based neurorehabilitation training method for patients based on digital twins, characterized in that, The method includes: Step S1: Collect multi-channel electroencephalogram (EEG) signals and functional near-infrared spectral signals from neurological patients, and form neural signal data frames after adaptive hierarchical preprocessing; Step S2: Using the neural signal data frames as the driving force, construct an active digital twin of the patient's nervous system based on a neural quality network model with multiple cortical column nodes; the active digital twin includes a short-term neural state and parameter joint estimation module and a long-term neural plasticity evolution module. Step S3: Based on the neural parameter evolution sequence output by the neural plasticity evolution module, quantify and predict the patient's neural plasticity window period; apply controlled virtual perturbation to the model parameter estimates output by the neural state and parameter joint estimation module within the neurophysiological constraint domain, and combine the spatiotemporal activation template of the motor imagination brain network to drive the positive simulation of the active digital twin to generate an labeled synthetic EEG dataset with individual neurophysiological characteristic constraints; Step S4: The motor intention decoding network is virtually pre-trained using the synthetic EEG dataset to generate individual adaptation weights. Based on the adaptation weights, the neural signal data frames are decoded to generate motor intentions. The decoded motor intention categories are converted into control commands to drive the external rehabilitation device to perform corresponding actions. Step S5: Based on the neural state estimate and model parameter estimate output by the neural state and parameter joint estimation module, and combined with the simulated and predicted response of the somatosensory cortex to different functional electrical stimulation parameter combinations in the active digital twin, the functional electrical stimulation parameter combination that makes the predicted somatosensory cortex activation response intensity tend towards the optimization target value is determined through optimization search. The parameter combination is sent to an external rehabilitation device to apply sensory feedback, and the measured somatosensory cortex response is fed back to the active digital twin for error correction. According to the neural plasticity window period, the training intensity and training frequency of subsequent training are dynamically adjusted.

2. The brain-computer interface-based neurorehabilitation training method for patients based on digital twins according to claim 1, characterized in that, Step S1 specifically includes: The EEG signal was sequentially subjected to cascaded bandpass filtering of 8-13Hz and 14-30Hz to extract the μ-rhythm and β-rhythm neural oscillation components. An online adaptive sparse component analysis (SCA) algorithm was used to update the artifact dictionary in real time with a sliding time window to adaptively identify and remove dynamic random artifacts. After z-score normalization of each channel, sliding window framing was performed with a preset frame length and inter-frame overlap rate. The functional near-infrared spectral signal was calculated based on the modified Beer-Lambert law to calculate the changes in oxyhemoglobin and deoxyhemoglobin concentrations and then bandpass filtered. The processed EEG data frames and functional near-infrared spectral data frames were merged into neural signal data frames after time alignment.

3. The brain-computer interface-based neurorehabilitation training method for patients based on digital twins according to claim 1, characterized in that, In step S2, the method for constructing an active digital twin of the patient's nervous system based on a neural quality network model with multiple cortical column nodes, driven by the neural signal data frame, specifically includes: Based on the Jansen-Rit neural quality model, an augmented state vector is jointly constructed by the state vectors of each node, the model parameters within the node, and the connection weights between nodes. Based on the unscented Kalman filter algorithm, the augmented state vector is subjected to online joint state-parameter estimation to obtain the posterior estimation sequence of short-term neural states and parameters. The median and interquartile range statistics of intra-node synaptic gain parameters and inter-node connection weight parameters are extracted from the posterior estimation sequence for each training segment to obtain candidate long-term parameter evolution sample points. For each sample point, its Mahalanobis distance anomaly score is calculated. After removing abnormal samples, they are included in the long-term parameter evolution time series. Trend fitting and extrapolation prediction are performed on the parameter evolution time series to achieve weekly-monthly tracking and prediction of neural plasticity evolution.

4. The brain-computer interface-based neurorehabilitation training method for patients based on digital twins according to claim 3, characterized in that, The method for online state-parameter joint estimation of the augmented state vector based on the unscented Kalman filter algorithm specifically includes: Using the unscented Kalman filter algorithm, the state variable vector and model parameter vector of the neural mass model are combined to form an augmented state vector. Based on the augmented state vector, a set of Sigma points is generated. The set of Sigma points is then predicted in one step using a discretized nonlinear state transition function and then weighted and summed to obtain the prior state estimate and prior covariance. Subsequently, the prior predicted Sigma points are mapped to the observation space through an observation function. The Kalman gain is calculated using the measured EEG signal at the current time as the observation, and the prior estimate is corrected posteriorly. This process is iterated to track the brain state.

5. The brain-computer interface-based neurorehabilitation training method for patients based on digital twins according to claim 1, characterized in that, In step S3, the method for quantifying and predicting the patient's neuroplasticity window specifically includes: Synaptic connection strength parameters are extracted from the model parameter estimates for each training segment to form a long-term parameter evolution time series. The time series is then subjected to trend fitting and extrapolation prediction using an autoregressive moving average (ARMA) model. When the rate of change of the excitation-inhibition ratio represented by the predicted excitatory synaptic gain parameters, inhibitory synaptic gain parameters, and corresponding connection strength coefficients exceeds a preset threshold, the patient is determined to have entered the neuroplasticity window period, and training reinforcement suggestions are issued to the clinician.

6. The brain-computer interface-based neurorehabilitation training method for patients based on digital twins according to claim 1, characterized in that, In step S3, the method of applying controlled virtual perturbations to the model parameter estimates output by the neural state and parameter joint estimation module within the neurophysiological constraint domain specifically includes: Monte Carlo sampling is performed within the neighborhood of the posterior distribution of the parameter estimation mean to generate multiple sets of perturbation parameter vectors. Each set of perturbation vectors follows a Gaussian distribution with the mean of the zero vector and the perturbation variance adaptively determined by the diagonal elements of the parameter posterior covariance matrix. The thalamic cortex input mean corresponding to the three motor category labels of left-hand motor imagery, right-hand motor imagery, and resting motion is used as the external driver. Based on the neurophysiological laws of motor imagery, the inter-node connection weight modulation mode of enhanced contralateral sensorimotor event-related desynchronization (ERD) and enhanced ipsilateral event-related synchronization (ERS) is superimposed to drive the positive simulation of the active digital twin, generating a multi-channel synthetic EEG time series carrying motor category labels. After screening, a synthetic EEG dataset is formed. The screening process includes: calculating the deviation between the synthetic EEG dataset and the current real EEG data frames of patients in terms of power spectral density in the μ / β band and Frechet feature distance; when the deviation exceeds a preset threshold, shrinking the perturbation variance scale and resampling until the dataset passes the distribution consistency check and is included in the dataset.

7. The brain-computer interface-based neurorehabilitation training method for patients based on digital twins according to claim 1, characterized in that, In step S4, the method for decoding motion intent based on the adaptation weights of the neural signal data frame specifically includes: The confidence score of each decoding output is calculated, and only neural signal data frames with confidence scores exceeding a preset gating threshold and their corresponding decoding labels are included in the fine-tuning sample pool. The samples in the fine-tuning sample pool are used as supervision signals to perform small-batch incremental updates on the individual adaptation weights. At the same time, the individual adaptation weights generated by the virtual pre-training are used as regularization anchors to apply parameter drift penalties, thereby realizing the dynamic evolution of individual adaptation weights with the patient's neural state.

8. The brain-computer interface-based neurorehabilitation training method for patients based on digital twins according to claim 7, characterized in that, The method for realizing the dynamic evolution of individual adaptation weights with the patient's neural state specifically includes: During the patient's formal rehabilitation training, the confidence level of the motor intention decoding results is continuously evaluated. Samples corresponding to decoding results where the component with the highest value in the classification probability vector exceeds a preset threshold are pseudo-labeled. Once the preset threshold is reached, an online fine-tuning of the individual-specific adaptation layer is triggered, so that the decoding network continuously converges to the patient's current neural state.

9. The brain-computer interface-based neurorehabilitation training method for patients based on digital twins according to claim 1, characterized in that, In step S5, the method of feeding back the measured somatosensory cortex response to the active digital twin for error correction specifically includes: The peak amplitude and latency of somatosensory evoked potentials in the corresponding leads of the somatosensory cortex after functional electrical stimulation were extracted as measured response indices. The measured response indices were subtracted from the twin predicted response indices to form a residual vector. Using the residual vector as the observation, a heavy-tailed noise robust Kalman filter mechanism was used to perform posterior correction on the synaptic gain parameters and cortical excitability parameters of the somatosensory cortex sub-model. The corrected parameters were then incorporated into the neural parameter evolution sequence.

10. The brain-computer interface-based neurorehabilitation training method for patients based on digital twins according to claim 1, characterized in that, In step S5, the method for dynamically adjusting the training intensity and frequency of subsequent training based on the neural plasticity window specifically includes: Based on the mean ratio of excitatory to inhibitory synaptic gain estimates for each cortical node in the model parameter estimates, when the rate of change of the excitatory-inhibition ratio is lower than the plasticity activation threshold, the current training intensity and frequency, as well as the previous functional electrical stimulation parameters, remain unchanged; when the rate of change of the excitatory-inhibition ratio exceeds the plasticity activation threshold, the duration of a single training session is increased to a preset threshold and the training interval is shortened to a preset minimum safe interval, while triggering a full parameter space re-optimization of the functional electrical stimulation parameters; when the mean decoding confidence or the rate of prolongation of somatosensory evoked potential latency exceeds a preset fatigue warning threshold, the training intensity is reduced and the rest interval is extended, while a fatigue warning is sent to the clinician.

Images