An Automatic Identiﬁcation Method for the Blink Artifacts in the Magnetoencephalography with Machine Learning

: Magnetoencephalography (MEG) detects very weak magnetic ﬁelds originating from the neurons so as to study human brain functions. The original detected MEG data always include interference generated by blinks, which can be called blink artifacts. Blink artifacts could cover the MEG signal we are interested in, and therefore need to be removed. Commonly used artifact cleaning algorithms are signal space projection (SSP) and independent component analysis (ICA). These algorithms need to locate the blink artifacts, which is typically done with the identiﬁcation of the blink signals in the electrooculogram (EOG). The EOG needs to be measured by electrodes placed near the eye. In this work, a new algorithm is proposed for automatic and on-the-ﬂy identiﬁcation of the blink artifacts from the original detected MEG data based on machine learning; speciﬁcally, the artiﬁcial neural network (ANN). Seven hundred and one blink artifacts contained in eight MEG signal data sets are harnessed to verify the effect of the proposed blink artifacts identiﬁcation algorithm. The results show that the method can recognize the blink artifacts from the original detected MEG data, providing a feasible MEG data-processing approach that can potentially be implemented automatically and simultaneously with MEG data measurement.


Introduction
Magnetoencephalography (MEG) is a technique employing sensitive sensors to detect the weak magnetic field signal generated by the neurons of the brain without invasion or damage [1][2][3]. The detected data can be used to analyze the vital information of the brain. Currently, MEG is not only applied in advanced brain function research, such as auditory, vision, sensation, and even recognition [4][5][6][7][8][9], but is also widely used in epilepsy surgery and the diagnoses of Parkinson's disease, psychosis, and other brain function diseases [10][11][12][13][14].
Since the magnetic signal emitted from brain activities is feeble, the magnetic shielding room (MSR) is needed during the MEG recordings. In addition, other methods are also applied to stabilize the magnetic field with different types of sensors [15][16][17][18][19][20]. Nevertheless, if the original detected signal is directly observed, little information can be obtained owing to the existence of various types of noise [2,21]. Among all kinds of noise, the biological signals, such as blinks, eye movements, and heartbeats, are unique and the most disturbing for the MEG [22]. Interferences caused by these biological signals are called artifacts, among which the blink artifacts are the most significant and must be removed. In order to remove blink artifacts and observe the weak signal triggered by the brain activities more precisely [23,24], the original data need to be processed. Signal space projection (SSP) [25][26][27][28] and independent component analysis (ICA) [29,30] are widely used to remove the blinks and other artifacts in the MEG. However, these two methods may need extra electrodes to measure the blink signals in the electrooculogram (EOG) and decide the occurrence time of the blink artifacts in the MEG for practical use. A professional staff is also needed to identify the blink artifacts. To make up for the shortcomings of these two algorithms, some decent automatic processes for artifacts rejection based on statistics have been proposed by many research groups [31][32][33][34][35]. Machine learning algorithms are also introduced into the MEG data processing [36,37].
Some algorithms [38,39] have been proposed to classify different artifacts in the MEG and EEG (electroencephalography) based on machine learning. Most of these algorithms require data preprocessing, such as ICA, to extract the feature vectors of the artifacts. The preprocessing essentially reduces the dimension of the original detected MEG signal and then recombines these signals. The dimension reduction distorts the spatial dimension of the detected signal, which not only results in the loss of information of the original signal, but also adds some redundant messages into the MEG signal segments that are not influenced by the blink artifacts. In this work, we use the artificial neural network (ANN) to analyze the blink artifacts from the original detected MEG signal, and make full use of the information of their spatial distribution features, so as to realize its automatic identification.
With the development of computer science and information theory, the ANN is proposed as a powerful computing tool [40][41][42]. The ANN is a model built from the human brain neural network and can form different structures according to different connection methods, aiming to solve classification problems. It is a persuasive tool to extract data characteristics and eliminate noise, which plays a major role in various fields, such as automatic speech recognition (ASR) [43,44], computer image processing [45], and biomedical imaging [46,47].
The objective of the present study is to identify blink artifacts from the original detected MEG data. To address the aim, the original detected MEG data are color-coded into 2D images to display the spatial distribution of the blink artifacts, which are characteristic. These 2D images are taken as the input of the GoogLeNet [48] to automatically identify the blink artifacts. GoogLeNet is a convolutional neural network (CNN), which is a typical ANN and can recognize image features, indicating that no extra electrodes and professional confirmation are needed for artifacts identification and providing a novel perspective for the analysis of artifacts in the MEG. The total solution is shown to have the potential to assist the SSP and the ICA algorithm to realize the automatic and on-the-fly rejection for artifacts.

Materials and Methods
Data processing for the original detected MEG signal aims to remove noise, including line frequency noise, blink artifacts, eye movement artifacts, cardiac artifacts, etc. Blink artifacts are the noise that has the greatest influence on the original detected MEG signal. In order to remove the blink artifacts with the SSP and the ICA algorithm, the occurrence time of the blink artifacts should firstly be decided from the blink signals in the EOG. In this work, an ANN is utilized to recognize the spatial distribution characteristics of the blink artifacts, providing a new solution to automatically identify the blink artifacts from the original detected MEG data in real time. The process includes three sequential steps: obtaining the data, training the network, and testing the network.

Data Collection
The eight data sets used in this work come from Brainstorm website [49]. Four data sets are used to train the ANN, and all data sets are used to test the trained ANN. These data are recorded with a SQUID (superconducting quantum interference device) MEG system, which includes 274 sensors. The device is produced by the CTF corporation, Canada. The distribution of its sensors is shown in Figure 1 and the sensors involved in the data recording are given in Table 1.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 3 of 12 system, which includes 274 sensors. The device is produced by the CTF corporation, Canada. The distribution of its sensors is shown in Figure 1 and the sensors involved in the data recording are given in Table 1.  Eight MEG signal data sets are recorded by 274 magnetic field sensors. The 1st and 2nd data sets have a time duration of 360 s, the remaining six data sets have a time duration of 300 s, and all data sets have a sampling rate of 600 Hz. These data sets are collected from 4 different subjects and every two data sets are collected from the same subject, as is shown in Table 2. In total, 701 blink artifacts are contained in these data sets. It is a large enough data set for the artificial neural network (ANN) to extract spatial distribution characteristics.  Eight MEG signal data sets are recorded by 274 magnetic field sensors. The 1st and 2nd data sets have a time duration of 360 s, the remaining six data sets have a time duration of 300 s, and all data sets have a sampling rate of 600 Hz. These data sets are collected from 4 different subjects and every two data sets are collected from the same subject, as is shown in Table 2. In total, 701 blink artifacts are contained in these data sets. It is a large enough data set for the artificial neural network (ANN) to extract spatial distribution characteristics.
Filters are used to preprocess the original detected MEG signal data. A second-order infinite impulse response (IIR) notch filter with 3 dB bandwidth of 2 Hz is used to clear line frequency noise of 60, 120, and 180 Hz. An even-order linear phase finite impulse response (FIR) low-pass 200 Hz filter with 60 dB stopband attenuation is performed on the signal to remove other noise and leave the MEG signal. In addition, the blink signal in the EOG is employed to mark the blink artifacts in the MEG. In Figure 2a, according to the obvious blink signals provided by the vertical eye movement signal (VEOG), we can mark the time point of the emergence of the blink artifacts in the MEG signal, as is shown in Figure 2b. These time points allow us to accurately intercept the blink artifacts in the MEG as training data. In addition, when the ANN is used to identify the blink artifacts, we can also use these time points to compare the identification results. This step is often referred to as data calibration in machine learning. It should be noted that in Brainstorm, the blink signals in the EOG are used to determine the blink artifacts in the MEG. The specific algorithm is used to detect the rising edge of the blink signals to get its occurrence time. The detection algorithm needs to set fixed parameters including threshold, duration, frequency, etc. Nevertheless, because the blink signals are not always the same, coupled with the interference of other sources of noise, this method is not always effective. The new identification method can determine whether the sampled data are the blink artifact from the original detected MEG data. The center of all sampling time points identified as the blink artifacts is considered to be the time when each blink artifact appears, as is shown in the shadow of Figure 2. Table 2. The training data set and the testing data set can be obtained from the 8 MEG signal data sets. Every two MEG signal data sets are collected from the same subject. The 1st, 3rd, 5th, and 7th data sets are used to generate 2D images for training and testing. The 2nd, 4th, 6th, and 8th data sets are used to generate 2D images for testing only.

Training Data Set
Testing Data Set

Subject 1
The 1st data set 1200 blink 2D images 2000 non-blink 2D images 3600 2D images The 2nd data set For testing only 3600 2D images

Subject 2
The 3rd data set 2820 blink 2D images 3000 non-blink 2D images 3000 2D images The 4th data set For testing only 3000 2D images

Subject 3
The 5th data set 3900 blink 2D images 4000 non-blink 2D images 3000 2D images The 6th data set For testing only 3000 2D images

Training Data Set
In order to directly identify the blink artifacts from the original detected MEG data, we need to train the ANN with the calibrated MEG signal first. In this work, the Goog-LeNet is utilized to recognize the spatial distribution characteristics of the blink artifacts, which are highlighted by drawing 2D images. These 2D images are our training data.
Since each sampling data point is recorded by 274 sensors located at different locations on the head surface, it comprises a matrix with dimensions of 274×1. Each sampling

Training Data Set
In order to directly identify the blink artifacts from the original detected MEG data, we need to train the ANN with the calibrated MEG signal first. In this work, the GoogLeNet is utilized to recognize the spatial distribution characteristics of the blink artifacts, which are highlighted by drawing 2D images. These 2D images are our training data.
Since each sampling data point is recorded by 274 sensors located at different locations on the head surface, it comprises a matrix with dimensions of 274 × 1. Each sampling data point can be used to draw a 3D image of the signal strength at different locations on the head surface, as is shown in Figure 3a. One second of MEG data can be used to draw 600 images. However, 3D images are difficult for the GoogLeNet to recognize. The locations of the sensors distributed on a 3D surface can be topographically mapped onto a 2D plane. The signal amplitude measured with each sensor is subsequently color-coded as 2D images, as is shown in Figure 3b. These 2D images can be the input of the GoogLeNet. According to the calibrated MEG signal, a total of 66,780 2D images can be obtained with 8 data sets for training and testing of the ANN. The 1st, 3rd, 5th, and 7th data sets are employed to draw the training and validation 2D images. All data sets are employed to draw the testing 2D images.

Training Data Set
In order to directly identify the blink artifacts from the original detected MEG data, we need to train the ANN with the calibrated MEG signal first. In this work, the Goog-LeNet is utilized to recognize the spatial distribution characteristics of the blink artifacts, which are highlighted by drawing 2D images. These 2D images are our training data.
Since each sampling data point is recorded by 274 sensors located at different locations on the head surface, it comprises a matrix with dimensions of 274×1. Each sampling data point can be used to draw a 3D image of the signal strength at different locations on the head surface, as is shown in Figure 3a. One second of MEG data can be used to draw 600 images. However, 3D images are difficult for the GoogLeNet to recognize. The locations of the sensors distributed on a 3D surface can be topographically mapped onto a 2D plane. The signal amplitude measured with each sensor is subsequently color-coded as 2D images, as is shown in Figure 3b. These 2D images can be the input of the GoogLeNet. According to the calibrated MEG signal, a total of 66,780 2D images can be obtained with 8 data sets for training and testing of the ANN. The 1st, 3rd, 5th, and 7th data sets are employed to draw the training and validation 2D images. All data sets are employed to draw the testing 2D images. The 1st data set contains 20 blink artifacts, each of which lasts for approximately 0.4 s. In order to obtain the spatial information of the blink artifacts and expand the training data set, we capture a 0.1 s signal segment of the center of the blink artifacts and label these data as blink data, which means that we have a length of 0.1 (seconds) × 20 (blinks) × 600 (sampling rate) blink data. Since it is recorded with 274 sensors, the training data has the form of a matrix with dimensions of 274 × 1200. For signal segments with no artifacts, we randomly select 2000 sampled data for training and label them as non-blink data. Thus, we get a data matrix with dimensions of 274 × 3200, which can be used to draw 3200 2D images to train the ANN. Four 2D images of blink data and four 2D images of non-blink data are randomly selected from 3200 2D images, as is shown in Figure 4. Similarly, the 3rd, 5th, and 7th data sets include 71, 65, and 211 blink artifacts, respectively. We can get 2820, 3900, and 12,660 blink 2D images from these 3 data sets, respectively. At the same time, 3000, 4000, and 12,000 non-blink 2D images can be randomly selected from these data sets, respectively. These 2D images are then employed to train the ANN. The training data set is shown in Table 2.
blink data are randomly selected from 3200 2D images, as is shown in Figure 4. Similarly, the 3rd, 5th, and 7th data sets include 71, 65, and 211 blink artifacts, respectively. We can get 2820, 3900, and 12,660 blink 2D images from these 3 data sets, respectively. At the same time, 3000, 4000, and 12,000 non-blink 2D images can be randomly selected from these data sets, respectively. These 2D images are then employed to train the ANN. The training data set is shown in Table 2.

Artificial Neural Network
After obtaining the training data, the pre-trained network GoogLeNet is adopted to recognize the blink artifacts. It is a 144-layer CNN. Each layer of the network architecture can be regarded as a filter. The initial layer is primarily used to identify common features of the image, such as blobs, edges, and colors. The subsequent layers focus on more specific features to classify categories. GoogLeNet is pre-trained to classify images into 1000 object categories. According to the experiments of EEG data processing using the ANN [50], we need to readjust and retrain GoogLeNet for blink artifacts recognition, and specific adjustment of its 3 layers is made. We use MATLAB [51] to carry out the data processing.
To prevent over-fitting, a dropout layer is used in the original net. The dropout layer randomly sets input elements to zero with a given probability, and the default value is 0.5. Because we have fewer categories, in order to prevent overfitting, we replace the final dropout layer in the network with a new dropout layer with the probability of 0.6. Convolutional layers of the network extract image features. These features are employed to classify the input image by the last connected layer and classification layer. These two

Artificial Neural Network
After obtaining the training data, the pre-trained network GoogLeNet is adopted to recognize the blink artifacts. It is a 144-layer CNN. Each layer of the network architecture can be regarded as a filter. The initial layer is primarily used to identify common features of the image, such as blobs, edges, and colors. The subsequent layers focus on more specific features to classify categories. GoogLeNet is pre-trained to classify images into 1000 object categories. According to the experiments of EEG data processing using the ANN [50], we need to readjust and retrain GoogLeNet for blink artifacts recognition, and specific adjustment of its 3 layers is made. We use MATLAB [51] to carry out the data processing.
To prevent over-fitting, a dropout layer is used in the original net. The dropout layer randomly sets input elements to zero with a given probability, and the default value is 0.5. Because we have fewer categories, in order to prevent overfitting, we replace the final dropout layer in the network with a new dropout layer with the probability of 0.6. Convolutional layers of the network extract image features. These features are employed to classify the input image by the last connected layer and classification layer. These two layers of the original GoogLeNet contain information on how to combine these features into class probabilities, a loss value, and predicted labels. To classify the blink and non-blink images, we replace these two layers with new layers adapted to the data. The last connected layer is replaced by a new fully-connected layer with the number of filters equal to the number of classes (blink and non-blink). The classification layer is replaced by a new one without class labels, and it is automatically set as the output class during the network training.
Finally, we set the initial learning rate as 0.0001, the epoch as 10, and set a random seed as the default value in MATLAB [52]. The learning rate determines the variation range of parameters in the ANN and is a decreasing function of the epoch. The epoch represents the number of times the ANN is trained with the same set of training data. We use 80% of the images in training data for training and the remainder for validation. The training data comes from the first data set, which is also adopted to evaluate the recognition efficiency. The results are described in the section Results.

Testing
All 8 data sets are used to verify the efficiency of the method for blink artifacts identification. These data sets are taken as the testing data source. During 360 s (300 s) of data collection, 60 sampling data can be obtained every 0.1 s. Only one piece of sampling data is randomly selected and made into 2D images for identification. We can obtain 3600 2D images from the 1st and 2nd data sets and 3000 2D images from the remaining 6 data sets to be the testing data, as is shown in Table 2. The 2D image of the testing data set is then input into the trained ANN to identify whether it is a blink artifact. The identification results can be plotted in time for comparison with the MEG and the EOG signal. Detailed identification results are described in the Supplementary Materials.
In addition, in order to evaluate the identification results, we perform statistical calculation on the final identification results. The detection accuracy of each data set is obtained. The accuracy can be calculated as follows: False alarm represents that there are no blink artifacts in the original MEG signal data, but it is identified as blink artifacts by ANN. Missing alarm represents that the blink artifacts exist in the original MEG signal data, but they are not identified by ANN.

Results
In the first data set, there are practically 20 blinks [49]. The ANN identifies 20 blinks in total, with an accuracy of 100%. For the second data set, that contains 20 blinks, the ANN identifies 23 blinks. There are three false alarms at 124.7, 149.3, and 300.6 s. Compared with the 100% accuracy for the first data set, it can be thought that the false alarms are caused by insufficient training data. On the one hand, the fact that it can identify all blinks implies that the method is effective. On the other hand, the false alarms indicate that it does not learn enough data features. The false alarms mainly result from the MEG signal which has similar characteristics to the blink artifacts. The non-blink 2D images we randomly selected may not be enough to cover all the MEG signal features, making it difficult for ANN to distinguish the MEG signals that similar to the blink artifacts. For the third and fourth data sets, that include 47 and 80 blink artifacts, there are one and five false alarms, respectively, and zero missing alarms. The fifth and sixth data sets contain 65 and 78 blink artifacts, respectively. The identification results show there are zero and three false alarms, respectively, and one missing alarm in each. At last, 211 and 180 blink artifacts are included in the seventh and eighth data sets, respectively. The ANN has 18 and 30 false detections and 4 and 11 missing detections, respectively, for these two data sets. The identification results for the eight data sets are summarized in Table 3.
In addition, test accuracies of the first, third, fifth, and seventh data sets tend to be 100%, because the testing data from these data sets are similar to the training data set. For the other four data sets, we get accuracies of 86.96%, 94.12%, 95.06%, and 80.48% for the second, fourth, sixth, and eighth data sets, respectively. The eighth data set has distinctly worse results than other data sets, because the eighth MEG signal data set, measured from the fourth subject, contains more environmental noise (please see Supplementary Materials for more details). Since we only filter the original MEG signal data, different measurement environments will have a great impact on the identification results.

Discussion
In this work, blink artifacts identification is achieved by the ANN. Eight data sets collected from four different subjects are used to verify the recognition effect of the ANN. The first, third, fifth, and seventh data sets are employed to train the ANN. All eight data sets are used to test the ANN. The accuracy of the identification results shows the method can recognize the blink artifacts in the MEG signal. Nevertheless, the false alarms and the missing alarms appear. The missing alarms mainly result from insufficient training data. The MEG signals are signals with high inter-individual variability. Different subjects may produce different MEG signals. If we used the MEG data from one subject to train the network, and then used the trained network to identify the blink artifacts measured from another subject, the accuracy would drop dramatically. Considering that, using the ANN to recognize the blink artifacts generated by different subjects still needs further research. In addition, due to the different states of the subject during the experiment and the different measuring environment, the blink signals produced by the same subject may vary. The false alarms may be bound to the characteristics of the MEG signals and not to the characteristics of the blink artifacts. Some MEG signals may exhibit similar characteristics to the blink artifacts, resulting in misjudgment by the ANN. If various training data sets could be used to train the network, the recognition accuracy might be further improved, which needs to be proved by further experiments.
When recognizing spatial distribution characteristics of blink artifacts using the ANN, the identification algorithm does not introduce additional electrodes, and it can be implemented automatically and simultaneously with data measurement. These advantages endow the method with the ability, when combined with the SSP and the ICA algorithm, to realize the automatic and on-the-fly rejection for blink artifacts. Nevertheless, volumes of calibrated data are required to train the ANN to improve its accuracy. In this work, 66,780 2D images generated from eight MEG signal data sets are used for training and testing GoogLeNet. If we want to apply this method to practical problems, various types of artifacts are needed to train the ANN, which are missing in this experiment and need to be further supplemented.
In addition, if other artifacts data, such as eye movement artifacts and cardiac artifacts, are collected and employed for feature extraction, the ANN has the potential to classify and identify these artifacts in the MEG. The eye movement and cardiac artifacts are regarded as noise in the MEG and need to be removed, but for some other typical applications, such as the research on relationships between the eye movements and human brain activities, our method provides thoughts to extract the information from these artifacts. In addition, the machine learning algorithm has the potential to be migrated and applied to signal source estimation in the future. In magnetic source imaging (MSI), a more suitable network could be employed to classify the signals emitted from different parts of the brain, to establish the map between the location of the signal sources and the detected signal. Furthermore, the ANN should be capable of making a preliminary diagnosis by identifying the characteristics of the MSI. If such work can be realized through machine learning, it would hopefully reduce the burden of doctors and improve the efficiency of diagnosis and treatment.

Conclusions
In this work, a method based on the ANN for the identification of the blink artifacts is proposed, and 66,780 2D images generated from eight MEG signal data sets are used to train and test the convolutional artificial neural network, in which GoogLeNet is used as a pre-trained input. The identification based on the ANN enables us to identify the blink artifacts directly from the original detected MEG signal, rather than the EOG detected by electrodes. This method has the potential to be combined with the traditional SSP and ICA algorithm to realize the automatic rejection for blink artifacts.
The spatial distribution characteristics of blink artifacts can be recognized by ANN in this method. MEG signals with similar characteristics, such as cardiac artifacts, auditory evoked magnetic fields and visual evoked magnetic fields, can also be identified by this