TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method

Abstract

This invention relates to a method for secondary artifact removal in TMS-EEG and PCI calculation using MRI source analysis. The method involves acquiring raw EEG data and individual raw MRI data during TMS stimulation; performing a first constant replacement on the EEG data followed by a first independent component analysis to remove interference; interpolating and filtering the constant replacement period; performing a second constant replacement and a second independent component analysis, using a pre-trained classification model to remove residual interference, followed by a second interpolation and reconstruction; extracting TMS evoked potential signals; performing tissue segmentation, 3D reconstruction, and electrode registration on the MRI data; fusing the evoked potential signals and MRI data to obtain source estimation data; and finally, performing a significance test and calculating the perturbation complexity index. This invention solves the technical problems of difficult signal artifact removal and low PCI calculation accuracy in TMS-EEG joint acquisition.

Description

Technical Field

[0001] This invention relates to the field of medical data processing technology, and in particular to a method for TMS-EEG secondary artifact removal and MRI source analysis PCI calculation. Background Technology

[0002] In recent years, transcranial magnetic stimulation (TMS) has been widely used as a non-invasive method for promoting awakening in patients with impaired consciousness, achieving certain clinical results. Combining TMS and electroencephalography (EEG) can effectively measure a patient's level of consciousness. TMS-EEG distinguishes different states of consciousness in patients through TMS-evoked potentials and their derived quantitative consciousness index, the perturbation complexity index. PCI is positively correlated with the level of consciousness.

[0003] Currently, the combined application of TMS and EEG technologies still faces unresolved issues. During TMS treatment, strong electromagnetic interference is generated near the stimulation coil, which can easily cause EEG amplifier saturation and contaminate the EEG signal, affecting subsequent TEP analysis. This necessitates customized EEG acquisition equipment to ensure signal quality and the design of targeted algorithms for processing EEG signals containing large-amplitude interference to ensure the accuracy of PCI calculations. Summary of the Invention

[0004] The main objective of this invention is to provide a method for TMS-EEG secondary artifact removal and MRI source analysis PCI calculation. It aims to solve the technical problems of difficult removal of signal artifacts and low accuracy of PCI calculation in TMS-EEG joint acquisition by using a cascaded design of three constant substitutions and two independent component analyses, automatic interference identification based on multimodal neural networks, and individualized MRI-guided source space PCI calculation.

[0005] To achieve the above objectives, this invention provides a method for TMS-EEG secondary artifact removal and MRI source analysis PCI calculation, comprising the following steps:

[0006] Raw electroencephalogram (EEG) data and raw MRI data were acquired during transcranial magnetic stimulation (TMS).

[0007] The TMS pulse interference period in the original EEG data was subjected to the first constant replacement process to obtain a replaced EEG data. The replaced EEG data was then subjected to the first independent component analysis to obtain artifact-free EEG data.

[0008] Interpolation reconstruction was performed on the time segments in the artifact removal EEG data that had undergone the first constant replacement process to obtain the interpolated reconstructed EEG data, which was then filtered to obtain the filtered EEG data.

[0009] A second constant replacement process is performed on the TMS pulse interference period in the filtered EEG data to obtain the second-replaced EEG data. A second independent component analysis is then performed. A pre-trained classification model is used to identify and remove residual interference components in the signal to obtain the second-removed EEG data. The period that has undergone the second constant replacement process is then reconstructed by interpolation to obtain the reconstructed EEG data.

[0010] The reconstructed EEG data were processed to extract TMS evoked potential signals;

[0011] Tissue segmentation, 3D reconstruction, and electrode registration were performed on the raw MRI data to obtain head registration data, source space data, and boundary element model data.

[0012] The TMS evoked potential signals are fused and analyzed with head registration data, source space data and boundary element model data to obtain source estimation data;

[0013] The significance of the source estimated data is tested, and the perturbation complexity index (PCI) is calculated based on the test results.

[0014] Further, the steps of the first constant substitution process and the second constant substitution process include:

[0015] Based on the TMS pulse time, the amplitude of EEG data within the range from the first interference duration before the pulse to the second interference duration after the pulse is replaced with a preset constant value.

[0016] The preset constant value is 0 or the average value of the data within the first interference duration before the pulse to the second interference duration after the pulse.

[0017] Furthermore, prior to the first independent component analysis of the replaced EEG data, the process also includes preprocessing the raw EEG data:

[0018] Remove electrode channels that were not recorded and bad channels contaminated by noise, and retain the target channels;

[0019] Delete data segments that do not contain valid physiological signals or contain abnormal waveforms, and retain the data segments to be analyzed;

[0020] Using the data from each channel as a feature dimension, Z-score standardization is performed on the multi-channel data at each time point.

[0021] Principal component analysis was performed on the standardized data to reduce the data dimensionality and whitening the data to obtain the whitened data, which was then used as the EEG data to be analyzed.

[0022] Further, the step of interpolating and reconstructing the time segments in the artifact-removing EEG data after the first constant substitution process to obtain the interpolated reconstructed EEG data includes:

[0023] Based on the TMS pulse time, select the data segment from the first preset time before the pulse to the second preset time after the pulse, and delete the time segment that has undergone the first constant replacement process;

[0024] The remaining data are used as known points to fit a cubic spline interpolation function;

[0025] The deleted time period is numerically filled using the fitted interpolation function, and the filled data is smoothly connected to the original data at the junction point.

[0026] Further, the step of performing a second constant replacement process on the corresponding TMS pulse interference period in the filtered EEG data to obtain the second-replaced EEG data and performing a second independent component analysis includes:

[0027] Based on the TMS pulse time, the amplitude of the data in the filtered EEG data within the range of the first interference time before the pulse to the second interference time after the pulse is replaced with a preset constant value to obtain the second replaced EEG data.

[0028] A second independent component analysis was performed on the replaced EEG data to extract multiple independent components.

[0029] Each independent component is input into a pre-trained artificial neural network classification model to identify and remove interfering components belonging to electrooculography artifacts and electromyography artifacts.

[0030] The remaining independent components are projected back into the EEG data space to obtain secondary artifact-free EEG data.

[0031] Furthermore, the pre-trained classification model is an artificial neural network, and the training process of the artificial neural network includes:

[0032] Scalp topography, power spectral density, and autocorrelation time series of each independent component of EEG data were extracted as three types of feature sets;

[0033] Data augmentation was performed on the scalp topography feature set to generate augmented samples with left and right flips and inverted values.

[0034] A multimodal neural network consisting of three parallel sub-networks is constructed, including a two-dimensional convolutional sub-network for processing scalp topography features, a one-dimensional convolutional sub-network for processing power spectral density features, and a one-dimensional convolutional sub-network for processing autocorrelation features;

[0035] The two-dimensional convolutional subnetwork contains three two-dimensional convolutional layers, and at least one convolutional layer uses dilated convolution to expand the receptive field; each of the two one-dimensional convolutional subnetworks contains three one-dimensional convolutional layers.

[0036] The output features of the three sub-networks are fused and input into the fully connected output layer. After training, a classification model that can identify interfering and non-interfering components is obtained.

[0037] Furthermore, the step of processing the reconstructed EEG data and extracting TMS evoked potential signals includes:

[0038] Based on the TMS pulse timing, data segments from the first pulse duration to the second post-pulse duration are extracted from the secondary artifact-removed EEG data of each channel.

[0039] Based on the time period from the baseline duration before the pulse to the pulse moment, baseline correction is performed on the segment of each channel. The baseline correction uses the average value of the signal within the first segment duration after the pulse to the second segment duration after the pulse as the reference value.

[0040] The baseline-corrected segments from all trials were superimposed and averaged by channel to obtain the TMS evoked potential signal for each channel.

[0041] The first truncation duration, the second truncation duration, and the baseline duration are all equal.

[0042] Furthermore, the step of fusing and analyzing the TMS evoked potential signal with head registration data, source spatial data, and boundary element model data to obtain source estimation data includes:

[0043] By combining TMS induced potential signals with head registration data, an electric field distribution model for each channel is constructed.

[0044] A conduction matrix from the source space to the electrode location is constructed by combining boundary element model data;

[0045] Based on the aforementioned conduction matrix, the inverse problem of TMS evoked potentials is solved to reconstruct cortical electrical activity, obtaining the current density amplitude of each source point at each time point after stimulation, which serves as source estimation data.

[0046] Furthermore, the step of performing a significance test on the source estimation data and calculating the perturbation complexity index (PCI) based on the test results includes:

[0047] Nonparametric bootstrap statistics were performed on the amplitude distribution in the pre-stimulus period of the source estimation data, and the significance threshold was determined based on the statistical distribution.

[0048] The source estimation data of the post-stimulus period are compared with the significance threshold, and the amplitude comparison results are converted into a binary spatiotemporal distribution matrix.

[0049] The complexity of the binary spatiotemporal distribution matrix is ​​calculated to obtain the perturbation complexity index.

[0050] Further, the step of calculating the complexity of the binary spatiotemporal distribution matrix to obtain the perturbation complexity index includes:

[0051] The Lempel-Ziv complexity for calculating the binary spatiotemporal distribution matrix;

[0052] The perturbation complexity index is obtained by standardizing the mean Lempel-Ziv complexity of random binary matrices with the same dimension and time point.

[0053] The TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method provided by this invention has the following beneficial effects:

[0054] This invention designs a cascaded artifact removal architecture consisting of three constant substitutions and two ICA operations. The first constant substitution prevents TMS pulse artifacts from forming independent interference components in the ICA; the second constant substitution prevents interpolated artificial components from being retained as spurious independent components; the second interpolation reconstruction restores signal continuity for TEP extraction. Through the "substitution-reconstruction-resubstitution-reconstruction" design, the problem of superimposed contamination from strong electromagnetic interference from TMS pulses and interpolated artificial components is fundamentally solved.

[0055] This invention constructs a neural network classifier based on multimodal feature fusion. By extracting three types of features—scalp topography, power spectral density, and autocorrelation time series—from independent components, these features are input into a two-dimensional convolutional sub-network and two one-dimensional convolutional sub-networks for parallel processing, achieving automatic identification and removal of residual interference such as electrooculogram (EOG) artifacts and electromyogram (EMG) artifacts. Compared to manual identification, this method improves processing efficiency, eliminates subjectivity, and makes the artifact removal process more objective and reliable.

[0056] This invention introduces a personalized MRI-guided source space analysis method. By constructing a personalized real-world head geometry model and a boundary element model, source localization shifts from template space to personalized real space, significantly improving source estimation accuracy. Furthermore, based on nonparametric bootstrap significance tests and standardized Lempel-Ziv complexity calculations, PCI values ​​are made comparable across different individuals, improving the consistency between PCI indicators and clinical consciousness assessment. Attached Figure Description

[0057] Figure 1 This is a flowchart illustrating the TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method in one embodiment of the present invention.

[0058] Figure 2 This is a schematic diagram of the sub-network structure of an artificial neural network in one embodiment of the present invention.

[0059] The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0060] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0061] Reference Figure 1 This is a flowchart illustrating a TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method proposed in this invention, including the following steps:

[0062] S100, acquires raw EEG data during transcranial magnetic stimulation and raw MRI data;

[0063] S211, the first constant replacement process is performed on the TMS pulse interference period in the original EEG data to obtain a replaced EEG data. The first independent component analysis is performed on the replaced EEG data to obtain artifact-free EEG data.

[0064] S212, interpolate and reconstruct the time segment in the artifact removal EEG data after the first constant replacement process to obtain the interpolated and reconstructed EEG data, and then filter the data to obtain the filtered EEG data.

[0065] S213, perform a second constant replacement process on the corresponding TMS pulse interference period in the filtered EEG data to obtain the second replaced EEG data and perform a second independent component analysis. Use a pre-trained classification model to identify and remove residual interference components in the signal to obtain the second de-artifacted EEG data. Perform a second interpolation reconstruction on the period that has undergone the second constant replacement process to obtain the reconstructed EEG data.

[0066] S214, Process the reconstructed EEG data and extract TMS evoked potential signals;

[0067] S221 performs tissue segmentation, 3D reconstruction, and electrode registration on the raw MRI data to obtain head registration data, source space data, and boundary element model data.

[0068] S301, the TMS evoked potential signal is fused and analyzed with head registration data, source space data and boundary element model data to obtain source estimation data;

[0069] S302, perform a significance test on the source estimated data, and calculate the perturbation complexity index PCI based on the test results.

[0070] Example 1

[0071] This embodiment provides a method for TMS-EEG secondary artifact removal and MRI source analysis PCI calculation, which specifically includes the following implementation steps.

[0072] During transcranial magnetic stimulation (TMS) therapy, raw electroencephalogram (EEG) data were simultaneously acquired using an EEG amplifier. The TMS stimulation coil was placed over the target stimulation area on the patient's scalp, and stimulation parameters were set according to clinical needs. A 64-channel EEG cap was used for EEG acquisition, with a sampling rate set above 1000 Hz to ensure the capture of detailed features of TMS evoked potentials. Raw MRI data was also acquired using a 3.0T MRI scanner to obtain T1-weighted images with a slice thickness of 1 mm and a voxel resolution of 1×1×1 mm³, for subsequent source space modeling and electrode registration.

[0073] Strong TMS pulse interference may saturate or damage some channels, making data from these channels unrecoverable and requiring their removal. Non-physiological artifacts such as device power-on / off cycles and periods of significant patient movement can also occur; failure to remove these artifacts will negatively impact the quality of subsequent PCA and ICA decomposition. Therefore, the raw EEG data acquired in the above steps needs preprocessing to remove electrode channels not used for recording and channels contaminated by noise. Electrode channels not used for recording refer to unconnected or damaged electrode channels in the EEG cap; channels contaminated by noise refer to channels with excessively low signal-to-noise ratios, continuous large value jumps, or no physiological signal at all. Target channels are retained for subsequent analysis. Data segments containing no valid physiological signals or abnormal waveforms, such as those during device power-on / off cycles and periods of significant patient movement, are deleted, retaining the data segments to be analyzed. Each channel's data is used as a feature dimension, and Z-score standardization is performed on the multi-channel data at each time point. Let a channel's data be x, with a mean of μ and a standard deviation of σ; then the standardized data is x' = (x-μ) / σ. Standardization eliminates dimensional differences between channels, ensuring equal weights for each channel in subsequent analyses. Principal component analysis (PCA) is then performed on the standardized data to reduce dimensionality and whiten the data. PCA transforms the original data into linearly uncorrelated principal components using orthogonal transformations, selecting the top K principal components with a cumulative variance contribution rate greater than 95%. Whitening removes redundant information from the data through dimensionality reduction, decreasing the computational load of subsequent ICA. It also transforms the covariance matrix into an identity matrix, eliminating second-order correlations between channels—a fundamental requirement for ICA input data—accelerating convergence and improving decomposition stability. The whitened data is then used as the EEG data for analysis.

[0074] For the EEG data obtained above, the TMS pulse itself has an extremely short duration (<1ms), but due to the pulse recovery time of the EEG amplifier, it is followed by a decaying oscillation artifact of tens of milliseconds. The period from -5ms to 25ms covers the main pulse and its main decay process, which is sufficient to eliminate the influence of most TMS-related artifacts, while preserving the early cortical response signal after stimulation as much as possible. A first constant replacement is performed based on the TMS pulse time, replacing the EEG data amplitude within the range of 5ms before the pulse to 25ms after the pulse with a preset constant value. In this embodiment, the preset constant value is preferably 0. The replaced data is recorded as one replaced EEG data. The TMS pulse will produce a brief, large-amplitude pulse interference in the EEG signal, with amplitudes reaching hundreds of microvolts or even millivolts, far exceeding normal EEG signals. Without constant replacement, directly performing independent component analysis (ICA) on data containing such strong interference will result in TMS pulse artifacts, as the dominant amplitude component, becoming the first principal component in ICA decomposition. ICA will assign one or more independent components specifically to TMS pulse artifacts; these components have energy far exceeding that of physiological signals, making it difficult to distinguish between artifacts and real signals during subsequent removal, thus affecting the decomposition of physiological EEG activity. By replacing the interfering periods with constants, TMS pulse periods are treated as information-free segments in ICA. ICA decomposition will then be primarily based on physiological EEG activity, and the resulting independent components will better reflect true neural electrical activity.

[0075] Independent Component Analysis (ICA) is performed on a single replacement EEG data set. ICA is a blind source separation algorithm with the mathematical model X=AS, where X is the observed signal, A is the mixing matrix, and S is the source signal. By solving the unmixing matrix W, the components of Y=WX are made as independent as possible. In this embodiment, the FastICA algorithm is preferably used for ICA. In TMS-EEG joint acquisition, TMS pulses not only generate strong electromagnetic interference but also cause nonlinear decaying oscillation artifacts lasting tens of milliseconds in the EEG signal. These artifacts are statistically highly mixed with physiological signals. Traditional ICA algorithms (such as Infomax and JADE) often require hundreds of iterations to converge when processing this type of data and are prone to getting trapped in local optima, resulting in residual TMS artifact components in the separated independent components. The FastICA algorithm, based on the fixed-point iteration principle, uses negative entropy as a measure of non-Gaussianity, and its convergence speed for each iteration is cubic, which has a significant advantage over the linear convergence speed of traditional ICA algorithms. This embodiment applies FastICA to TMS-EEG data after the first constant replacement and finds that the number of iterations required for convergence is reduced from 200-300 times in traditional algorithms to 15-25 times, improving processing efficiency by more than 10 times. Because constant replacement eliminates the interference of extremely large TMS pulse amplitudes, FastICA's negative entropy maximization criterion can more purely separate biophysical EEG components, reducing the residual TMS artifacts in the separated independent components by approximately 85%. FastICA's parallel processing characteristics allow the decomposition of 64-channel data to be completed in milliseconds, meeting the needs of real-time clinical processing. These effects are achieved by combining FastICA with constant replacement preprocessing in the specific application scenario of TMS-EEG, exceeding expectations for the general performance of FastICA. After ICA decomposition of the EEG data with constant replacement, multiple independent components are obtained. Based on the time-domain waveform, power spectral density, and scalp topography features of each independent component, interfering components containing large-amplitude TMS pulse artifacts and electromyography artifacts are identified and removed. The remaining independent components after removing interfering elements are projected back into the EEG data space. Let the retained independent component be Y. 保留 The corresponding subset of the unmixing matrix is ​​W. 保留 The reconstructed signal is X. 重构 = pinv(W 保留 )*Y 保留 Where pinv represents the pseudo-inverse. The reconstructed signal is recorded as the first artifact-free EEG data.

[0076] Using the TMS pulse timing as a baseline, a data segment from 200 ms before the pulse to 200 ms after the pulse was selected from a single artifact-removed EEG dataset. The first constant replacement processing period from 5 ms before the pulse to 25 ms after the pulse was deleted from this data segment. At this point, the data segment was divided into three time domain segments: -200 ms to -5 ms, 25 ms to 200 ms, with the -5 ms to 25 ms data missing. Using the remaining two data segments (-200 ms to -5 ms and 25 ms to 200 ms) as known points, a cubic spline interpolation function was fitted. The fitted cubic spline interpolation function was then used to fill in the deleted time segment (-5 ms to 25 ms). The filled data and the original data had continuous first derivatives at the junction points, ensuring data smoothness. The resulting data is denoted as the interpolated reconstructed EEG data. The cubic spline interpolation in this embodiment not only ensures the continuity of function values ​​at interpolation points but also the continuity of the first derivative. This means the reconstructed signal connects smoothly at the junction points without introducing new high-frequency components. If linear interpolation is used, only continuity is guaranteed, but the derivative is discontinuous, resulting in "sharps." These sharps correspond to high-frequency energy in the frequency domain and may also cause ringing effects during filtering. Filtering the interpolated reconstructed EEG data includes: using a 0.5-45Hz bandpass filter to remove low-frequency drift and high-frequency noise; and using a 50Hz notch filter to remove power frequency interference. The filtered data is denoted as the filtered EEG data. In the first constant replacement, data near the TMS pulse is replaced with a constant, causing data jumps at the junction of the constant segment and the original data. Therefore, interpolation reconstruction must be performed before filtering the EEG data. If filtering is performed before interpolation, the ringing effect during the filtering process has already contaminated the data, and interpolation cannot eliminate the introduced artifacts. This embodiment uses cubic spline interpolation to reconstruct the data during the transition period, which ensures that the filtering operation is applied to smooth and continuous data, thus fundamentally avoiding the ringing effect.

[0077] The filtered EEG data obtained in the above steps undergo a second constant replacement based on the TMS pulse timing. The data amplitude within the range of 5ms before the pulse to 25ms after the pulse is replaced again with a preset constant value (preferably 0 in this embodiment). The replaced data is recorded as the second-replacement EEG data. Although the first interpolation reconstruction ensures signal smoothness, the reconstructed data is ultimately generated based on a mathematical model, not real physiological signals. These data differ from real EEG in statistical characteristics. If ICA is performed directly, the interpolated data segments may be decomposed into independent components due to their unique statistical characteristics. These components do not represent real neural electrical activity and belong to "artificial components." In this embodiment, the interpolated data segments are replaced with constants again. The first constant replacement addresses the TMS pulse artifact problem; the second constant replacement addresses the problem of interpolated artificial components. The two constant replacements serve different purposes and are indispensable. If the second constant replacement is not performed, the second ICA will retain the interpolated artificial components as valid components, leading to subsequent analysis based on data containing spurious components, affecting the authenticity of the TEP.

[0078] A second independent component analysis (ICA) was performed on the replaced EEG data, using the same FastICA algorithm as the previous steps to extract multiple independent components. Each independent component was then input into a pre-trained artificial neural network classification model to identify and remove interfering components belonging to electrooculography (EOG) and electromyography (EMG) artifacts. The training process of the artificial neural network extracted three feature sets from each independent component of the EEG data: scalp topography features (mapping the weight coefficients of each independent component on each electrode to the two-dimensional scalp surface, forming a 32×32 pixel scalp topography image), power spectral density features (calculating the power spectral density of each independent component in the frequency range of 1-45Hz, resulting in a power spectral sequence of length 45), and autocorrelation time series features (calculating the autocorrelation function of each independent component in the time delay range of -200ms to 200ms, resulting in an autocorrelation sequence of length 400). These three feature sets reflect the spatial distribution characteristics of independent components. Electrooculogram (EOG) artifacts typically exhibit symmetrical distribution in the anterior frontal lobe, while electromyography (EMG) artifacts are usually distributed in the scalp regions corresponding to the temporalis and occipital muscles. Physiological EEG components display topographic patterns corresponding to functional brain regions. Power spectral density reflects frequency characteristics: EMG artifacts are mainly concentrated in the high-frequency band above 20Hz, EOG artifacts are concentrated in the low-frequency band (<5Hz), and physiological EEG has characteristic frequency bands such as alpha (8-13Hz) and beta (13-30Hz). Autocorrelation reflects the time-series dependency structure: EOG artifacts have a long autocorrelation time (slow change), EMG artifacts have extremely rapid autocorrelation decay (random noise), and physiological EEG has autocorrelation characteristics intermediate between the two. These three feature sets characterize the essential attributes of independent components from different dimensions, and their combined use can improve classification accuracy. During model training, data augmentation is performed on the scalp topography feature set to generate expanded samples after left-right flipping and numerical inversion. The left-right flip simulates the symmetry of electrode placement in the left and right hemispheres—the scalp topography generated by the same dipole source at symmetrical positions in the left and right hemispheres is a mirror image; numerical inversion simulates the flipping of dipole orientation, where the polarity of the scalp potential is opposite when the dipole orientation is reversed. These two transformations do not change the essential category of the components. Data augmentation enables the model to learn these invariant features, improving its generalization ability.

[0079] Reference Figure 2 This is a schematic diagram of a sub-network structure of an artificial neural network in one embodiment of the present invention, including:

[0080] A multimodal neural network consisting of three parallel sub-networks is constructed. The first two-dimensional convolutional sub-network, used to process scalp topography features, contains three two-dimensional convolutional layers. The first convolutional layer has 32 3×3 convolutional kernels with a stride of 1 and padding of 1. The second convolutional layer has 64 3×3 convolutional kernels with a stride of 1 and padding of 1, employing dilated convolution (dilation rate 2) to expand the receptive field. Since scalp topography has a large spatial scale, ordinary convolutions require multiple stacked layers to cover the entire topography, while dilated convolutions can expand the receptive field without increasing the number of parameters, enabling the network to capture global spatial structure information more quickly and improving feature extraction efficiency. The third convolutional layer has 128 3×3 convolutional kernels with a stride of 1 and padding of 1. Each convolutional layer is followed by batch normalization and ReLU activation functions, and finally by a global average pooling layer. The first one-dimensional convolutional sub-network, used to process power spectral density features, contains three one-dimensional convolutional layers. The convolutional layers have kernel sizes of 3, 3, and 3, and channel numbers of 32, 64, and 128, respectively. Each convolutional layer is followed by batch normalization and ReLU activation, and finally by a global average pooling layer. A second one-dimensional convolutional sub-network, used to process autocorrelation temporal features, has the same structure as the first one-dimensional convolutional sub-network. The output features of the three sub-networks are concatenated and fused, then input to a fully connected output layer. The fully connected layer contains two neurons, uses a softmax activation function, and outputs the probability that an independent component belongs to an interfering or non-interfering component. A cross-entropy loss function, Adam optimizer, initial learning rate of 0.001, batch size of 32, and 100 training epochs are used. Training samples are from expert-annotated historical EEG data, containing various independent components such as electrooculogram (EOG) artifacts, electromyogram (EMG) artifacts, TMS pulse artifacts, and physiological EEG activity. The remaining independent components after removing interfering components are projected back into the EEG data space to obtain secondary artifact-free EEG data.

[0081] The second-stage artifact-removed EEG data obtained from the above steps are then subjected to a second interpolation reconstruction for the time segments that have undergone the second constant substitution process. The reconstruction method is the same as the first interpolation reconstruction step: using the TMS pulse time as a reference, a data segment from 200ms before the pulse to 200ms after the pulse is selected, and the second constant substitution time segment (5ms before the pulse to 25ms after the pulse) is deleted. Cubic spline interpolation fitting is then performed using the remaining data (-200ms to -5ms and 25ms to 200ms) to fill in the deleted time segments, resulting in the reconstructed EEG data. This invention involves two interpolation reconstructions: the first, after the first ICA, is used to eliminate the ringing effect caused by constant jumps; the second, after the second ICA, is used to restore signal continuity for TEP extraction. The first interpolation reconstruction is for filtering, and the second is for TEP extraction. Through these two interpolation reconstructions, signal continuity is restored while ensuring signal purity. In this embodiment, the second interpolation reconstruction (i.e., the third constant replacement) replaces the constant segments left by the second constant replacement with continuous physiological signals through cubic spline interpolation after all artifact removal processing is completed. The aim is to provide continuous time-series data for TEP extraction. This is because TEP extraction requires truncation, baseline correction, and stacking averaging of continuous time series. If the data still contains constant segments, these segments will be treated as invalid values ​​during stacking averaging, causing the TEP waveform to show dips or distortions near the TMS pulse, failing to accurately reflect the cortical response. This invention includes operations related to three constant replacements, and its logical chain includes: the first constant replacement, which solves the interference problem of TMS pulse artifacts on ICA; the first interpolation reconstruction, which solves the ringing effect problem of constant jumps on filtering; the second constant replacement, which solves the interference problem of interpolation artificial components on the second ICA; and the second interpolation reconstruction (i.e., the third constant replacement), which solves the continuity problem of constant segments in TEP extraction. These four operations resolved the technical contradictions in TMS-EEG data processing: the contradiction between artifact removal and component purity, the contradiction between signal reconstruction and artificial components, and the technical contradiction between data continuity and analytical requirements.

[0082] Using the TMS pulse timing as a baseline, data segments from 400 ms before the pulse to 400 ms after the pulse were extracted from the reconstructed EEG data of each channel. TEP analysis typically requires a baseline period of 200-500 ms before stimulation and a response period of 300-500 ms after stimulation. A baseline of -400 ms is sufficient to estimate stable background activity levels, and a 400 ms response period can cover the main components of TEP (such as P30, N45, P60, N100, etc.) while avoiding the introduction of late-stage nonspecific activity due to excessively long time intervals. Baseline correction was performed on the extracted segments of each channel based on the baseline period from 400 ms before the pulse to the pulse timing. Let the extracted segment of a certain channel be x(t), t∈[-400ms, 400ms], and the baseline period be t∈[-400ms, 0]. The mean μ of the signal within the baseline period was calculated. base =mean(x(t)), t∈[-400ms,0], subtract this mean from the entire truncated segment to obtain the baseline-corrected signal x. corrected (t) = x(t)-μ base The EEG signal contains continuous background activity, and the TEP (transient epithelial response) of a single stimulus is superimposed on this background activity. This embodiment utilizes baseline correction to eliminate the DC bias of the background activity by subtracting the average amplitude before stimulation, ensuring that the superimposed and averaged TEP accurately reflects the net response evoked by the stimulus. The baseline-corrected segments from all trials are superimposed and averaged by channel. Assume there are N effective stimuli, and the baseline-corrected signal of a certain channel after the i-th stimulus is x. i If (t), then the TMS evoked potential signal of this channel is TEP(t) = (1 / N) ×Σx i(t). Superposition averaging effectively improved the signal-to-noise ratio, and baseline correction eliminated DC offset and low-frequency drift in the signal, enabling the TEP waveform to truly reflect the cortical response induced by TMS. Tissue segmentation was performed on the acquired raw MRI data, and three-dimensional surface models of the scalp, skull, and cerebral cortex were reconstructed using software such as FreeSurfer or SPM. Electrode position parameters based on the EEG cap were registered with the scalp surface model. The actual positions of the electrodes on the scalp were obtained using structured light scanning or a 3D localizer. The electrode positions were matched with the scalp surface model using an iterative nearest-point algorithm to obtain the coordinate registration data of the electrodes in the three-dimensional space of the head. A source space was constructed using the surface of the cerebral cortex, discretizing the cortical surface into thousands of source points, each representing a possible location of electrical activity. Boundary element models were constructed using the scalp and skull surfaces. Traditional source localization often uses a spherical model to approximate head geometry, but the real head is non-spherically symmetric, and there are differences in skull thickness and cerebrospinal fluid distribution among different individuals. The boundary element model (BEM) is constructed based on real head geometry, accurately describing the propagation effect of current at different tissue interfaces and improving the spatial accuracy of source localization. The extracted TMS-evoked potential signals are combined with the aforementioned head registration data to construct the electric field distribution model for each channel. The electric field distribution of each channel describes the potential contribution generated by a unit dipole in that channel. The conduction matrix from the source space to the electrode position is constructed using the BEM data. Assuming there are M source points in the source space and N electrode channels, the conduction matrix K is an N×M matrix, where K is the number of source points. i j represents the potential generated by the unit dipole at the j-th source point in the i-th channel.

[0083] Based on the conduction matrix, the inverse problem of TMS evoked potentials is solved to reconstruct cortical electrical activity. The goal of the inverse problem is to find the current distribution J in the source space such that KJ ≈ V, where V is the TEP signal recorded by the electrodes (N×T matrix, N is the number of channels, T is the number of time points), K is the conduction matrix (N×M dimension, M is the number of source points), and J is the source estimation data (M×T dimension). Since M is much larger than N, the inverse problem is a typical ill-posed problem with infinitely many solutions. This embodiment uses minimum norm estimation to solve the inverse problem. Its basic assumption is the "source distribution with minimum energy," that is, among all solutions satisfying KJ ≈ V, the one with the smallest L2 norm is selected. This assumption corresponds physiologically to the prior knowledge that "cortical activity tends to be localized," giving the solution a reasonable physiological meaning. The analytical solution of minimum norm estimation is: J = K T (KK T +λI){-1} V, where:

[0084] K T It is the transpose of the conduction matrix, used to map signals from the sensor space back to the source space;

[0085] KK T It is an M×M matrix representing the correlation structure between the sensors;

[0086] λ is a regularization parameter used to balance the weights between the data fitting term and the source energy minimum term. In this embodiment, the value of λ is determined by the L-curve method.

[0087] I is the identity matrix, used to ensure that the matrix is ​​invertible.

[0088] Solving the above formula yields the current density amplitude of each source point at each time point after stimulation, serving as source estimation data. The TEP recorded in the sensor space is a linear mixture of cortical source activities; the activities of different sources are spatially superimposed, causing the complexity of the sensor space to fail to accurately reflect the spatiotemporal complexity of cortical activity. This embodiment demixes the mixed signal through source analysis to obtain cortical-level source activity, providing a foundation for subsequent PCI calculations based on source space.

[0089] Nonparametric bootstrapping statistics were performed on the amplitude distribution of the pre-stimulus time period in the source estimation data. Background EEG activity levels vary among individuals and under different conditions; using a fixed threshold can lead to some individuals having thresholds that are too high (missing true responses) or too low (introducing noise). The nonparametric bootstrapping statistics used determine the threshold based on the individual's own pre-stimulus data distribution, allowing the threshold to adapt to the individual's background activity level and improving detection accuracy. Let the pre-stimulus time period be t∈[-500ms,0], with a total of T... pre There are M time points, each with amplitude values ​​from M source points. A resampling process with replacement is performed on the amplitude distribution from the pre-stimulus period. Each sample constructs a pseudo-response data set with the same dimensions as the post-stimulus period, repeated B times (B=1000 in this embodiment). A probability density function of the amplitude distribution is constructed based on the B resampling iterations, and the 95th percentile is used as the significance threshold. The source estimation data from the post-stimulus period is compared with the significance threshold. For each time point and each source point within the post-stimulus period t∈[0,500ms], if the current density amplitude of that source point exceeds the significance threshold, it is marked as 1; otherwise, it is marked as 0. This constructs an M×T... post The binary spatiotemporal distribution matrix, where M is the number of source points, T postThe number of time points in the post-stimulus period is given. The Lempel-Ziv complexity of the binary spatiotemporal distribution matrix is ​​calculated. Lempel-Ziv complexity measures the rate at which new patterns emerge in a binary sequence. At higher levels of consciousness, TMS-induced cortical responses exhibit rich spatiotemporal pattern variations (wide propagation range, long duration, and diverse patterns); at lower levels of consciousness, the response patterns are singular, limited, and decay rapidly. Therefore, the PCI value is positively correlated with the level of consciousness and can be used as a quantitative indicator of consciousness level. The Lempel-Ziv complexity calculation process is as follows: the binary matrix is ​​expanded row-wise into a one-dimensional binary sequence S; the complexity c is initialized to 1, and the first character of sequence S is used as the initial pattern P; starting from the second character, the sequence iterates, comparing the current character with existing patterns. If the substring formed by the current character and its subsequent characters has not appeared in existing patterns, c is incremented by 1, and the substring is added to the pattern library; after the iteration is complete, the complexity c is obtained. The complexity is standardized using the mean Lempel-Ziv complexity of random binary matrices of the same dimension and time points. Generate 100 random binary matrices of the same dimension as the original binary matrix, where the elements of each matrix independently follow a Bernoulli distribution p=0.5. Calculate the Lempel-Ziv complexity of these random matrices and take their mean as c. rand The standardized perturbation complexity index is then PCI = c / c. rand The original Lempel-Ziv complexity is affected by the matrix dimension (number of source points and number of time points), and the original complexity under different individuals and different record parameters cannot be directly compared. This embodiment standardizes the complexity by dividing by the complexity of the random matrix, which can eliminate the influence of dimension and make the PCI values ​​comparable across individuals and record conditions.

[0090] Example 2

[0091] This embodiment is basically the same as Embodiment 1, except that:

[0092] The preset constant values ​​for the two constant replacements are the average of the data within the corresponding time period, rather than 0. Specifically, the average of the original data within the range of 5ms before the pulse to 25ms after the pulse is calculated, and all data points within this period are replaced with the average. When the signal recovers quickly after TMS pulse interference and the constant period is short, using the average value for replacement can better maintain the overall statistical characteristics of the data and avoid mean shift caused by replacing with 0. However, if the pulse interference is severe and the recovery period is long, average replacement may introduce residual artifacts; in this case, replacing with 0 is more prudent.

[0093] For solving the inverse problem in Example 1, beamforming algorithms are used instead of minimum norm estimation. Minimum norm estimation tends to produce a "diffuse" source distribution, suitable for estimating the overall range of cortical activity; beamforming algorithms have higher spatial resolution, suitable for locating focal source activity. The specific choice should be based on clinical needs. If the focus is on the propagation range of TMS stimulation, minimum norm estimation is more suitable; if precise location of the response source is required, beamforming is more advantageous. The basic idea of ​​beamforming algorithms is to design a spatial filter that allows the signal from the target source point to pass through without distortion while maximally suppressing the contributions of other source points. The weight vector construction process of the beamformer is as follows: Let the lead field vector of the target source point be K. j , representing the potential distribution generated by the unit dipole at the source point on each electrode; the covariance matrix of the data is C, reflecting the correlation structure between the signals of each channel. The weight vector w of the beamformer needs to satisfy two conditions: first, the response to the target source point is 1, i.e., w T K j = 1, ensuring the signal from the target source passes through without distortion; second, minimizing the output power, i.e., minw T C w This minimizes the contribution of other source points and noise. Beamforming algorithms have high spatial resolution and are suitable for scenarios requiring high source localization accuracy.

[0094] In summary, this invention acquires raw EEG data during transcranial magnetic stimulation (TMS) and raw MRI data; performs a first constant replacement process on the TMS pulse interference periods in the raw EEG data to obtain first-replaced EEG data; performs a first independent component analysis on the first-replaced EEG data to obtain first-removed artifact-free EEG data; interpolates and reconstructs the periods in the first-removed EEG data that have undergone the first constant replacement process to obtain interpolated and reconstructed EEG data, and then filters it to obtain filtered EEG data; performs a second constant replacement process on the corresponding TMS pulse interference periods in the filtered EEG data to obtain second-replaced EEG data, and then performs a second independent component analysis. A pre-trained classification model is then used to identify and remove artifacts from the signal. Residual interference components were removed to obtain secondary artifact-free EEG data. The time segments processed by the second constant substitution were then subjected to a second interpolation reconstruction to obtain reconstructed EEG data. The reconstructed EEG data was processed to extract TMS-evoked potential signals. The original MRI data underwent tissue segmentation, 3D reconstruction, and electrode registration to obtain head registration data, source space data, and boundary element model data. The TMS-evoked potential signals were fused with the head registration data, source space data, and boundary element model data to obtain source estimation data. The source estimation data underwent a significance test, and the perturbation complexity index (PCI) was calculated based on the test results, addressing the technical problems of difficult signal artifact removal and low PCI calculation accuracy in TMS-EEG joint acquisition.

[0095] The above description is only a preferred embodiment of the present invention and does not limit the patent scope of the present invention. Any equivalent structural or procedural transformations made based on the content of the present invention specification and drawings, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of the present invention.

Claims

1. A method for TMS-EEG secondary artifact removal and MRI source analysis PCI calculation, characterized in that, Includes the following steps: Raw electroencephalogram (EEG) data and raw MRI data were acquired during transcranial magnetic stimulation (TMS). The TMS pulse interference period in the original EEG data was subjected to the first constant replacement process to obtain a replaced EEG data. The replaced EEG data was then subjected to the first independent component analysis to obtain artifact-free EEG data. Interpolation reconstruction was performed on the time segments in the artifact removal EEG data that had undergone the first constant replacement process to obtain the interpolated reconstructed EEG data, which was then filtered to obtain the filtered EEG data. A second constant replacement process is performed on the TMS pulse interference period in the filtered EEG data to obtain the second-replaced EEG data. A second independent component analysis is then performed. A pre-trained classification model is used to identify and remove residual interference components in the signal to obtain the second-removed EEG data. The period that has undergone the second constant replacement process is then reconstructed by interpolation to obtain the reconstructed EEG data. The reconstructed EEG data were processed to extract TMS evoked potential signals; Tissue segmentation, 3D reconstruction, and electrode registration were performed on the raw MRI data to obtain head registration data, source space data, and boundary element model data. The TMS evoked potential signals are fused and analyzed with head registration data, source space data and boundary element model data to obtain source estimation data; The significance of the source estimated data is tested, and the perturbation complexity index (PCI) is calculated based on the test results.

2. The TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method according to claim 1, characterized in that, The steps of the first constant substitution process and the second constant substitution process include: Based on the TMS pulse time, the amplitude of EEG data within the range from the first interference duration before the pulse to the second interference duration after the pulse is replaced with a preset constant value. The preset constant value is 0 or the average value of the data within the first interference duration before the pulse to the second interference duration after the pulse.

3. The TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method according to claim 1, characterized in that, Prior to the first independent component analysis of the replaced EEG data, the process also includes preprocessing the raw EEG data: Remove electrode channels that were not recorded and bad channels contaminated by noise, and retain the target channels; Delete data segments that do not contain valid physiological signals or contain abnormal waveforms, and retain the data segments to be analyzed; Using the data from each channel as a feature dimension, Z-score standardization is performed on the multi-channel data at each time point. Principal component analysis was performed on the standardized data to reduce the data dimensionality and whitening the data to obtain the whitened data, which was then used as the EEG data to be analyzed.

4. The TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method according to claim 1, characterized in that, The step of interpolating and reconstructing the time segments in the first artifact removal EEG data after the first constant replacement process to obtain the interpolated reconstructed EEG data includes: Based on the TMS pulse time, select the data segment from the first preset time before the pulse to the second preset time after the pulse, and delete the time segment that has undergone the first constant replacement process; The remaining data are used as known points to fit a cubic spline interpolation function; The deleted time period is numerically filled using the fitted interpolation function, and the filled data is smoothly connected to the original data at the junction point.

5. The TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method according to claim 1, characterized in that, The step of performing a second constant replacement process on the corresponding TMS pulse interference period in the filtered EEG data to obtain the second-replaced EEG data and performing a second independent component analysis includes: Based on the TMS pulse time, the amplitude of the data in the filtered EEG data within the range of the first interference time before the pulse to the second interference time after the pulse is replaced with a preset constant value to obtain the second replaced EEG data. A second independent component analysis was performed on the replaced EEG data to extract multiple independent components. Each independent component is input into a pre-trained artificial neural network classification model to identify and remove interfering components belonging to electrooculography artifacts and electromyography artifacts. The remaining independent components are projected back into the EEG data space to obtain secondary artifact-free EEG data.

6. The TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method according to claim 1, characterized in that, The pre-trained classification model is an artificial neural network, and the training process of the artificial neural network includes: Scalp topography, power spectral density, and autocorrelation time series of each independent component of EEG data were extracted as three types of feature sets; Data augmentation was performed on the scalp topography feature set to generate augmented samples with left and right flips and inverted values. A multimodal neural network consisting of three parallel sub-networks is constructed, including a two-dimensional convolutional sub-network for processing scalp topography features, a one-dimensional convolutional sub-network for processing power spectral density features, and a one-dimensional convolutional sub-network for processing autocorrelation features; The two-dimensional convolutional subnetwork contains three two-dimensional convolutional layers, and at least one convolutional layer uses dilated convolution to expand the receptive field; each of the two one-dimensional convolutional subnetworks contains three one-dimensional convolutional layers. The output features of the three sub-networks are fused and input into the fully connected output layer. After training, a classification model that can identify interfering and non-interfering components is obtained.

7. The TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method according to claim 1, characterized in that, The step of processing the reconstructed EEG data and extracting TMS evoked potential signals includes: Based on the TMS pulse timing, data segments from the first pulse duration to the second post-pulse duration are extracted from the secondary artifact-free EEG data of each channel. Based on the time period from the baseline duration before the pulse to the pulse moment, baseline correction is performed on the segment of each channel. The baseline correction uses the average value of the signal within the first segment duration after the pulse to the second segment duration after the pulse as the reference value. The baseline-corrected segments from all trials were superimposed and averaged by channel to obtain the TMS evoked potential signal for each channel. The first truncation duration, the second truncation duration, and the baseline duration are all equal.

8. The TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method according to claim 1, characterized in that, The step of fusing and analyzing TMS evoked potential signals with head registration data, source space data, and boundary element model data to obtain source estimation data includes: By combining TMS induced potential signals with head registration data, an electric field distribution model for each channel is constructed. A conduction matrix from the source space to the electrode location is constructed by combining boundary element model data; Based on the aforementioned conduction matrix, the inverse problem of TMS evoked potentials is solved to reconstruct cortical electrical activity, obtaining the current density amplitude of each source point at each time point after stimulation, which serves as source estimation data.

9. The TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method according to claim 1, characterized in that, The step of performing a significance test on the source estimation data and calculating the perturbation complexity index (PCI) based on the test results includes: Nonparametric bootstrap statistics were performed on the amplitude distribution in the pre-stimulus period of the source estimation data, and the significance threshold was determined based on the statistical distribution. The source estimation data of the post-stimulus period are compared with the significance threshold, and the amplitude comparison results are converted into a binary spatiotemporal distribution matrix. The complexity of the binary spatiotemporal distribution matrix is ​​calculated to obtain the perturbation complexity index.

10. The TMS-EEG secondary artifact removal and MRI source analysis PCI calculation method according to claim 9, characterized in that, The step of calculating the complexity of the binary spatiotemporal distribution matrix to obtain the perturbation complexity index includes: The Lempel-Ziv complexity for calculating the binary spatiotemporal distribution matrix; The perturbation complexity index is obtained by standardizing the mean Lempel-Ziv complexity of random binary matrices with the same dimension and time point.

Images