A method for magnetoencephalic source localization based on nearest neighbor multivariate linear regression model

Abstract

The application provides a brain magnetic source positioning method based on a nearest neighbor multivariate linear regression model. The method comprises the following steps: constructing a brain model based on the nearest neighbor multivariate linear regression model and representing the space-time correlation between brain source activities; constructing a source covariance matrix in the model based on a graph Laplacian smoother, and representing the space correlation of the covariance hierarchy between the source activities; and making reasonable assumptions for the model based on Bayesian inference and convex optimization theory, and deriving a formula for iterative solution of the model parameters. The source positioning algorithm realized by the application can obtain more accurate source estimation, including accurate positioning of the source position, the range size and the source time waveform estimation, and provides a solid technical foundation for medical applications of the brain magnet.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of biomedical engineering, and particularly relates to a magnetoencephalography source localization method based on nearest neighbor multivariate linear regression model. BACKGROUND

[0002] Magnetoencephalography (MEG) is a non-invasive functional imaging technique for exploring human brain activity, which has attracted much attention in neuroscience and clinical research in recent decades. Compared with some other functional imaging techniques (such as electroencephalogram and near-infrared imaging), MEG can more accurately locate brain activity due to its high temporal and spatial resolution. Magnetoencephalography source localization refers to reconstructing potential brain activity from the magnetoencephalography signals obtained at the sensor. There are two main methods for MEG source localization: dipole fitting method and distributed source imaging method. The dipole fitting method assumes that there are a limited number of dipoles in the brain, and then estimates the parameters of these dipoles by least squares fitting. In the distributed source imaging method, the cortical surface is divided into a certain number of triangular meshes, and each mesh represents a dipole. Compared with the dipole fitting method, the imaging method determines the number of sources in advance, and only needs to solve a linear problem. Since the number of unknown parameters far exceeds the equation, there are still many solutions to this linear problem. Proper prior constraints should be added to obtain a unique solution.

[0003] Traditional distributed source imaging methods usually use L2 norm constraints, such as minimum norm estimation (MNE). MNE constrains the L2 norm of the current source to obtain the solution with the minimum energy. In the level of sensor signals, the magnetic field strength generated by shallow sources is larger than that generated by deep sources, so MNE is biased towards shallow sources. Variants of MNE have appeared in some existing literature, including weighted MNE (wMNE) and low resolution brain electromagnetic tomography (LORETA). Although the methods based on L2 norm are computationally efficient, their estimated sources are diffuse, covering most of the cerebral cortex outside the actual cortical activity area, and the result is that many sources cannot be spatially separated.

[0004] Accurate estimation of source range is necessary for studying brain activity. Estimating the brain tissue causing epilepsy and mapping brain function areas are important applications. In order to achieve accurate source range estimation, the present application proposes a nearest neighbor multivariate linear regression algorithm to achieve source range estimation. A graph Laplacian spatial smoother is used to model the source covariance, which is constantly updated using convex optimization theory. SUMMARY

[0005] The technical problems solved by the present application are that the traditional distributed source imaging method cannot accurately locate the size of the source range.

[0006] The method for locating the brain magnetic source based on the nearest neighbor multivariate linear regression model comprises the following steps:

[0007] Step 1: In a good magnetic shielding environment, first obtain the empty room signal data before the experiment starts to calculate the noise covariance matrix Σ ε , and then carry out the experiment to obtain high-quality brain magnetic signal data B and pre-process the obtained brain magnetic signal data B.

[0008] Step 2: First, use the nuclear magnetic resonance instrument to collect the MRI data of the subject. Segment the MRI image data to segment the subject's head into: scalp, skull, cerebrospinal fluid, gray matter and white matter, etc. Then create a single-shell model head model, and then register the head model coordinate system and the brain magnetic map sensor array coordinate system and further calculate the guide field matrix L.

[0009] Step 3: input the pre-processed brain magnetic signal data B, noise covariance matrix Σ ε and guide field matrix L into the nearest neighbor multivariate linear regression model for source positioning.

[0010] Step 3.1: first initialize the model parameters, obtain the graph Laplacian space smoother G according to formula 4 using the subject's MRI data, and then substitute the brain magnetic signal B, noise covariance matrix Σ ε and guide field matrix L obtained from step 1 and step 2 into the nearest neighbor multivariate linear regression model.

[0011] Step 3.2: first update the source signal estimate S using formula 5, then update the parameters a and γ using formulas 7 and 8 respectively, and then calculate the size of the loss function according to formula 6.

[0012] Step 3.3: judge whether the size of the loss function calculated in step 3.2 converges. If it does not converge, repeat step 3.2, and if it converges, the estimated value S of the source signal is the final source estimate.

[0013] The source positioning algorithm implemented by the present application can obtain more accurate source estimates, including accurate positioning of the source position, range size and source time waveform estimates, and provides a solid technical foundation for the medical application of brain magnetic. BRIEF DESCRIPTION OF DRAWINGS

[0014] Figure 1 Magnetic source localization flow chart.

[0015] Figure 2 Single-shell model head model.

[0016] Figure 3 Magnetic source localization principle flow chart based on nearest neighbor multivariate linear regression model. DETAILED DESCRIPTION

[0017] The present application proposes a kind of magnetic source localization method based on nearest neighbor multivariate linear regression model, the following algorithm principle and implementation steps are specifically explained by combining with Figure 1 、 Figure 2 And Figure 3

[0018] Before introducing specific operation steps, the relationship between the brain magnetic signal data B, the guiding field L derived from head model, source localization result S and noise ε is introduced briefly as follows:

[0019] B=LS+ε (1)

[0020] Wherein, the brain magnetic signal represented by acquisition is N c Represent channel number, N t Table time adopts the number of points, source signal is represented as N s It is assumed that the number of sources in distributed source model is represented, guiding field Mathematical representation of the magnetic source activity generated by neuronal activity under real head model conduction rule, noise is represented as

[0021] The specific steps of the present application are shown as follows Figure 1 The specific steps of the present application are shown as follows

[0022] Step 1: in the magnetic shielding room or magnetic shielding barrel and other magnetic shielding environment, use sensor (such as superconducting quantum interference device (SQUID) and optical pumping atomic magnetometer (QuSpin)) to measure empty room signal and brain magnetic signal data B.Before experiment, empty room signal needs to be measured to calculate noise covariance matrix Σ ε , carry out brain magnetic experiment again, obtain good quality brain magnetic signal data B, carry out some necessary pretreatment to the obtained brain magnetic signal data B, including identifying and rejecting bad channel and using band-pass filter and superposition average to remove noise.After the above pretreatment, brain magnetic signal data can be used for subsequent source localization.

[0023] ​Step 2: First, the MRI data of the subject is collected using a nuclear magnetic resonance instrument. Then the MRI data of the subject is imported into the Freesurfer software, and the MRI image data is segmented using the software to divide the subject's head into: scalp, skull, cerebrospinal fluid, gray matter and white matter, etc. Then a single-shell model head model is created using Matlab software, as shown in Figure 2 , only the conductivity of the cortex needs to be set. After building a suitable head model, the head model coordinate system and the sensor array coordinate system are registered and the guide field matrix L is further calculated.

[0024] Step 3: The preprocessed magnetoencephalographic signal data B, noise covariance matrix ε and guide field matrix L are input into the nearest neighbor multivariate linear regression model for source localization, and the specific process of the algorithm is as shown in Figure 3 .

[0025] Before introducing the specific operation step 3, the framework and principle of the nearest neighbor multivariate linear regression model are introduced:

[0026] First, the concept of nearest neighbor is introduced, as shown in Figure 2 . The cortex part in the subject's MRI data can be segmented by triangular mesh to obtain a series of point cloud data for easy reconstruction. When MRI is segmented, sources with close distance and in the same tissue are generally connected to form a connection relationship, and the concept of nearest neighbor is derived therefrom. In Figure 2 , the left is the segmented triangular mesh, and the right is a local enlarged view showing the central blue source and its surrounding 7 green sources with a connection relationship and close distance, which we call the nearest neighbor of the blue source.

[0027] There is a temporal and spatial relationship between the activities of neurons in the brain, so the most intuitive is that there is a temporal and spatial relationship between the source and its nearest neighbor. The present application uses a nearest neighbor multivariate linear regression model to describe this temporal and spatial relationship, as follows:

[0028]

[0029] where λ is a constant between 0 and 1, the present application takes 0.95, a n characterizes the different contributions of the source itself and its nearest neighbors at the previous time, and the common distributed source model contains nearly ten thousand sources, so a n also has nearly ten thousand, which will lead to overfitting of the model, so it is simplified to a, i.e. the a n corresponding to different sources are the same, N(n) represents the nearest neighbor set of the nth source, d n,i is a proportional coefficient representing the influence degree of different nearest neighbors on the source activity, generally d n,i = 1 / Nn N n The number of nearest neighbors, w n,t It is an unknown disturbance used to correct model bias. It is generally considered to be the same in all sources and is simplified to w. t The above formula can be simplified to:

[0030]

[0031] In this matrix, F represents the spatiotemporal relationship between sources, I is the identity matrix, Q is the weighted adjacency matrix, which represents the correlation between each source and its nearest neighbor, M is the adjacency matrix, and diag(h) represents the operation of transforming vector h into a diagonal matrix, where the diagonal elements are the values ​​of vector h. To better solve for the source localization results, we assume that a follows a Gaussian distribution with mean μ and variance β, and the unknown perturbation w t It also follows a Gaussian distribution with a mean of 0 and a variance of Σ. w The most common assumption is the variance Σ. w It is a diagonal matrix, and to achieve better source range estimation, this invention uses graphical Laplace to construct the variance Σ. w The details are as follows:

[0032]

[0033] Where G is the Graph Laplacian space smoother. N s This represents the number of source dipoles. We set the variance matrix. γ=diag(γ i ), where we assume γ i Independent of each other, G i It is the i-th column of G.

[0034] With the unknown perturbation in place, we can iteratively solve for the source estimate using the following formula:

[0035]

[0036] in, This represents the known source estimate from the previous time step; therefore, the source at time step 0 needs to be known in advance for the iteration to continue. The value of N t This represents the number of sampling points. However, F and Σ in the above formula... w It is not known, but F and Σ w The loss function is completely determined by parameters a and γ, respectively, therefore we need to solve for the values ​​of parameters a and γ. To do this, we can obtain the loss function through Bayesian inference. as follows:

[0037]

[0038] where,∑ b = L∑ w L T +∑ ε , r = B - LFS -1 , |·| is the determinant symbol, log represents taking the logarithm, trace is the operation of taking the trace of a matrix, is a constant independent of parameters a and γ. We can get the expression of the iteration of parameter a by setting the derivative of the loss function to 0 as follows:

[0039]

[0040] Matrices and are transition matrices when calculating a and have no actual meaning, and are the t-th columns of matrices and respectively. Due to the influence of the logarithm and the determinant, we cannot get the explicit expression of the iteration of parameter γ by setting the derivative of the loss function to 0. Since the loss function is a convex function of parameter γ, we use convex analysis to solve this problem and get the expression of the iteration of parameter γ as follows:

[0041]

[0042] In the above formula, is the set of the i-th source and its nearest neighbors, (·) i,j represents the i-th row and j-th column element of a matrix, (S-FS -1 ) i is the i-th column of matrix S-FS-1, and ||·||2 represents the two-norm of a matrix.

[0043] Next, the specific implementation of step 3 in the present application is described: Figure 3

[0044] Step 3.1: First, initialize parameters a = 0.5 and γ i = 1, then set the source s0 at time 0 as a 0 vector, use the MRI data of the subject to obtain the graph Laplacian space smoother G according to formula 4, and then substitute the brain magnetic signal B, the noise covariance matrix ∑ ε and the guide field matrix L obtained from steps 1 and 2 into the nearest neighbor multivariate linear regression model.

[0045] Step 3.2: First, update the source signal estimate S using formula 5, then update parameters a and γ using formulas 7 and 8 respectively, and then calculate the size of the loss function according to formula 6.​

[0046] Step 3.3: judging whether the loss function calculated in step 3.2 converges or not. If not, repeat step 3.2, if yes, the estimated value S of the source signal is the final estimated value of the source signal.

[0047] The contents not described in detail in the specification of the present application belong to the prior art known to the person skilled in the art.

[0048] The above embodiments are provided only for the purpose of describing the present application, and are not intended to limit the scope of the present application. The scope of the present application is defined by the appended claims. Various equivalent substitutions and modifications made without departing from the spirit and principles of the present application should be encompassed within the scope of the present application.​

Claims

1. A method of magnetoencephalography source localization based on nearest neighbor multivariate linear regression model, characterized in that, The method comprises the following steps: Step 1: calculate the noise covariance matrix using the empty room signal data before the start of the magnetoencephalography experiment Obtain magnetoencephalography signal data b after the start of the magnetoencephalography experiment, and pre-process the data; Step 2: Collecting MRI data of the subject, constructing a single-shell head model by using the MRI data, registering the data in the head model coordinate system and the sensor array coordinate system into the same coordinate system, and performing registration and further calculation to obtain a guide field matrix L; Step 3: Inputting the preprocessed magnetoencephalogram data, noise covariance matrix and guide field matrix into a nearest neighbor multivariate linear regression model for source localization; The step 3 comprises: Step 3.1: Initialization of parameters and Setting the source at time 0 a zero vector, the graph Laplacian space smoother is obtained from the MRI data according to equation (1), （1） where G is the graph Laplacian space smoother, , I is the identity matrix, M is the adjacency matrix, The brain magnetic signal data b obtained from step 1 and step 2 respectively, noise covariance matrix And the guide field matrix L are substituted into the nearest neighbor multivariate linear regression model, Step 3.2: updating the source signal estimation value S by using formula (2), （2） wherein, is the value of the source S at time t, denotes the known last time step source estimate, is the source covariance matrix assumed by the algorithm model, wherein it is assumed that are mutually independent, is the i-th column of G, is the brain magnetic sensor data at time t, F matrix characterizes the spatio-temporal connections between sources, denotes the number of sampling points; Reuse the formula (3) and (4) to update the parameters and , （3） and are parameters assumed by the algorithm model are the mean and variance of is a constant between 0 and 1, and are the t-th columns of matrices and respectively, trace is the operation of taking the trace of a matrix, and are transition matrices for calculating parameters and have no practical meaning, Q is a weighted adjacency matrix that can represent the correlation between each source and its nearest neighbor. （4） In the above formula is the set of the ith source and its nearest neighbors, denotes the element of the ith row, jth column of the matrix, is the matrix of the ith column, denotes taking the two-norm of the matrix, denotes the number of sampling points; Then the size of the loss function is calculated according to formula (5) ; （5） wherein denotes the number of sampling points, , , , is the determinant sign, log denotes the logarithm, trace is the operation of taking the trace of a matrix, is a constant parameter and are constants independent of both Step 3.3: Determine whether the loss function calculated in step 3.2 converges or not. If not, repeat step 3.

2. If yes, the estimated value S of the source signal is the final estimated value of the source. Step 3.3: Determine whether the loss function calculated in step 3.2 converges or not. If not, repeat step 3.

2. If yes, the estimated value S of the source signal is the final estimated value of the source.

2. The method of magnetoencephalographic source localization according to claim 1, wherein, The step 1 comprises: In a highly magnetically shielded environment, the brain magnetic signal data of the subject is measured by using a sensor, and the obtained brain magnetic signal data is preprocessed, wherein the preprocessing work comprises identifying and removing bad channels and removing noise.

3. The method of magnetoencephalographic source localization according to claim 2, wherein, The step 2 comprises: The MRI data of the subject is collected by using a nuclear magnetic resonance instrument, the MRI data of the subject is imported into a segmentation software for segmentation, the head cortex is segmented into scalp, skull, cerebrospinal fluid, gray matter and white matter, the conductivity of each cortex is set, a single-shell head model is constructed, the data in the head model coordinate system and the sensor array coordinate system are registered into the same coordinate system, and further calculation is performed to obtain a guide field matrix.

4. The method of magnetoencephalographic source localization of claim 1, wherein, The step 3.1 comprises: In the nearest neighbor multivariate linear regression model, the nearest neighbor multivariate linear regression model based on formula (6) is used to construct a brain model and represent the spatiotemporal correlation between brain source activities; （6） wherein, denotes the source activity at time t-1, denotes the nearest neighbor of the source activity at time t-1; denotes the unknown disturbance, is a constant between 0 and 1, taken as 0.95, characterizes the contribution of the source itself at the previous time, denotes the nearest neighbor set of the nth source, is a proportionality coefficient, characterizing the influence degree of different nearest neighbors on the source activity, taken as , denotes the number of nearest neighbors, is the unknown disturbance, used to correct the model bias; The source covariance matrix in the model constructed based on the graph Laplacian smoother in formula (1) is used to represent the spatial correlation of the covariance hierarchy between source activities; Based on the Bayesian inference and convex optimization theory, reasonable assumptions are made for the model, and the above formulas (3) and (4) for iterative solution of the model parameters are derived.

Images