A magnetoencephalogram tracing method based on deep neural network

Abstract

The present invention discloses a magnetoencephalogram (MEG) source tracing method based on a deep neural network, which relates to the field of brain source imaging. MEG data are obtained through a pre-set head model, and the MEG data are traced at a single moment through a neural network framework LassoNet to obtain a source dipole vector corresponding to each moment. A plurality of source dipole vectors are combined in chronological order to obtain a source signal sequence, and the source signal sequence is input into a trained source tracing model, and a reconstructed source tracing result is output. Based on the reconstructed source tracing result and the head model, a source dipole matrix corresponding to the MEG signal is obtained, and the source position and source strength are calculated according to the source dipole matrix. The MEG source tracing method obtains a more accurate reconstructed source tracing result.

Description

Technical Field

[0001] The present invention relates to the field of brain source imaging, and in particular to a magnetoencephalogram (MEG) source tracing method based on a deep neural network. Background Art

[0002] Magnetoencephalography (MEG) and electroencephalography (EEG) have a temporal resolution of milliseconds, but the number of sensors used in existing MEG and EEG measurement systems is far smaller than the number of equivalent source electric dipoles in the brain. Therefore, traditional algorithms usually artificially add constraints to obtain a unique source position. These constraints limit the applicable scenarios of traditional algorithms, and the robustness of traditional algorithms is usually poor.

[0003] The problem of magnetoencephalography-based tracing can usually be defined as:

[0004] ;

[0005] in, For the The measurement signal of magnetoencephalography at each moment, represents the number of MEG channels, represents the real vector space; is the guidance field matrix, represents the number of source space partitions, 3 represents the three directions of the source dipole, and the guidance field matrix reflects the mapping relationship between the magnetoencephalogram or electroencephalogram sensor and the electric dipole moment of the source at any position; represents a dipole in source space; Representative The noise at the moment. The number of sources in the above problem ( ) is usually much larger than the number of sensors ( ), so solving the source is usually solving an ill-posed problem, and it is difficult to solve the unique source.

[0006] In the norm solution method used for magnetoencephalography tracing, the solution process is defined as:

[0007] ;

[0008] in, express The magnetic brain signals measured at all times; is the guiding field matrix; for The dipole vector at the moment is the solution target; the above formula is a typical traditional method using the norm to solve. Adding a norm as a regularization term can obtain a sparse solution. and Indicates different moments, Represents different regularization parameters, which need to be determined manually.

[0009] Traditional source tracing methods have developed rapidly in the past few decades. Dipole fitting was a widely used method for source tracing in the early days. In 1992, Mosher et al. used the MUSIC (Multiple Signal Classification) algorithm based on subspace decomposition for source tracing. With the development of functional magnetic resonance imaging technology, it has become relatively easy to obtain the structure of the human brain, the technology of generating guidance fields based on forward models has become mature, and the related technology of solving the source using paradigm optimization has developed rapidly. The minimum norm estimation method has been widely used. The number of related methods based on the minimum norm estimation method has gradually increased, such as dSPM (dynamic Statistic Parametric Mapping) using Bayesian maximum a posteriori estimation, sLORETA (standardized low resolution brain electromagnetic tomography) based on standardized current density images, and wMNE (weighted Minimum Norm Estimation) that can estimate sources at different depths.

[0010] In recent years, the development of artificial intelligence technology has brought new development directions to source tracing algorithms. In 2020, Alexandra et al. (Deep Learning for Non-invasive Cortical Potential Imaging) used U-net to reconstruct the source space; in 2022, Sun et al. (Deep neural networks constrained by neural mass models improve electrophysiological source imaging of spatiotemporal brain dynamics) used feedforward neural networks and recurrent neural networks for EEG source tracing; in 2023, Jiao et al. (MMDF-ESI: Multi-Modal Deep Fusion ofEEG and MEG for Brain Source Imaging) fused the features of MEG and EEG and used convolutional neural networks for source tracing.

[0011] Traditional algorithms are mainly used for single-source tracing, and are not very effective for multi-source tracing. Existing deep learning-based algorithms have solved some of the problems, but the following problems still exist:

[0012] The differences in head models between different groups of people are not fully considered, and the same head model is used for model training, resulting in poor model robustness and a large gap in the recall rates of different examples. The source solution process mainly uses convolutional neural networks, recurrent neural networks or feedforward neural networks, etc. These models are used for source tracing calculations and usually require a large number of parameters to retain features, retaining many non-source interference features.

[0013] Although using long short-term memory artificial neural networks (LSTM) or recurrent neural networks (RNN) to learn signal timing features can solve the problem of result jumps at different times in traditional algorithms, RNN and LSTM have weak processing capabilities for global dependencies. Summary of the invention

[0014] Based on the technical problems existing in the background technology, the present invention proposes a magnetoencephalogram tracing method based on deep neural network to obtain more accurate reconstruction tracing results.

[0015] The present invention proposes a magnetoencephalogram (MEG) tracing method based on a deep neural network. The method obtains MEG data through a pre-set head model, performs single-moment tracing on the MEG data through the neural network framework LassoNet, obtains the source dipole vector corresponding to each moment, combines multiple source dipole vectors in chronological order to obtain a source signal sequence, inputs the source signal sequence into a trained tracing model, outputs a reconstructed tracing result, obtains the source dipole matrix corresponding to the MEG signal based on the reconstructed tracing result and the head model, and calculates the source position and source strength according to the source dipole matrix.

[0016] The training process of the traceability model is as follows:

[0017] Step 1: Establish different head models according to the head structure magnetic resonance imaging of people of different age groups, use the head models to establish the guidance field matrix required for tracing, calculate the simulation source matrix based on the improved neural group model, and multiply the guidance field matrix with the simulation source matrix to obtain the simulated brain magnetic signal;

[0018] Step 2: Use the neural network framework LassoNet to trace the source of each moment of the simulated EEG signal and output the single-moment source dipole vector;

[0019] Step 3: Input the single-time source dipole vector and the simulated source matrix into the Transformer architecture. The encoder module in the Transformer architecture encodes the single-time source dipole vector in the time dimension and the source space dimension in stages. The source dipole strength is calculated based on the single-time source dipole vector and time encoding is performed. The time encoding result is sequentially passed through the time self-attention module and the time convolutional network to obtain the time feature matrix. The time feature matrix and the source space dimension encoding result are feature fused. The fused feature matrix is ​​sequentially passed through the self-attention module and the feedforward neural network. Then the output of the encoder module is used as the key K and value V of the multi-head attention module in the decoder module.

[0020] Step 4: The decoder module in the Transformer architecture performs time-series coding on the analog source matrix, and performs source sampling on the analog source matrix after time-series coding. The sampled analog source matrix is ​​respectively calculated and queried. ,key ,value , and after the inverse fast Fourier transform module, the query ,key ,value Transform to obtain query Q;

[0021] Step 5. Based on the key K, value V output from step 3 and query Q output from step 4, calculate the multi-head attention feature, input the feature map output by the multi-head attention module into the pre-trained expert system to reconstruct the source signal, output the decoding attention feature, use the residual connection to fuse the multi-head attention feature and the decoding attention feature, and then perform layer normalization. The normalized result continues to be input into the next decoder module. If it reaches the last decoder module, the normalized result of the last layer is linearly transformed to obtain the reconstructed tracing result.

[0022] Furthermore, in step 1, different head models are established for head structure magnetic resonance imaging of people of different age groups, and the process is as follows:

[0023] According to the age distribution, the brain modeling age groups were divided into: 0 to 3 years old, 3 to 14.4 years old, 14.4 to 28.7 years old, 28.7 to 40 years old, 40 to 55 years old, and 55 to 80 years old;

[0024] According to age groups, different cerebral cortical source models are established, and the number of regional divisions and the location of the regional center of each cerebral cortical source model are determined;

[0025] The cerebral cortex is divided into different areas according to the head model, and the electrical signals fed back from the cerebral cortex are converted into the average pulse density of the emitted action potential to obtain the simulated source signals of each partition. The simulated source signals of each partition are spliced ​​according to the timing process of the source activity to obtain the simulated source matrix.

[0026] Furthermore, in step 2, the optimization goal of the neural network framework LassoNet is:

[0027] ;

[0028] in, express Simulated brain magnetic signals at all times, Indicates The guiding field matrix for each age group, To be solved The source dipole vector at time instant, is the regularization parameter, is the constraint vector, are the model parameters that need to be updated.

[0029] Furthermore, step three is as follows:

[0030] Calculate the source intensity of each dipole based on the single-time source dipole vector to obtain a source intensity sequence;

[0031] The source intensity sequence is temporally encoded, and the temporal encoding results are sequentially input into the temporal self-attention module and the temporal convolutional network TCNs for encoding to obtain the temporal feature matrix;

[0032] Source space dimension encoding: Perform feature transformation on the single-time source dipole vector to obtain the transformed source dipole sequence, and convert the first to the second vector of each vector of the transformed source dipole sequence into Elements as The source signal component in the direction to Elements as The source signal component in the direction to Elements as Directional source signal component;

[0033] The transformed source dipole sequence is 3D position encoded in the neighborhood space dimension, and the temporal dimension encoding result and the 3D position encoding result are feature fused. The fused features are passed through the self-attention module and the feedforward neural network in turn, and the output of the encoder module is used as the key K and value V of the multi-head attention module in the decoder module.

[0034] Further, in step five, the expert system includes six network models and a gating unit. The six network models are set in parallel. The gating unit is used to select the network model. The selected network model is used for source signal correction. The pre-training of the expert system is performed in the Transformer architecture. The Transformer architecture only retains the multi-head attention module in the decoder module, the residual connection and the layer normalization module in the Transformer architecture during the expert system pre-training process. The pre-training process of the expert system is as follows:

[0035] Use step 1 to obtain the simulated source matrix and simulated brain magnetic signal used for expert system training, and use step 2 to obtain the single-time source dipole vector;

[0036] Different weight update strategies are used during pre-training and formal training: during pre-training, the weights of all network models and gating units are updated; during formal training, the weights of the selected network models and gating units are updated;

[0037] The output of the multi-head attention module is input into the expert system, the mean square error between the simulated source matrix and the temporary traceability results output by each network model in the expert system is calculated and arranged in ascending order, and the weights in the gate control unit and the network model are updated according to the mean square error result;

[0038] The gating unit output is used to select the probability of each network model. The smaller the mean square error, the higher the probability of the gating unit output, that is, the higher the probability of selecting the corresponding network model;

[0039] During formal training, the first three network models are selected based on the output probabilities of the gating units and arranged in descending order. The weighted sum of the temporary traceability results output by the selected network models and the output probabilities of the gating units is used as the reconstructed traceability results output by the expert system in the current round of training. The expert system is adjusted based on the reconstructed traceability results output each time.

[0040] Furthermore, the weighted sum of the temporary traceability result output by the selected network model and the output probability of the gate control unit is used as the temporary reconstructed traceability result. , temporarily reconstruct the traceability results The calculation process is as follows:

[0041] ;

[0042] in, is the output probability of the gate unit, which indicates the probability of the network model being selected and the weight of the output result of the selected network model. They represent the temporary tracing results output by the selected network model respectively. Respectively represent the index of the selected network model, For the analog source matrix, is the fusion feature coefficient, the initial state of pre-training , the network model is gradually updated during training, .

[0043] Furthermore, in step 4, the specific generation process of query Q is as follows:

[0044] The matrix after sampling the source Computing attention queries ,key ,value ;

[0045] Query ,key Perform discrete Fourier transform and inverse fast Fourier transform, calculate query ,key The frequency domain characteristics of the autocorrelation coefficient matrix are then calculated ;

[0046] Based on the autocorrelation coefficient matrix Sum Calculate the attention matrix, add the attention matrix to the mean of the simulated source matrix intensity to get the first attention feature, and use the first attention feature as the query Q.

[0047] Furthermore, the autocorrelation coefficient matrix The calculation is as follows:

[0048] ;

[0049] in, represents the inverse fast Fourier transform, represents the discrete Fourier transform, represents the matrix transpose, Represents matrix product.

[0050] Furthermore, in the step 1, the improved neural group model is based on the single cortical column and double cortical column models proposed by Jansen-Rit, and adjusts the excitatory feedback transmission process between different sources based on the occurrence time of brain stimulating electrophysiological activities;

[0051] Each partition is composed of multiple cortical columns. First, the central cortical column is selected, and the central cortical column and the adjacent cortical columns are divided into the same area; the same cortical column contains multiple pyramidal cells and other cells, and other cells will transmit excitatory or inhibitory feedback to the pyramidal cells receiving feedback; other cortical columns in the same area will transmit excitatory feedback to the central cortical column;

[0052] In single-source simulation, only the feedback information of the central cortical column and the adjacent cortical columns needs to be considered. Based on prior physiological knowledge, the excitatory feedback of the cortical columns in other areas is adjusted to a fixed value.

[0053] When simulating multiple source points: the central cortical column of each source point calculates the electrical pulse feedback of other source points; for each selected source point, if other source points have stimulating electrophysiological activities synchronously with the central source point, the selected source point will receive excitatory feedback with a higher weight than other source points; if the stimulating electrophysiological activities of other source points and the central source point occur asynchronously, the weight will be attenuated according to the time interval between the stimulating electrophysiological activities of the selected source point and other source points. The longer the interval, the greater the attenuation coefficient.

[0054] Furthermore, for the pyramidal cells receiving feedback in the same cortical column, specifically: the electrical signal fed back by the cerebral cortex is converted into the average pulse density of the emitted action potential by a conversion function, and the conversion function is set as follows according to the period of brain electrophysiological activity:

[0055] Before stimulating neural activity, the hyperbolic tangent function is used to simulate the slow rise and fall of potential changes, and the conversion function It is expressed as:

[0056] ;

[0057] At the beginning of the stimulating neural activity, the potential mutation is approximated by the inverse tangent function, and the conversion function It is expressed as:

[0058] ;

[0059] The square wave burst activity when the sudden change potential continues uses the activation function to ensure that the potential continues with fluctuations at the peak value, and the conversion function It is expressed as:

[0060] ;

[0061] The fluctuation process after the potential mutation It is expressed as:

[0062] ;

[0063] in, Determine the maximum firing rate of a nerve cell, represents the postsynaptic potential at 50% of the maximum discharge rate, It plays a role in regulating the gradient of the model. represents feedback from cells other than pyramidal cells;

[0064] Excitatory feedback It is expressed as:

[0065] ;

[0066] Inhibitory feedback It is expressed as:

[0067] ;

[0068] in, and Determines the maximum amplitude of excitatory and inhibitory postsynaptic potentials; and represents the sum of the inverse of the passive membrane time constant and other spatially distributed delays in the dendritic network, Indicates time, Represents a natural constant.

[0069] The advantages of the magnetoencephalogram (MEG) tracing method based on a deep neural network provided by the present invention are: the MEG tracing method based on a deep neural network provided in the structure of the present invention, the embodiment improves the Jansen-Rit algorithm when simulating training data, so that the simulation results are more in line with neurophysiological characteristics, based on the division of age groups, a more suitable brain model is used for different tracing objects, the single-time tracing uses LassoNet based on the multiplier alternating direction method to select features, reduces unnecessary features, and accelerates model convergence, the source reconstruction module is based on the Transformer architecture, and a cross-multidimensional attention mechanism based on the temporal dimension and the source space dimension is used in the Transformer encoder module, and a temporal convolutional network is used in the temporal dimension to process the temporal attention, and the spatial dimension attention integrates the temporal dimension attention features, and the inverse fast Fourier transform (IFFT) is used in the Transformer decoder module to calculate the autocorrelation coefficient, which reduces the amount of calculation, and the expert system is used to correct the source intensity to obtain a more accurate reconstruction and tracing result. BRIEF DESCRIPTION OF THE DRAWINGS

[0070] Figure 1 It is a structural schematic diagram of the present invention;

[0071] Figure 2 It is a structural diagram of the traceability model;

[0072] Figure 3 It is the Transformer structure diagram;

[0073] Figure 4 Training flowchart for the traceability model;

[0074] Figure 5 This is the inference flowchart of the trained provenance tracing model. DETAILED DESCRIPTION

[0075] Below, the technical solution of the present invention is described in detail through specific embodiments. Many specific details are set forth in the following description to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific implementation disclosed below.

[0076] Figures 1 to 5 As shown in the figure, the invention proposes a method for tracing the origin of magnetoencephalogram based on deep neural network. The reasoning and execution process after the tracing model is trained are as follows: Figure 5 As shown, specifically: the magnetoencephalogram data is acquired through a preset head model, the magnetoencephalogram data is traced at a single moment through the neural network framework LassoNet to obtain the source dipole vector corresponding to each moment, multiple source dipole vectors are combined in chronological order to obtain a source signal sequence, the source signal sequence is input into the trained source reconstruction model, and the reconstruction tracing result is output. The reconstruction tracing result may contain one or more source signals of different intensities. Based on the reconstruction tracing result and the head model, the source dipole matrix corresponding to the magnetoencephalogram signal is obtained, and the source position and source intensity are calculated according to the source dipole matrix.

[0077] Figures 1 to 4 As shown in the figure, the training process of the traceability model is as follows:

[0078] Step 1: Based on the dynamic changes of cortical gray matter, cortical white matter and overall brain volume of people of different age groups during human growth, corresponding head models are established for people of different age groups. The head models are used to establish the guidance field matrix required for tracing the source. The simulation source matrix is ​​calculated based on the improved neural group model, and the guidance field matrix is ​​multiplied by the simulation source matrix to obtain the simulated brain magnetic signal;

[0079] Based on the research results of the paper "Brain charts for the human lifespan", a relatively accurate age group division method can be obtained. The study summarized the head structure magnetic resonance imaging of 101,457 volunteers from around the world, whose ages ranged from 115 days after conception to 100 years old. The study performed GAMLSS model fitting on MRI measurement samples of total cortical gray matter volume (Grey Matter Volume, GMV), total white matter volume (White Matter Volume, WMV), total subcortical gray matter volume (subcortical Grey Matter Volume, sGMV) and total cerebrospinal fluid volume (Cerebrospinal FluidVolume, CSF). Based on the conclusions in the paper, this embodiment divides the age distribution of head modeling into: 0 to 3 years old, 3 to 14.4 years old, 14.4 to 28.7 years old, 28.7 to 40 years old, 40 to 55 years old, and 55 to 80 years old. This embodiment combines the calculation of relevant parameters in the three-shell Boundary Element Method model to adjust the modeling process of the above-mentioned people of different age groups.

[0080] The modeling of the traceability problem usually considers GMV, WMV and overall brain volume indicators. In this embodiment, the weights of sGMV and CSF are adjusted to 0.1, and the overall brain volume is added as a new parameter. After model calculation, the GMV, WMV and overall brain volume of children aged 0 to 3 years old will have a large increase. In this embodiment, children aged 0 to 3 years old are further divided, with 6 months as the time separation point. In each time interval, head models are established for two children respectively, and then the average value of 12 different head model data is calculated, and the new head model is used as the head modeling of the age group. The 3 to 14.4-year-old population has similar dynamic changes as the 0 to 3-year-old children, but the changes are slowed down, the time interval is expanded to 3 years, and the average value is finally calculated. There is no obvious change in the overall brain volume of the 14.4 to 55-year-old population, and the GMV and WMV volumes have slow dynamic changes with age. In this embodiment, 6 samples are taken in each age group, and the average value is calculated by modeling. The GMV, WMV and overall brain volume of the 55 to 80-year-old population have a certain degree of dynamic reduction. In this embodiment, the time interval is adjusted to 5 years, and the average value is calculated by modeling respectively.

[0081] This embodiment is based on the three-layer BEM model and the above process to establish different source space models, set the cerebral cortex dipoles to 11948, and thus set the brain partition. The setting of brain partition needs to take into account the robustness and versatility of the model. This embodiment adopts the method of using the same number of partitions for 6 age groups, which can reduce the number of parameters used in the model. Since the head volume, shape and other characteristics of people of different age groups are different, the area of ​​the head model division area of ​​different people is different. The number of dipoles in each area is set to about 11, and the source space is divided into 1024 areas. The number of partitions used in this embodiment can be adjusted according to actual conditions. Different partition strategies can be adopted in different tracing tasks. For example, using relevant prior knowledge, the source area can be limited to a smaller space. If the model needs to have different partition numbers for people of different age groups, different model training labels can be set for each age group. This setting will increase certain time and space costs. In this embodiment, the magnetoencephalogram uses a 306-channel measuring instrument. Based on the above settings, a guidance field matrix for people of different age groups can be established ,in, Indicates the number of channels, =1024 represents the number of source space partitions, and 3 represents the three directions of the source dipole.

[0082] Compared with the traditional method, the modeling method according to different age groups in this embodiment can adapt to the changes in the traceability object, ensure the robustness of the model, and further improve the accuracy of model traceability through targeted modeling. At the same time, the method of dividing by age groups reduces the amount of calculation compared to modeling by individual.

[0083] The improved neural group model of this embodiment, namely the improved Jansen-Rit algorithm, is based on the single cortical column and double cortical column models proposed by Jansen-Rit, and adjusts the excitatory feedback transmission process between different sources based on the occurrence time of brain stimulating electrophysiological activities;

[0084] Each partition is composed of multiple cortical columns. First, the central cortical column is selected, and the central cortical column and its adjacent cortical columns are divided into the same area; the same cortical column contains multiple pyramidal cells and other cells (other pyramidal cells, stellate cells or basket cells, etc.), which will transmit excitatory or inhibitory feedback to the pyramidal cells receiving feedback; other cortical columns in the same area will transmit excitatory feedback to the central cortical column;

[0085] When simulating a single source point, only the feedback information within the central cortical column and the adjacent cortical columns needs to be considered. Based on prior physiological knowledge, the excitatory feedback of the cortical columns in other areas is adjusted to a fixed value.

[0086] When simulating multiple source points, the central cortical column of each source point needs to calculate the electrical pulse feedback of other source points; for each selected source point, if other source points and the central source point have stimulating electrophysiological activities synchronously, the selected source point will receive excitatory feedback with a higher weight than other source points; if the stimulating electrophysiological activities of other source points and the central source point occur asynchronously, the weight will be attenuated according to the time interval between the stimulating electrophysiological activities of the selected source point and other source points. The longer the interval, the greater the attenuation coefficient.

[0087] This embodiment combines the characteristics of brain electrophysiological activities, especially the characteristics of brain neurological diseases such as epilepsy, and integrates the potential fluctuation characteristics before, during, during and after neural activity in the cerebral cortex to modify the Jansen-Rit algorithm.

[0088] For the pyramidal cells receiving feedback in the same cortical column: the electrical signal fed back from the cerebral cortex is converted into the average pulse density of the emitted action potential by a conversion function, and the conversion function is set as follows according to the period of brain electrophysiological activity:

[0089] Before the stimulating neural activity, the hyperbolic tangent function was used to simulate the slow rise and slow fall of the potential change, and the conversion function was expressed as:

[0090] ;

[0091] At the beginning of the stimulating neural activity, the potential mutation is approximated by the inverse tangent function, and the conversion function is expressed as:

[0092] ;

[0093] The square wave burst activity when the mutation potential continues uses the SoftSign function (an activation function) to ensure that the potential continues with fluctuations at the peak. The conversion function is expressed as:

[0094] ;

[0095] The fluctuation process after the potential mutation is expressed as:

[0096] ;

[0097] in, Determine the maximum firing rate of a nerve cell, represents the postsynaptic potential at 50% of the maximum discharge rate, It plays a role in regulating the gradient of the model. represents feedback from cells other than pyramidal cells;

[0098] Excitatory feedback It is expressed as:

[0099] ;

[0100] Inhibitory feedback It is expressed as:

[0101] ;

[0102] in, and Determines the maximum amplitude of excitatory and inhibitory postsynaptic potentials; and represents the sum of the inverse of the passive membrane time constant and other spatially distributed delays in the dendritic network, Indicates time, Represents a natural constant.

[0103] Different source activity potentials can be obtained by characterizing the connection constants of the interactions between pyramidal cells and excitatory and inhibitory interneurons (other cells). By changing the value of the constant, brain activity simulation signals with different characteristics can be obtained.

[0104] Based on the above process, the analog source signal of each partition can be obtained. The analog source signal at a certain moment is Indicates that the analog source matrix is ​​spliced ​​in time sequence ,in represents the real number space. In this embodiment, Indicates the number of source space partitions, Indicates the time sampling rate, that is, the duration of the simulated EEG signal, which can be adjusted according to the specific scenario. In this embodiment, =1024, =600.

[0105] Based on the above steps, we have obtained: head models of people of different age groups, including the location of cortical partitions, the adjacent relationships between different partitions and the dipole information contained in each partition; a guidance field matrix established based on the magnetoencephalography instrument and the head model, and a simulation source matrix based on the head model.

[0106] Based on the guidance field matrix and the simulated source matrix, the simulated EEG signal is calculated for single-time source tracing. The construction of the simulated EEG signal can be expressed as:

[0107] ;

[0108] in, It is the guiding field matrix for people of different age groups. Indicates the number of channels of the EEG or MEG measuring instrument. This embodiment uses a MEG acquisition instrument with 306 channels. Indicates the number of source space partitions. This embodiment uses 1024 partitions. Indicates the time sampling rate, which can be adjusted according to the specific scenario; To simulate the magnetic brain signal for the input of the tracing model.

[0109] Step 2: Use the neural network framework LassoNet to trace the source of each moment of the simulated EEG signal and output the single-moment source dipole vector;

[0110] Based on "LassoNet: A Neural Network with Feature Sparsity", this embodiment uses LassoNet to solve the single-time source signal, reduces the number of parameters in the calculation process and obtains a sparse solution.

[0111] LassoNet is a neural network framework with global feature selection, which achieves feature sparsity by adding a residual connection and allows features to participate in hidden layer feature value calculation only when the input feature represented by its residual connection is active.

[0112] LassoNet is often used to solve and optimize ill-posed problems. Based on the least absolute shrinkage and selection operator method, LassoNet is mainly used for variable selection and regularization of regression analysis and has strong generalization ability. That is:

[0113] ;

[0114] in, Represents the variable vector to be optimized, corresponding to the source vector in the tracing problem , represents a real number, represents the response variable vector, corresponding to the observation vector in the tracing problem, i.e., the brain magnetic vector, Represents the characteristic matrix, corresponding to the guidance field matrix in the tracing problem , represents the regularization parameter, represents the two-norm, Represents a norm.

[0115] LassoNet has a good fitting effect on single-time source dipoles and can screen non-source weak features during the solution process.

[0116] This embodiment modifies the optimization objective of the original LassoNet, and uses the Alternating Direction Method of Multipliers (ADMM) to divide the optimization term and the regularization term into two parts, constructs an augmented Lagrangian function, introduces dual residuals, increases the model's learnable parameters, reduces artificially determined constraints, and the LassoNet optimized based on the ADMM method converges faster.

[0117] In this embodiment, the original LassoNet optimization process is modified as follows:

[0118] ;

[0119] in, express The input with age label at all times is the vector corresponding to the simulated EEG signal. To be solved The source dipole vector at time instant, is a constraint vector satisfying , are the model parameters that need to be updated, is a regularization parameter with an initial value of 0.1, which decreases as the number of training iterations increases.

[0120] For the above problem, LassoNet is used to solve the source dipole vector at each moment The solution result is 600 source dipole vectors, and the single-time source dipole vector is obtained based on the neural network framework LassoNet.

[0121] Considering only the results at a single moment, LassoNet is used to output the source dipoles of each partition at each moment. Compared with the method of tracing the source using convolutional neural networks or deep neural networks, LassoNet can screen useful features, reduce the parameters of the calculation process, and obtain sparse solutions;

[0122] Step 3: Input the single-time source dipole vector and the simulated source matrix into the Transformer architecture. The encoder module in the Transformer architecture encodes the single-time source dipole vector in the temporal dimension and the source space dimension in stages. The source dipole strength is calculated based on the single-time source dipole vector and temporal encoding is performed. The temporal encoding result is sequentially passed through the temporal self-attention module and the temporal convolutional network (TCNs) to obtain the temporal feature matrix. The temporal feature matrix and the source space dimension encoding result are feature fused. The fused feature matrix is ​​sequentially passed through the self-attention module and the feedforward neural network. The output of the encoder module is used as the key K and value V of the multi-head attention module in the decoder module.

[0123] Perform temporal encoding and 3D position encoding on the single-time source dipole vector (3D position encoding is the source space dimension encoding), use cross attention on the single-time source dipole vector in the temporal dimension and the spatial dimension in the encoder module, and replace the traditional feedforward neural network with temporal convolutional networks (TCNs) after the temporal attention module;

[0124] In the encoder module:

[0125] (a1) Calculate the source strength of each dipole based on the single-time source dipole vector, encode the source strength of each dipole in the time dimension, and obtain a source strength sequence encoding vector;

[0126] When encoding in the time series dimension, only the intensity change of the source is considered, and the single-time source dipole vector is obtained by step 2 , the strength of each source dipole can be found:

[0127] ;

[0128] in is the source strength for each dipole, Represent the three components of each source dipole in three-dimensional space.

[0129] Based on the above steps, for each single moment source dipole vector Calculate the source intensity vector .

[0130] (a2) The sequence of source intensity vectors is input into the temporal attention mechanism module and the temporal convolutional network TCNs in sequence to obtain the temporal dimension encoding result;

[0131] The Transformer encoder module inputs the source intensity sequence into the temporal attention mechanism module in sequence based on the cross-attention mechanism:

[0132] ;

[0133] ;

[0134] in, represents the timing code, A temporal encoding vector representing the source strength, Represents the attention mechanism, through the attention mechanism Processing is performed to obtain the source intensity sequence attention vector ; , is the time series embedding dimension. In this embodiment, =1024, , =128 is the attention head dimension.

[0135] Use temporal convolutional networks (TCNs) to replace the feedforward neural network in Transfomer, and use the temporal convolutional network to process the source intensity sequence attention vector to obtain the temporal dimension encoding result:

[0136] ;

[0137] in Represents the encoding result of the time series dimension, represents a temporal convolutional network, = 3072 is the dimension of the source dipole vector space.

[0138] (a3) Perform characteristic transformation on the single-time source dipole vector to obtain the transformed source dipole sequence, and replace the first to the last vector in each vector of the transformed source dipole sequence with the Elements as The source signal component in the direction to Elements as The source signal component in the direction to Elements as Directional source signal component;

[0139] (a4) performing 3D position encoding on the transformed single-moment source dipole sequence in the neighborhood space dimension to obtain a spatial sequence encoding vector, and using the spatial sequence encoding vector as the 3D position encoding result;

[0140] In the neighborhood space dimension, the transformed single-time source location sequence is Perform 3D position encoding:

[0141] ;

[0142] ;

[0143] ;

[0144] ;

[0145] ;

[0146] ;

[0147] in represents the three components of the source dipole in the source space, Indicates the number information of each component, and their value range is ,in is the dimension of the source dipole vector space.

[0148] (a4) The temporal dimension encoding result and the 3D position encoding result are feature fused. The fused features are passed through the self-attention module and the feedforward neural network in turn, and the output of the encoder is used as the key K and value V of the multi-head attention module in the decoder module:

[0149] ;

[0150] ;

[0151] in, represents the 3D position encoding, It represents the feature vector obtained by fusion of the temporal dimension encoding result and the 3D position encoding result. It represents the feature vector obtained by fusing the temporal dimension encoding result and the 3D position encoding result through the self-attention module in sequence. After passing through the fully connected layer and layer normalization module, it serves as the output of the encoder module.

[0152] Step 4: The decoder module in the Transformer architecture performs time-series coding on the analog source matrix, and performs source sampling on the analog source matrix after time-series coding. The sampled analog source matrix is ​​respectively calculated and queried. ,key ,value , and through the inverse transformation module of the fast Fourier transform, the above results are transformed to obtain the query Q;

[0153] In the encoder module: The decoder module uses discrete Fourier transform (DFT) and inverse fast Fourier transform (IFFT) to replace the masked multi-head self-attention module in the traditional decoder module, and replaces the feedforward neural network of the decoder module with an expert system with a gating unit.

[0154] During training, the decoder module inputs the analog source matrix , and use the source signal to calculate the source intensity matrix , encode it and get the encoded matrix Based on the sparsity of the source signal, the inverse transform of the fast Fourier transform is used to calculate the autocorrelation coefficient. Since the source signal matrix is ​​not periodic and the strong source signal is stored in the encoded matrix in a discrete manner First, use the length ( )'s window pair encoding matrix Perform source sampling, extract the source of each encoded matrix and the area affected by the main source, and use discrete Fourier transform on the sampled matrix.

[0155] The matrix after sampling of each source is obtained by sampling ,in is the number of samples, Less than , in this embodiment Usually around 10. Calculate the sampled matrix Value, calculation query ,key The autocorrelation coefficient matrix :

[0156] ;

[0157] in, , represents matrix product, represents the matrix transpose, represents the inverse fast Fourier transform, represents discrete Fourier transform.

[0158] Calculate the attention matrix:

[0159] ;

[0160] Attention Matrix When calculating, change the dimension from Reduce to .

[0161] Then add the original intensity mean to the attention matrix output, and restore the time dimension to , transform the above results to obtain query Q.

[0162] Step 5. Based on the key K, value V output from step 3 and query Q output from step 4, calculate the multi-head attention feature, input the feature map output by the multi-head attention module into the pre-trained expert system to reconstruct the source signal, output the decoding attention feature, use the residual connection to fuse the multi-head attention feature and the decoding attention feature, and then perform layer normalization. The normalized result continues to be input into the next decoder module. If the last decoder module is reached, the normalized result of the last layer is linearly transformed to obtain the reconstructed tracing result.

[0163] Pre-train an expert system to replace the feedforward neural network module in the Transformer decoder module. The expert system includes six network models and a gating unit. The six network models are set in parallel. The network models respectively use a time series convolutional network, a neural network framework LassoNet, and a feedforward neural network. Each network model has two different depths. In this embodiment, the time series convolutional network uses 3 hidden layers and 6 hidden layers, respectively, the neural network framework LassoNet uses 5 hidden layers and 10 hidden layers, respectively, and the feedforward neural network uses 4 hidden layers and 8 hidden layers. The gating unit uses a simple multi-layer perceptron. This embodiment uses a multi-layer perceptron with 4 hidden layers. In this embodiment, the specific settings of the six network models can be two time series convolutional networks, two neural network frameworks LassoNet, and two feedforward neural networks. In the actual training process, the six network models are selected by the gating unit, and the output results of the selected network models are used as the traceability results of the current round to adjust the expert system. The adjustment of the expert system includes the adjustment of the model parameters of the network model and the weight update of the gating unit.

[0164] The expert system is pre-trained in a supervised manner. The pre-training of the expert system is performed in the Transformer architecture. During the pre-training of the expert system, the Transformer architecture only retains the multi-head attention module in the decoder module, the residual connection and layer normalization modules in the Transformer architecture. That is, the modules closed during the pre-training of the expert system are: the temporal self-attention module of the temporal encoding module, the temporal convolutional network, the self-attention module of the 3D position encoding module, the feedforward neural network, and the IFFT module. The retained modules are: the residual connection and layer normalization modules in the overall Transformer architecture, and the multi-head self-attention module in the decoder module. These retained residual connections and layer normalization modules are used to ensure that the feature map dimensions in the pre-training process are consistent with the feature map dimensions in the formal training. The pre-training process of the expert system is as follows:

[0165] (b1) Use step 1 to obtain the simulated source matrix and simulated EEG signal used for expert system pre-training, and use step 2 to obtain the single-time source dipole vector;

[0166] (b2) inputting the single-time source dipole vector and the analog source matrix into the Transformer architecture, wherein the encoder module in the Transformer architecture performs temporal dimension encoding and 3D position encoding on the single-time source dipole to obtain a source intensity sequence encoding vector and a 3D position encoding result;

[0167] (b3) The source intensity sequence encoding vector and the 3D position encoding result are input into the module retained by the encoder module, and the output of the encoder module is used as the pre-trained key ,value ; Encode the source intensity sequence into the module retained by the decoder module and output the pre-trained query ; Query ,key ,value After passing through the multi-head attention module, the output of the multi-head attention module is used as the fusion feature;

[0168] (b4) Input the fusion features into the expert system, calculate the mean square error between the simulated source matrix and the temporary traceability results output by each network model in the expert system, and arrange them in ascending order (the smaller the mean square error, the better the result), and update the weights in the gating unit and the network model according to the mean square error result. The probability of selecting each network model output by the gating unit is as small as the mean square error, and the higher the probability of the gating unit output, that is, the higher the probability of selecting the corresponding network model. During formal training, the first three network models are selected according to the output probability of the gating unit and arranged in descending order. The weighted sum of the temporary traceability results output by the selected model and the probability of the gating unit output is used as the temporary reconstructed traceability result output during expert pre-training. ;

[0169] Temporary reconstruction of traceability results The calculation process is as follows:

[0170] ;

[0171] in, is the output probability of the gate unit, which indicates the probability of the network model being selected and the weight of the output result of the selected network model. They represent the temporary tracing results output by the selected network model respectively. Respectively represent the index of the selected network model, is the analog source signal matrix, is the fusion feature coefficient, the initial state of pre-training , .

[0172] The pre-trained expert system replaces the feedforward neural network module in the Transformer decoder module and is deployed in the traceability model to assist in the effective output of the traceability results.

[0173] Through steps one to five, this embodiment uses a staged attention mechanism based on the temporal dimension and the source space dimension in the Transformer encoder module, uses a temporal convolutional network to process the temporal attention in the temporal dimension, and the spatial dimension attention integrates the temporal dimension attention features. The inverse transform of the fast Fourier transform is used in the Transformer decoder module to calculate the autocorrelation coefficient, which reduces the amount of calculation, uses an expert system to correct the source, and obtains more accurate results.

[0174] When establishing the forward process, this embodiment fully balances the performance and computational complexity of the tracing model, and models the brains of people of different age groups separately; when generating training data, the Jansen-Rit model is improved, and different conversion functions are used for different stages of cortical source activity to simulate source signals that are more in line with neurophysiology; the reverse process fully considers the temporal information of the source and the neighborhood relationship of different sources, first using LassoNet to perform preliminary tracing of the dipoles, and then using Transformer to build dependencies in the temporal dimension and the adjacent source space dimension, and correcting the dipoles, thereby improving the output quality of the reconstructed tracing results.

[0175] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can make equivalent replacements or changes according to the technical scheme and inventive concept of the present invention within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims

1. A magnetoencephalogram tracing method based on deep neural network, It is characterized in that The magnetic brain data is obtained through a preset head model, and the magnetic brain data is traced at a single moment through the neural network framework LassoNet to obtain the source dipole vector corresponding to each moment. Multiple source dipole vectors are combined in time order to obtain a source signal sequence, which is input into the trained tracing model and the reconstructed tracing result is output. Based on the reconstructed tracing result and the head model, the source dipole matrix corresponding to the magnetic brain signal is obtained, and the source position and source intensity are calculated according to the source dipole matrix. The training process of the traceability model is as follows: Step 1: Establish different head models according to the head structure magnetic resonance imaging of people of different age groups, use the head models to establish the guidance field matrix required for tracing, calculate the simulation source matrix based on the improved neural group model, and multiply the guidance field matrix with the simulation source matrix to obtain the simulated brain magnetic signal; Step 2: Use the neural network framework LassoNet to trace the source of each moment of the simulated EEG signal and output the single-moment source dipole vector; Step 3: Input the single-time source dipole vector and the simulated source matrix into the Transformer architecture. The encoder module in the Transformer architecture encodes the single-time source dipole vector in the time dimension and the source space dimension in stages. The source dipole strength is calculated based on the single-time source dipole vector and time encoding is performed. The time encoding result is sequentially passed through the time self-attention module and the time convolutional network to obtain the time feature matrix. The time feature matrix and the source space dimension encoding result are feature fused. The fused feature matrix is ​​sequentially passed through the self-attention module and the feedforward neural network. Then the output of the encoder module is used as the key K and value V of the multi-head attention module in the decoder module. Step 4: The decoder module in the Transformer architecture performs time-series coding on the analog source matrix, and performs source sampling on the analog source matrix after time-series coding. The sampled analog source matrix is ​​respectively calculated and queried. ,key ,value , and after the inverse fast Fourier transform module, the query ,key ,value Transform to obtain query Q; Step 5: Based on the key K and value V output from step 3 and the query Q output from step 4, the multi-head attention feature is calculated, the feature map output by the multi-head attention module is input into the pre-trained expert system for source signal reconstruction, the decoding attention feature is output, the multi-head attention feature and the decoding attention feature are fused using residual connection, and then layer normalization is performed. The normalized result continues to be input into the next decoder module. If the last decoder module is reached, the normalized result of the last layer is linearly transformed to obtain the reconstructed tracing result. In step 1, different head models are established for the head structure magnetic resonance imaging of people of different age groups. The process is as follows: According to the age distribution, the brain modeling age groups were divided into: 0 to 3 years old, 3 to 14.4 years old, 14.4 to 28.7 years old, 28.7 to 40 years old, 40 to 55 years old, and 55 to 80 years old; According to age groups, different cerebral cortical source models are established, and the number of regional divisions and the location of the regional center of each cerebral cortical source model are determined; The cerebral cortex is divided into different areas according to the head model, and the electrical signals fed back from the cerebral cortex are converted into the average pulse density of the emitted action potential to obtain the simulated source signals of each partition, and the simulated source signals of each partition are spliced ​​according to the time sequence of the source activity to obtain the simulated source matrix; In addition, in the step 1, the improved neural group model is based on the single cortical column and double cortical column models proposed by Jansen-Rit, and adjusts the excitatory feedback transmission process between different sources based on the occurrence time of brain stimulating electrophysiological activities; Each partition is composed of multiple cortical columns. First, the central cortical column is selected, and the central cortical column and the adjacent cortical columns are divided into the same area; the same cortical column contains multiple pyramidal cells and other cells, and the other cells transmit excitatory or inhibitory feedback to the pyramidal cells receiving feedback; other cortical columns in the same area transmit excitatory feedback to the central cortical column; In single-source simulation, only the feedback information of the central cortical column and the adjacent cortical columns needs to be considered. Based on prior physiological knowledge, the excitatory feedback of the cortical columns in other areas is adjusted to a fixed value. When simulating multiple source points: the central cortical column of each source point calculates the electrical pulse feedback of other source points; for each selected source point, if other source points have stimulating electrophysiological activities synchronously with the central source point, the selected source point will receive excitatory feedback with a higher weight than other source points; if the stimulating electrophysiological activities of other source points and the central source point occur asynchronously, the weight will be attenuated according to the time interval between the stimulating electrophysiological activities of the selected source point and other source points. The longer the interval, the greater the attenuation coefficient.

2. The method for tracing the origin of magnetoencephalogram based on deep neural network according to claim 1, It is characterized in that In step 2, the optimization goal of the neural network framework LassoNet is: in, express Simulated brain magnetic signals at all times, Indicates The guiding field matrix for each age group, To be solved The source dipole vector at time instant, is the regularization parameter, is the constraint vector, are the model parameters that need to be updated, Represents matrix transpose.

3. The method for tracing the origin of magnetoencephalogram based on deep neural network according to claim 1, It is characterized in that Step 3 is as follows: Calculate the source intensity of each dipole based on the single-time source dipole vector to obtain a source intensity sequence; The source intensity sequence is temporally encoded, and the temporal encoding results are sequentially input into the temporal self-attention module and the temporal convolutional network for encoding to obtain the temporal feature matrix; Source space dimension encoding: Perform feature transformation on the single-time source dipole vector to obtain the transformed source dipole sequence, and convert the first to the second vector of each vector of the transformed source dipole sequence into Elements as The source signal component in the direction to Elements as The source signal component in the direction to Elements as Directional source signal component; The transformed source dipole sequence is 3D position encoded in the neighborhood space dimension, and the temporal dimension encoding result and the 3D position encoding result are feature fused. The fused features are passed through the self-attention module and the feedforward neural network in turn, and the output of the encoder module is used as the key K and value V of the multi-head attention module in the decoder module.

4. The method for tracing the origin of magnetoencephalogram based on deep neural network according to claim 1, It is characterized in that In step 5, the expert system includes six network models and a gating unit. The six network models are set in parallel. The gating unit is used to select the network model. The selected network model is used for source signal correction. The pre-training of the expert system is performed in the Transformer architecture. During the expert system pre-training process, the Transformer architecture only retains the multi-head attention module in the decoder module, the residual connection and the layer normalization module in the Transformer architecture. The pre-training process of the expert system is as follows: Use step 1 to obtain the simulated source matrix and simulated brain magnetic signal used for expert system training, and use step 2 to obtain the single-time source dipole vector; Different weight update strategies are used during pre-training and formal training: during pre-training, the weights of all network models and gating units are updated; During formal training, the weights of the selected network model and gating units are updated; The output of the multi-head attention module is input into the expert system, the mean square error between the simulated source matrix and the temporary traceability results output by each network model in the expert system is calculated and arranged in ascending order, and the weights in the gate control unit and the network model are updated according to the mean square error result; The gating unit output is used to select the probability of each network model. The smaller the mean square error, the higher the probability of the gating unit output, that is, the higher the probability of selecting the corresponding network model; During formal training, the first three network models are selected based on the output probabilities of the gating units and arranged in descending order. The weighted sum of the temporary traceability results output by the selected network models and the output probabilities of the gating units is used as the reconstructed traceability results output by the expert system in the current round of training. The expert system is adjusted based on the reconstructed traceability results output each time.

5. The method for tracing the origin of magnetoencephalogram based on deep neural network according to claim 4, It is characterized in that The weighted sum of the temporary traceability result output by the selected network model and the output probability of the gate control unit is used as the temporary reconstructed traceability result. , temporarily reconstruct the traceability results The calculation process is as follows: in, is the output probability of the gate unit, which indicates the probability of the network model being selected and the weight of the output result of the selected network model. They represent the temporary tracing results output by the selected network model respectively. Respectively represent the index of the selected network model, For the analog source matrix, is the fusion feature coefficient, the initial state of pre-training , the network model is gradually updated during training, .

6. The method for tracing the origin of magnetoencephalogram based on deep neural network according to claim 1, It is characterized in that In step 4, the specific generation process of query Q is as follows: The matrix after sampling the source Computing attention queries ,key ,value ; Query ,key Perform discrete Fourier transform and inverse fast Fourier transform, calculate query ,key The frequency domain characteristics of the autocorrelation coefficient matrix are then calculated ; Based on the autocorrelation coefficient matrix Sum Calculate the attention matrix, add the attention matrix to the mean of the simulated source matrix intensity to get the first attention feature, and use the first attention feature as the query Q.

7. The method for tracing the origin of magnetoencephalogram based on deep neural network according to claim 6, It is characterized in that in, Autocorrelation coefficient matrix The calculation is as follows: in, represents the inverse fast Fourier transform, represents the discrete Fourier transform, represents the matrix transpose, Represents matrix product.

8. The method for tracing the origin of magnetoencephalogram based on deep neural network according to claim 1, It is characterized in that For the pyramidal cells receiving feedback in the same cortical column, specifically: the electrical signal fed back by the cerebral cortex is converted into the average pulse density of the emitted action potential by a conversion function, and the conversion function is set as follows according to the period of brain electrophysiological activity: Before stimulating neural activity, the hyperbolic tangent function is used to simulate the slow rise and fall of potential changes, and the conversion function It is expressed as: At the beginning of the stimulating neural activity, the potential mutation is approximated by the inverse tangent function, and the conversion function It is expressed as: The square wave burst activity when the sudden change potential continues uses the activation function to ensure that the potential continues with fluctuations at the peak value, and the conversion function It is expressed as: The fluctuation process after the potential mutation It is expressed as: in, Determine the maximum firing rate of a nerve cell, represents the postsynaptic potential at 50% of the maximum discharge rate, It plays a role in regulating the gradient of the model. represents feedback from cells other than pyramidal cells; Excitatory feedback It is expressed as: Inhibitory feedback It is expressed as: in, and Determines the maximum amplitude of excitatory and inhibitory postsynaptic potentials; and represents the sum of the inverse of the passive membrane time constant and other spatially distributed delays in the dendritic network, Indicates time, Represents a natural constant.

Images