Explainable Machine Learning Methods for Classification of Brain States during Visual Perception

Abstract

The aim of this work is to find a good mathematical model for the classification of brain states during visual perception with a focus on the interpretability of the results. To achieve it, we use the deep learning models with different activation functions and optimization methods for their comparison and find the best model for the considered dataset of 31 EEG channels trials. To estimate the influence of different features on the classification process and make the method more interpretable, we use the SHAP library technique. We find that the best optimization method is Adagrad and the worst one is FTRL. In addition, we find that only Adagrad works well for both linear and tangent models. The results could be useful for EEG-based brain–computer interfaces (BCIs) in part for choosing the appropriate machine learning methods and features for the correct training of the BCI intelligent system.

1. Introduction

At present, the application of artificial intelligence methods to the analysis of biological and medical data is an important and actively researched scientific task [1,2,3]. Among the works that are actively conducted in this direction are the analysis of medical images (CT and X-rays scans, etc.) [4,5], diagnostics of the human cardiovascular system [6], genomic medicine [7], different brain activity neurovisualization methods [8,9,10], etc. Within the latter, of particular interest is the diagnosis of various brain conditions and pathologies based on neuroimaging data such as a fMRI, fNIRS, MEG, and EEG [11,12,13]. In particular, machine learning methods including deep learning have found their application to process the brain signals in broad ranges, such as mental workload, disease prediction, stroke prediction [14,15], classification of EEG/MEG signals [16,17,18], prediction of sleep stages [19,20], and finding EEG biomarkers [21,22]. In the case of applying machine learning methods to diagnose EEG/MEG features in real time, there is an opportunity to implement brain–computer interfaces for neurorehabilitation, control of human brain states and robotics [23].

One of the central points of the application of artificial intelligence methods to medical and biological tasks is the interpretability of such approaches [24,25,26,27]. This is important for the creation of various assistive medical decision support systems [28,29], when a medical professional must understand and interpret the decision obtained using artificial intelligence methods. In this regard, in neuroscience, it is of great interest to develop and analyze various approaches for the diagnosis of neuroimaging data that are interpretable. Especially important is the interpretation of those features on the basis of which we create a particular machine learning system for application in biology and medicine. In the tasks of creating brain–computer interfaces for rehabilitation, communication with neurological patients, or the control of some external device, such as manipulators or exoskeletons, the most commonly used method is electroencephalography (EEG) to record brain activity and detect its features for subsequent command generation [23].

Explainable artificial intelligence (XAI) algorithms are considered to follow the principle of explainability [30]. A concept of explainability does not have a joint definition yet, so, there is a number of interpretations [31]. One of them is that explainability in machine learning can be considered as “the collection of features of the interpretable domain, that have contributed for a given example to produce a decision (e.g., classification or regression)” [32]. So, analyzing the influence of the different features on classification process is necessary for the explainable machine learning methods.

For an EEG-based brain–computer interface (BCI), deep mastering interpretability can display how different factors contained in EEG affect the machine learning model selection. For instance, Bang at al. [33] conducted analysis in comparison sample-sensible interpretation via the layer-wise relevance propagation (LRP) approach between the two subjects and revelead the potential reasons that cause the worse overall performance of one among all of them. The LRP approach is widely utilized by Sturn et al. [34] to analyze the deep learning model designed for a motor imagery mission. They assign the elements leading to an incorrect classification of artifacts of visible interest and eye movements, which remain in the EEG channels of the occipital and frontal regions. Ozdenizci at al. [35] proposed to apply a hostile inference approach to research strong features from EEG throughout unique subjects. Through interpreting the results with the LRP approach, they confirmed that their proposed method allowed the model of consciousness on neuroscience options in electroencephalogram, while it is less affected by artifacts from bone electrodes. Cui et al. [36] used magnificence’s class activation map (CAM) method [37] to examine character classifications of single-channel EEG alerts accumulated from a sustainably using enterprise. They identified that the model had been discovered to reveal brain activity patterns, such as alpha spindles and theta bursts, in addition to features that resulted from electromyography (EMG) activities, as proof to distinguish between drowsy and alert EEG alerts. In another work, Cui et al. [38] proposed a completely unique interpretation technique through taking gain of the hidden state output through the use of the long short-term memory (LSTM) layer to interpret the CNN-LSTM version designed for driving force drowsiness recognition from single-channel EEG. The same authors currently mentioned a novel interpretation method [39] based on a combination of the CAM approach [37] and the CNN–fixation techniques [40] for a multichannel EEG sign class and discovered stable functions across distinctive subjects for the venture of driver drowsiness popularity. Using the interpretation method, they also analyzed the reasons behind a few incorrectly classified samples. Regardless of the development, it is unclear to what extent the results of the interpretation may depend and how they may reflect the model’s selections. It is also not well defined in existing works why a specific interpretive technique is chosen above the other. These studies motivate us to conduct quantitative evaluations and comparisons of these interpretive strategies to gain in-depth knowledge of models designed for classification perceptual brain states of the group of voluntaries using EEG signals.

It should be noted that the convolutional neural networks (CNNs) being one of the most popular ML methods are commonly used and demonstrate a great success in image classification, natural language processing, computer vision, etc. [41,42]. However, their application to brain signals for achieving good accuracy requires a very deep CNN, which leads to a large number of parameters and high computational costs [43,44]. At the same time, the multilayer perceptrons are widely used for classifying EEG signals and demonstrate usually high efficiency [45,46].

Earlier in our works, various approaches based on machine learning have been investigated for EEG/MEG-based classification of brain states during the perception of different visual stimuli by the subject, including ambiguous images—Necker cubes [47,48,49]. In particular, the use of features based on the physiology of brain processes has been shown to improve the efficiency of the classification of brain states corresponding to the perception of visual images [50,51]. In our previous paper [52], the influence of image contrast (on the example of the famous painting by Leonardo da Vinci “Mona Lisa”) on information processing in the brain was investigated and the coherent resonance effect and reorganization of functional connections in the brain at a certain degree of image contrast were shown by the electroencephalography data. In the present work, we used this experimental dataset to classify brain states corresponding to different degrees of contrast. We also used a similar experimental dataset with EEG data on subjects’ perception of an ambiguous Necker cube with varying degrees of contrast to compare with perception of an artistic painting.

In this paper, we apply different machine-learning methods for the classification of brain states during visual perception and focus on the interpretability of the results. To estimate the influence of different features on the classification process and make the method more interpretable, we use the SHAP’s library technique. As data, we use 31-channel EEG recordings with a 250 Hz sampling rate filtered in five frequency bands. We find that a simple deep-learning model gives 100% for almost every dataset, which indicates an overfitting problem. So, we introduce four models with different combinations of two activation functions and different optimization methods. We find that the best optimization method is adagrad and the worst one is ftril. In addition, we find that only adagrad works well for both linear and tangent models. So, the contribution of the study is the following.

• We analyze the complex EEG dataset by using machine-learning techniques and find which optimization method is suitable for our dataset;
• Apply SHAP for estimation of the influence of different features to make the ML model more interpretable;
• Find the best optimizer that works well for both linear and tangent models.

The paper is organized as follows. In Section 2, we describe the EEG datasets used for the analysis, research paradigm, mathematical models of the deep learning approach (activation functions, optimization methods and the method for feature importance estimation), and case studies with the description of the neural network’s structure. Section 3 contains the results of application of different models for the classification of image intensity. In Section 4, we compare the considered deep learning methods with each other. Finally, in Section 5, we draw the conclusions.

2. Model and Methods

2.1. Datasets Description

Our computational analysis is based on the experimental neurophysiological data obtained earlier in the special experiments on the visual perception of images with different contrast levels [52]. The experiment consisted of observing images of varying brightness from $I=0.1$ to $I=1.0$, as shown in Figure 1. All the experimental EEG data of electric brain activity were recorded for 31 channels at a sampling rate of 250 Hz using the amplifier BE Plus LTM, manufactured by EB Neuro S.p.a., Florence, Italy. A detailed description of the experiment can be found in the ref. [52]. We used data from 5 subjects for two types of images. Two datasets corresponding to the observation of the ambiguous Necker cube (Figure 1A) and the Mona Lisa painting (Figure 1B) were considered. The complete set of all observed stimuli was formed by changing the brightness I of the stimuli, as shown in Figure 1A,B. A schematic illustration of the used experimental protocols is shown in Figure 2. At the start and at the end of the experiment, background activity was registered for 120 s. Each image with intensity I was presented to the subject for 60 s. Between the presentations, there was 20 s of rest. This study was conducted in accordance with the Helsinki Declaration and was approved by the Ethics Committee of Kant Baltic Federal University. All EEG data used for the analysis can be found in the repository [53].

2.2. EEG Preprocessing

Experimenters performed EEG preprocessing to remove various registration artifacts. After experimental registration, the EEG signals were filtered with a fourth-order Butterworth $(1÷100)$ Hz bandpass filter and a 50 Hz notch filter. In addition, an independent component analysis (ICA) was performed to remove eye blinking and heartbeat artifacts. It should be noted that in this study, we did not conduct any experimental studies and only used previously recorded data that had already been cleared of artifacts by the above procedures and with which we did not perform any additional manipulations.

2.3. Research Paradigm

The schematic representation of the research paradigm and overall structure of the research is presented in Figure 3.

We begin our consideration with the structure of the data to be analyzed. We consider trials of the 31-channel EEG with a duration of 1 min with a sampling rate of 250 Hz for each of 10 brightnesses I, i.e., at each moment of time n for brightness I, a data vector

$${\mathbf{x}}_n^I={\left(x_{1,n}^I,x_{2,n}^I,\dots ,x_{i,n}^I,\dots ,x_{31,n}^I\right)}^T$$

(1)

is registered. Here, index $i=1,2,\dots 31$ corresponds to EEG channels {Fp1, Fp2, …O1, O2} (see Figure 1C), $x_i$ is the signal registered in the i-th channel.

The most frequent and usual practice of EEG analysis is to look at different frequency ranges [54,55,56]. Here, we use the standard EEG delta-($\delta \in (1,4)$ Hz), theta-($\theta \in (5,8)$ Hz), alpha-($\alpha \in (8,12)$ Hz), beta-($\beta \in (13,30)$ Hz), and gamma-($\gamma \in (31,45)$ Hz) frequency ranges [57] for feature extraction and performance evaluation in perceptual brain states classification.

To obtain the data for each of the above frequency range, we reassigned the EEG signals (1) to the total average, subtracted the mean, and filtered with a fourth-order Butterworth $(f^L,f^H)$-Hz bandpass filter, where $f^L$ and $f^H$ are the boundaries of the frequency domain of interest [58]. For example, for the $\alpha$-range $f^L=8$ Hz and $f^H=12$ Hz, and we obtained the alpha-band signal ${\mathbf{x}}_n^{I,\alpha }$. Similarly, we obtain signals ${\mathbf{x}}_n^{I,\delta ,\beta ,\gamma }$ for all other $\delta$-, $\beta$-, and $\gamma$- frequency bands of interest.

So, we characterize the brain state during the visual perception of the images of each brightness by 60 s × 250 samples/s = 15,000 number of features (1) for each frequency band. We apply two deep learning models to classify brain states for different image brightnesses I. We consider three spatial domains of features: (i) all the EEG channels in the left and right hemispheres of the brain, (ii) EEG channels in the left hemisphere only, and (iii) EEG channels in the right hemisphere only.

We consider two separate strategies for learning and analyzing significant features that substantially affect model learning. In the first case, learning was based on the above-described features with separate frequency ranges, for each of which a different classification model was created. The input data were the trait values ${\mathbf{x}}_n^{I,\delta ,\alpha ,\beta ,\gamma }$ and image brightness I. The second case used a single model that combined all the features ${\mathbf{x}}_n^{I,\delta ,\alpha ,\beta ,\gamma }$ to predict image brightness I. In this case, different deep learning models with various neuron activation functions and various types of optimizers were used.

All deep learning models were tested on both datasets collected during the perception of images of Necker’s cube and Mona Lisa paintings. We applied SHAP (Shapley additive explanations) technology for feature importance estimation.

2.4. Mathematical Model of the Deep Learning Approach

Two types of activation functions are used: $tanh$ for the internal layers and $softmax$ for the output. For a compiler, we use $rmsprop$ as the optimizer and $categorical\_crossentropy$ for loss.

The $tanh$ function $\psi$ is described as

$$\psi \left(z_i\right)=\frac{1-e^{-2z_i}}{1+e^{-2z_i}}.$$

(2)

Here, $z_i=x_i$×$w_i+b_i$, where $x_i$ is the input layer value, $b_i$ is the broadcasting value through the columns (bias), $w_i$ is the weight value for $i\in [1,L]$, and L is the number of layers.

The standard $softmax$ function $\sigma$ is defined as

$$\sigma {\left(\mathbf{z}\right)}_i=\frac{e^{{\mathbf{z}}_i}}{{\sum }_{j=1}^K{e}^{{\mathbf{z}}_j}},$$

(3)

where $\mathbf{z}$ is the input vector, and K is the number of classes in the multi-class classifier. In our case, $K=10$ is the number of intensities I (Figure 1).

We also use rectified linear unit (ReLu) in the convolutions layers due to its computational simplicity and representational sparsity:

$$f\left(x\right)=\left\{\begin{array}{ccc}0\hfill & \mathrm{for}\hfill & x<0,\hfill \\ x\hfill & \mathrm{for}\hfill & x\ge 0,\hfill \end{array}\right.$$

(4)

where x is the input layer value.

Next, we will describe the optimization methods we use in the paper. Here, we use the following general notation: $\theta \in {\mathbb{R}}^d$ is the model parameters we need to optimize, $J\left(\theta \right)$ is the objective function, ${\nabla }_{\theta }\equiv g$ is the gradient of the objective function with respect to the parameters $\theta$, $\eta =0.01$ is the learning rate determining the size of the steps, t is the number of step, and $ϵ={10}^{-8}$ is a smoothing that avoids division by zero.

2.4.1. Stochastic Gradient Descent (SGD) as Optimizer

SGD is used for parameter updates for each training set. For example, if the dataset has $x^n$ features and $y^n$ target or label values, then for each step t, we have the following update rule [59]:

$${\Theta }_{t+1}={\Theta }_t-\eta g(\Theta ;x^{\left(i\right)};y^{\left(i\right)}),$$

(5)

where $i=1,\dots ,n$.

2.4.2. Root Mean Square Propagation (RMSprop) as an Optimizer

Root Mean Square propagation (RMSprop) [59,60] is very similar to gradient descent with momentum; the only difference is that it includes the second-order momentum instead of the first-order one, plus a slight change on the parameters’ update:

$$\begin{array}{cc}\hfill & E_t=\beta E_{t+1}+(1-\beta )g_t^2,\hfill \\ \hfill & {\Theta }_{t+1}={\Theta }_t-{\displaystyle \frac{\eta }{\sqrt{E_t+ϵ}}}{g}_t,\hfill \end{array}$$

(6)

where E is the decaying average over past squared gradients, and $\beta =0.9$ is the momentum term.

2.4.3. Adaptive Gradient Algorithm (Adagrad) as Optimizer

The adaptive gradient algorithm (Adagrad) [61] is a gradient-based optimization technique that achieves just that: it adjusts the learning rate to the parameters, producing more substantial updates for uncommon parameters and modest changes for frequent ones. As a result, it is ideal for dealing with sparse data. Dean et al. [62] discovered that Adagrad increased the robustness of SGD and used it to train large-scale neural networks at Google [63], which learned to detect cats in Youtube videos [64], among other things. In addition, Pennington et al. [65] employed Adagrad to train GloVe word embeddings, because uncommon words need considerably bigger updates than frequent words, as mentioned in Equation (7).

$$\begin{array}{cc}\hfill & {\Gamma }_t={\Gamma }_{t-1}+g_t^2,\hfill \\ \hfill & {\Theta }_{t+1}={\Theta }_t+\frac{\eta }{\sqrt{{\Gamma }_t+ϵ}}\odot g_t\hfill \end{array}$$

(7)

where ${\Gamma }_t$ is the sum of the squares of all past gradients, and ⊙ is the element-wise vector multiplication between $E_t$ and $g_t$.

2.4.4. Extension of Adagrad (Adadelta) as Optimizer

Adadelta [66] is an extension of Adagrad with a decreasing learning rate. The update rule of this methods is as follows [67]:

$$\begin{array}{cc}\hfill & \Delta {\theta }_t=-\frac{\sqrt{E{[\Delta {\theta }^2]}_{t-1}+ϵ}}{\sqrt{E{\left[g^2\right]}_t+ϵ}}{g}_t\hfill \\ \hfill & {\theta }_{t+1}={\theta }_t+\Delta {\theta }_t\hfill \end{array}$$

(8)

where $E{\left[g^2\right]}_t=\beta E{\left[g^2\right]}_{t-1}+(1-\beta )g_t^2$ and $E{[\Delta {\theta }^2]}_{t-1}=\beta E{[\Delta {\theta }^2]}_{t-2}+(1-\beta )g_{t-1}^2$.

2.4.5. Adaptive Moment Estimation (Adam) as Optimizer

Adaptive Moment Estimation (Adam) [68] is a method that can update parameters such as Adadelta and RMSprop. Here, the updated rule is as follows:

$${\theta }_{t+1}={\theta }_t-\frac{\eta }{\sqrt{{\widehat{v}}_t}+ϵ}{\widehat{m}}_t$$

(9)

where bias-corrected moments estimate:

$$\begin{array}{c}\hfill {\widehat{m}}_t=\frac{{\beta }_1{m}_{t-1}+(1-{\beta }_1)g_t}{1-{\beta }_1^t}\\ \hfill {\widehat{v}}_t=\frac{{\beta }_2{v}_{t-1}+(1-{\beta }_2)g_t^2}{1-{\beta }_2^t}\end{array}$$

(10)

where ${\beta }_1=0.9$ and ${\beta }_2=0.999$ are decay rates.

2.4.6. Extension to the Adaptive Movement Estimation (AdaMax) as Optimizer

The $v_t$ influence in the Adam update rule scales the gradient inversely proportionally to the ${\ell }_2$ norm of the past gradients (via the $v_{t-1}$ term) and current gradient ${\left|g_t\right|}^2$:

$$v_t={\beta }_2{v}_{t-1}+\left(1-{\beta }_2\right){\left|g_t\right|}^2.$$

(11)

Then, we update to the ${\ell }_p$ norm with parameterization ${\beta }_2$ as ${\beta }_2^p$ [68]:

$$v_t={\beta }_2^p{v}_{t-1}+\left(1-{\beta }_2^p\right){\left|g_t\right|}^p.$$

(12)

Norms for large p values generally become unstable. However, ${\ell }_{\infty }$ also generally shows balanced behavior. That is why the authors propose AdaMax [68] and show that $v_t$ with ${\ell }_{\infty }$ converges to the more stable value:

$$\begin{array}{cc}\hfill u_t& ={\beta }_2^{\infty }{v}_{t-1}+\left(1-{\beta }_2^{\infty }\right){\left|g_t\right|}^{\infty }\hfill \\ \hfill & =max\left({\beta }_2·v_{t-1},\left|g_t\right|\right).\hfill \end{array}$$

(13)

Then, we obtain the AdaMax update rule:

$${\theta }_{t+1}={\theta }_t-\frac{\eta }{u_t}{\widehat{m}}_t.$$

(14)

2.4.7. Nesterov-Accelerated Adaptive Moment Estimation (Nadam) as Optimizer

Nadam (Nesterov-accelerated Adaptive Moment Estimation) [69] combines Adam and Nesterov accelerated gradient. Combining Equations (9) and (10) and noting that $\frac{{\beta }_1{m}_{t-1}}{1-{\beta }_1^t}={\widehat{m}}_{t-1}$ which can be replaced by ${\widehat{m}}_t$, we obtain the Nadam update rule:

$${\theta }_{t+1}={\theta }_t-\frac{\eta }{\sqrt{{\widehat{v}}_t}+ϵ}({\beta }_1{\widehat{m}}_t+\frac{(1-{\beta }_1)g_t}{1-{\beta }_1^t}).$$

(15)

2.4.8. Follow the Regularized Leader (FTRL) as Optimizer

FTRL [70] strikes a trade-off between the advantages and disadvantages of forward–backward splitting (FBS) [71] and regularized dual averaging (RDA) [72]. The update rule is as follows:

$${\theta }^{(t+1)}={argmin}_{\theta }\left\{G^{(1:t)}·\theta +{\lambda }_1{\parallel \theta \parallel }_1+\frac{{\lambda }_2}{2}{\parallel \theta \parallel }_2^2+\frac{1}{2}\sum _{s=1}^t{\sigma }^{\left(s\right)}{∥\theta -{\theta }^{\left(s\right)}∥}_2^2\right\}$$

(16)

where ${\mathrm{G}}^{(1:\mathrm{t})}$ is the average of previous sub-gradients. In FTRL, the learning rate is different for different dimensions. If the training data deem that the dimension i needs to take a wider step than dimension j, the below learning rate will accommodate such updates:

$$\sum _{s=1}^t{\sigma }^{\left(s\right)}={\eta }_i^{\left(t\right)}=\alpha /\left(\beta +\sqrt{\sum _{s=1}^t{\left(g_i^{\left(s\right)}\right)}^2}\right).$$

(17)

Here, $\alpha$ and $\beta$ are a non-negative and non-decreasing sequence, and $\lambda$ is the learning rate.

2.5. SHAP (SHapley Additive exPlanations) for Feature Importance Estimation

Current explicable machine learning such as SHAP [73] supported the stochastic game [74] for a proof of model predictions. The Shapley value comes from cooperative theory of games, as an example, Shapley regression values [75] or Shapley sampling values [76]. Shapley regression values are feature importances for linear models within the presence of multiple correlation. The technique desires model fitting for all attribute subsets $S\subseteq F$. Every attribute ought to receive a gain that represents its contribution to the model prediction finally. To work out this gain for a given attribute $a\in F$, the predictions of the 2 models are compared on the present example x, wherever $x_S\cup \left\{a\right\}$ denotes the values of the input options within the set $S\cup \left\{a\right\}$. To take under consideration the influence of all alternative attributes, the variations are computed for each set $S\subseteq F\setminus \left\{i\right\}$. The ultimate Shapley price is calculated via a weighted average of these differences:

$${\varphi }_i=\sum _{S\subseteq F\setminus \left\{a\right\}}\frac{\left|S\right|!\left(\right|F|-|S|-1)!}{\left|F\right|!}\left[f_{S\cup \left\{a\right\}}\left(x_{S\cup \left\{a\right\}}\right)-f_S\left(x_S\right)\right].$$

(18)

2.6. Case Study

We implement two case studies shown in

Table 1. Description of case study.

	Participant Based	Validation	Focus Point	Model Number	Optimizer
Case Study I	Yes	No	Precision; Equation (19)	Single	Single
Case Study II	No	Yes	f1-score; Equation (21)	Multiple	Multiple

. In the first study, we obtained the maximum accuracy, which is 100% for almost every dataset that indicates an overfitting problem. So, we introduced a second method that has a data verification step and multiple model optimizers. Here, the equations for estimating the result:

$$Precision=\frac{TruePositive}{TruePositive+FalsePositive}$$

(19)

$$Recall=\frac{TruePositive}{TruePositive+FalseNegative}$$

(20)

$$F_1=2\times \frac{Precision\times Recall}{Precision+Recall}$$

(21)

2.6.1. Case Study I

First, we implement a simple model of ANN with RMSprop optimizer (Section 2.4.2) and a structure described in

Table 2. Deep learning model for Case Study I. The number of the parameters is equal to the number of connections between the layers.

Layer	Shape	Parameters	Activation Function
Input	15,000	-
Hidden layer 1	2500	37,502,500	tanh: Equation (2)
Hidden layer 2	1000	2,501,000
Hidden layer 3	500	500,500
Hidden layer 4	200	100,200
Hidden layer 5	100	20,100
Hidden layer 6	50	5050
Hidden layer 7	25	1275
Hidden layer 8	15	390
Output	10	176	softmax: Equation (3)

. The number of connections in the input layer corresponds to the number of features (15,000), and the number of output connections is equal to the number of brightnesses (10). The number of parameters is equal to the number of connections between the layers. For example, for the input and first hidden layer, it is defined as (15,000 + 1) × 2500 = 37,502,500, where “+1” is the additional channel for the intensity. The exception is the number of the parameters for the output layer: (15 + 1) × (10 + 1) = 176. Using this model, we obtain 100% accuracy for both datasets.

2.6.2. Case Study II

Next, we use the models with verification, different activation functions and optimizers. As validation, we default to the validation built in the Keras Sequential model: we choose 10% of the original dataset as the data on which to evaluate the loss and any model metrics at the end of each epoch. The model is not trained on these data. These data are only used for tuning hyper-parameters to make the model eligible for working well on unknown data.

We consider four models, which differ from each other in activation function: tanh, relu or their combinations. The structure of the models is described in

Table 3. Updated deep learning models for Case Study II. The ReLu activation function is described by Equation (4), tanh—Equation (2), softmax—Equation (3).

Layer	Shape	Param	Activation Function
			Model 1	Model 2	Model 3	Model 4
Input	15,000	-	tanh	relu	relu	tanh
Hidden layer 1	2500	37,502,500	tanh	relu	tanh	relu
Hidden layer 2	1000	2,501,000	tanh	relu	relu	tanh
Hidden layer 3	500	500,500	tanh	relu	tanh	relu
Hidden layer 4	200	100,200	tanh	relu	relu	tanh
Hidden layer 5	100	20,100	tanh	relu	tanh	relu
Hidden layer 6	50	5050	tanh	relu	relu	tanh
Hidden layer 7	25	1275	tanh	relu	tanh	relu
Hidden layer 8	15	390	tanh	relu	relu	tanh
Output	10	160	softmax

. The layers are the same as for Case I. So, the number of the parameters is also the same as for Case 1 but without an exception for the output layer: (15 + 1) × 10 = 160.

2.6.3. Computing System

The configuration of the computing system we used to perform ML follows:

• RAM: 503 GB;

• CPU: Intel(R) Xeon(R) CPU E5-2698 v4 @ 2.20 GHz;

• OS: Ubuntu 18.04.5 LTS 64 bit;

• GPU: NULL.

3. Results

In the paper, we use deep-learning neural networks with different structures, activation functions and optimizers for classification of the intensity (brightness I) of the image which was observed by the participant during the experiment. We have two datasets in which the Necker cube and Mona Lisa portrait were used as images, respectively. Each dataset includes 5 participants, 31 EEG channels, 10 intensities, and 15,000 data points for each of them. The dataset description is presented in more detail in Section 2.1.

We implement two case studies using different neural networks, which are described in Section 2.6. In both case studies, our goal is the classification of 10 intensities. Each of them is described below.

3.1. Case Study I

In this study, we decided to divide our dataset by frequency bands: $\delta$, $\theta$, $\alpha$, $\beta$, and $\gamma$ [77]. We divided our dataset into three parts based on EEG channels: (i) All channels, (ii) Channels from left hemisphere, (ii) Channels from right hemisphere, and we use the model shown in Table 2 on all of them [77]. The results of classification are shown in

Table 4. The results of classification by the deep learning model for Case Study I for Mona Lisa dataset when we use (first row) all EEG channels, (second row) the channels from only the left hemisphere of the brain and (third row) the channels from only the right one. The first value for each ”cell” (intersection of channels, participant and frequency band) is accuracy (* means 100%), and the second one is the intensity, which gives the maximum influence on the classification.

Channels	Participant	$\mathit{\delta }$	$\mathit{\theta }$	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\gamma }$
All EEG channels	1	*|0.1	*|0.1	*|0.1	98%|0.1	97%|0.1
	2	*|0.1	*|0.1	*|0.1	98%|0.1	95%|0.1
	3	*|0.1	98%|0.1	98%|0.1	94%|0.1>>	90%|0.1
	4	97%|0.1	97%|0.1	98%|0.1	98%|0.1	89%|0.1
	5	*|0.1	*|0.1	92%|0.1	92%|0.1	92%|0.1
Channels from left brain hemisphere	1	*|0.7	*|0.7	*|0.6	97%|0.5	97%|0.8
	2	*|0.7	*|0.7	*|0.7	*|0.9	*|0.7
	3	*|0.7	97%|0.9	97%|0.4	94%|0.9	94%|0.1
	4	*|0.7	*|0.2	94%|0.9	97%|0.9	94%|0.8
	5	*|0.7	*0.7	97%|0.9	*|0.7	*|0.9
Channels from right brain hemisphere	1	*|0.6	*|0.9	97%|0.6	*|0.6	90%|0.8
	2	*|0.6	*|0.6	*|0.6	*|0.6	93%|0.6
	3	*|0.6	*|0.6	*|0.6	*|0.6	90%|0.8
	4	*|0.6	*|0.6	*|0.6	97%|0.6	97%|0.6
	5	*|0.6	*|0.5	97%|0.5	97%|0.6	90%|0.6

^* Means 100% Accuracy.

and

Table 5. The results of classification by deep learning model for Case Study I for Necker cube dataset when we use (first row) all EEG channels, (second row) the channels from only the left hemisphere of the brain, and (third row) the channels from only the right one. The first value for each “cell” is accuracy (* means 100%), and the second one is the intensity, which gives the maximum influence on the classification.

Channels	Participant	$\mathit{\delta }$	$\mathit{\theta }$	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\gamma }$
All EEG channels	1	*|0.1	*|0.1	*|0.1	*|0.1	97%|0.1
	2	*|0.1	*|0.1	*|0.1	*|0.1	98%|0.1
	3	*|0.1	*|0.1	98%|0.1	98%|0.1	90%|0.1
	4	*|0.1	97%|0.2	97%|0.2	98%|0.1	*|0.1
	5	*|0.1	*|0.1	92%|0.2	*|0.1	97%|0.1
Channels from left brain hemisphere	1	*|0.7	*|0.7	*|0.7	*|0.7	91%|0.1
	2	*|0.7	*|0.8	*|0.3	*|0.9	*|0.8
	3	*|1.0	*|0.4	97%|0.9	97%|0.7	91%|0.5
	4	*|0.9	97%|0.4	*|0.9	97%|0.4	*|0.9
	5	*|0.7	*|0.7	*|0.9	*|0.9	*|0.9
Channels from right brain hemisphere	1	*|0.6	97%|0.2	97%|0.6	*|0.5	90%|0.6
	2	*|0.6	*|0.6	*|0.6	*|0.6	97%|0.6
	3	*|0.6	97%|0.6	*|0.6	97%|0.8	90%|0.6
	4	*|0.1	97%|0.6	*|0.6	*|0.6	*|0.6
	5	97%|0.6	*|0.6	97%|0.6	86%|0.2	93%|0.6

^* Means 100% Accuracy.

regarding the Mona Lisa and Necker cube datasets, respectively. The results of classification for both datasets are presented in Table 4 and Table 5. Here and further, the first value for each “cell” (intersection of channels, participant and frequency band) is accuracy defined by Equation (19) (* means 100%), and the second one is the intensity which gives the maximum influence on the classification. For estimating the influence, we use SHAP (SHapley Additive exPlanations) as described in Section 2.5. In tables, we show only the intensity with maximum influence, but SHAP allows estimating the influence of all intensities I for each time point. An example is shown in Figure A1 in Appendix A.

To estimate the quality of the DL model, we plot ROC, accuracy and loss curves. The ROC curves for all intensities are shown in Figure 4a. This plot is typical for every subset we investigated in Case Study I: AUC = 1.0 for every intensity, frequency range and dataset. Figure 5a illustrates the dependencies of accuracy and loss on every epoch. As one can see, the accuracy reaches 1.0 on the second epoch already, while the loss exponentially decreases with each increasing epoch.

Using this model, we obtained the maximum accuracy, which is 100% for almost every dataset that indicates an overfitting problem. So, we introduce the second method, which has a data verification step and multiple model optimizers.

The computation time for Case I is 155–162 s per model (1 dataset, 1 frequency) for all EEG channels and 74 to 82 s per model for the EEG channels from only the left or right hemisphere.

3.2. Case Study II

As it was mentioned above, using the first model for classification of the image’s intensity has an overfitting problem. To avoid it, we add data validation to our deep learning model. In addition, we estimate the f1-score (Equation (21)) instead of precision, as it was for model 1. The advantages of using the f1-score is that it can equally distribute all the classifier’s values, which cannot be done by using precision. Suppose for example that only 6 of 10 intensities were classified. Precision will work with only those 6 classified intensities, while the f1-score will use all 10 of them.

In Case Study II, we use four models (see Table 3) with different combinations of two activation functions (tanh and relu) to cross-check which model will work better (linear or tangent). We also use different optimizers to check how they perform for different models. The SHAP’s library technique [78] gives the proper explanation of the classification result of any deep learning model [77,79]. So, we can easily identify intensity with a maximum influence on classification (shown in Figure A1).

The ROC curves for all intensities are shown in Figure 4b for nadam optimizer. The curves are slightly different from the ones for Case I: not all AUCs are equal to 1.0, but all of them are higher then 0.99. As one can see, TPR is lower then 1.0 for low FPR for the number of classes, but these deviations are small and in a very short range, which indicates the high quality of the classifier. We should note that the illustrated ROC plot corresponds to the case with the highest deviation among all ROC plots. So, in every case, we have AUC ⪖ 0.99.

We estimate the accuracy and loss of the model during training and validation processes (Figure 5b,c). As one can see, the accuracy for both processes are rapidly growing and reaches its maximum value of 1.0 and 0.97, respectively, for the 10th epoch. Loss, in opposite, rapidly decreases and reaches zero value for the training process. However, for validation, it decreases up to 0.2 and then slowly increases.

Figure 6a represents the normal dependencies of them for each epoch, which are expected to be obtained. However, in some cases, they could differ from it: instead of the accuracy being increased, it decreases for a short epoch range and returns normal after it. An example of such irregular behavior is shown in Figure 6b.

Based on that, we define five behavior activities for estimation of the applicability of the model to the dataset. The description is provided in

Table 6. Behavior activity for estimation of the applicability of the model to the dataset for Case Study II. N means accuracy vs. epoch dependence corresponds to the normal one (Figure 6a), NN means irregular behavior (Figure 6b), ’-’ means accuracy is zero.

Model Accuracy Loss	Validation Accuracy Loss	Behavior
N	N	B1
NN	NN	B2
NN	N	B3
N	NN	B4
-	-	B5

. We compare the dependencies of accuracy and loss on the epoch to the normal ones (Figure 6a) during train and validation processes and obtain one of four results: during both processes, the dependencies correspond to the normal ones (B1), they both have irregular behavior (B2), or only one of them has irregular behavior (B3, B4). In addition, we introduce B5 behavior for the case of zero accuracy. So, B1 is the best behavior while B5 is the worst.

As it was mentioned above, in Case Study II, we use different models with different optimizers for image intensity classification by EEG signals filtered in five frequency bands. For each of them, we estimate behavior (Table 6), f1-score (Equation (21)), and the intensity with maximal influence. The computation time for Case II depends strongly on the activation function and optimization method; it varies from 288 to 6559 s per model. The average computational time for Case II is 2274 s. In the next subsections, we describe the results for each optimizer method.

3.2.1. Root Mean Square Propagation (rmsprop)

First, we use root mean square propagation (rmsprop) as the optimizer.

Table 7. The results of using the models with rmsprop optimizer for intensity classification.

Model	Dataset	Behavior	Frequency	f1-Score	Intensity
Model 1	Mona Lisa	B1	$\theta$	0.97	0.2
		B4	$\delta ,\alpha ,\beta ,\gamma$	0.97, 0.99, 0.95, 0.87	0.2, 0.4, 1, 0.9
	Necker Cube	B1	$\theta ,\alpha$	0.99, 0.97	0.2, 0.1
		B4	$\delta ,\beta ,\gamma$	0.99, 0.97, 0.93	0.2, 0.6, 0.1
Model 2	Mona Lisa	B1	$\alpha$ *	0.98	0.5
		B2	$\gamma$	0.27	0.8
		B5	$\delta ,\theta ,\beta$	-	-
	Necker Cube	B2	$\beta$	0.23	0.8
		B4	$\alpha$	0.25	0.9
		B5	$\delta ,\theta ,\gamma$	-	-
Model 3	Mona Lisa Necker Cube	B5	$\delta ,\theta ,\alpha ,\beta ,\gamma$	-	-
Model 4	Mona Lisa	B1	$\alpha$ *	0.99	0.1
		B4	$\delta ,\theta ,\beta ,\gamma$	0.97, 0.98, 0.95, 0.9	0.8, 0.5, 0.4, 0.2
	Necker Cube	B1	$\theta$	0.98	0.8
		B4	$\delta ,\alpha ,\beta ,\gamma$	0.99, 0.95, 0.95, 0.93	0.7, 0.1, 0.6, 0.7

* Refers to irregular behavior in accuracy vs. epoch graph shown in Figure 6.

shows the results of application of the models with rmsprop for classification. As one can see, models 2 and 3 do not perform well since they cannot classify the intensity for three and five frequency bands, respectively, for both datasets. Other models demonstrate high accuracy (f1-score > 0.9) and excellent (B1) and good (B4) behavior. Note that most frequencies for both datasets have the same behavior for models 1 and 4.

3.2.2. Adaptive Moment Estimation (Adam)

Next, we use Adaptive Moment Estimation (adam) as an optimizer.

Table 8. The results of using the models with Adam optimizer for intensity classification.

Model	Dataset	Behavior	Frequency	f1-Score	Intensity
Model 1	Mona Lisa	B4	$\delta ,\theta ,\alpha ,\beta ,\gamma$	0.99, 0.98, 0.96, 0.93, 0.85	0.7, 0.1, 1, 0.2, 1
	Necker Cube	B1	$\theta ,\gamma$	0.96, 0.86	0.5, 0.5
		B4	$\delta ,\alpha ,\beta$	0.98, 0.97, 0.9	0.2, 0.2, 0.6
Model 2	Mona Lisa	B2	$\beta$	0.32	0.3
		B5	$\delta ,\theta ,\alpha ,\gamma$	-	-
	Necker Cube	B5	$\delta ,\theta ,\alpha ,\beta ,\gamma$	-	-
Model 3	Mona Lisa Necker Cube	B5	$\delta ,\theta ,\alpha ,\beta ,\gamma$	-	-
Model 4	Mona Lisa	B1	$\theta$ *, $\alpha$ *	0.95, 0.97	0.6, 0.2
		B2	$\beta$	0.91	0.4
		B4	$\delta ,\gamma$	0.97, 0.9	0.7, 0.8
	Necker Cube	B1	$\theta$	0.96	0.1
		B4	$\delta ,\alpha ,\beta ,\gamma$	0.99, 0.97, 0.92, 0.83	0.6, 1, 0.1, 0.2

* Refers to irregular behavior in accuracy vs. epoch graph shown in Figure 6.

shows the results of application of the models with Adam for classification. As for the previous optimizer, models 2 and 3 do not perform well since they cannot classify intensity for four and five frequency bands, respectively, for both datasets. Other models demonstrate accuracy in the [0.83 0.99] range and excellent (B1) and good (B4) behavior. Note that the maximal f1-score is achieved for $\delta$ frequency in every dataset.

3.2.3. Nesterov-Accelerated Adaptive Moment Estimation (Nadam)

Then, we use nesterov-accelerated adaptive moment estimation (Nadam) as an optimizer.

Table 9. The results of using the models with Nadam optimizer for intensity classification.

Model	Dataset	Behavior	Frequency	f1-Score	Intensity
Model 1	Mona Lisa	B1	$\alpha$	0.96	0.1
		B4	$\delta ,\theta ,\beta ,\gamma$	0.99, 0.98, 0.93, 0.88	0.9, 0.2, 0.6, 1
	Necker Cube	B4	$\delta ,\theta ,\alpha ,\beta ,\gamma$	0.99, 0.97, 0.98, 0.94, 0.84	0.4, 0.4, 0.4, 0.7, 0.2
Model 2	Mona Lisa	B2	$\beta$	0.11	0.5
		B5	$\delta ,\theta ,\alpha ,\gamma$	-	-
	Necker Cube	B1	$\delta$	0.97	1
		B4	$\theta$	0.95	0.8
		B5	$\alpha ,\beta ,\gamma$	-	-
Model 3	Mona Lisa Necker Cube	B5	$\delta ,\theta ,\alpha ,\beta ,\gamma$	-	-
Model 4	Mona Lisa	B1	$\delta ,\alpha$	0.98, 0.97	0.2, 0.9
		B4	$\theta ,\beta ,\gamma$	0.99, 0.97, 0.86	0.4, 0.4, 0.1
	Necker Cube	B1	$\delta ,\alpha$	0.97, 0.94	0.4, 0.2
		B4	$\theta ,\beta ,\gamma$	0.97, 0.91, 0.88	0.3, 0.2, 0.3

shows the results of application the models with Nadam for classification. Again, models 2 and 3 do not perform well since they cannot classify intensity for four and five frequency bands, respectively, for both datasets. Other models demonstrate accuracy in [0.84 0.99] range as well as excellent (B1) and good (B4) behavior.

3.2.4. Extension to the Adaptive Movement Estimation (AdaMax)

The next optimizer we use is AdaMax, which is an extension to the Adaptive Movement Estimation.

Table 10. The results of using the models with AdaMax optimizer for intensity classification.

Model	Dataset	Behavior	Frequency	f1-Score	Intensity
Model 1	Mona Lisa	B1	$\delta ,\theta ,\alpha ,\beta ,\gamma$	0.99, 0.97, 0.98, 0.97, 0.90	0.2, 0.7, 0.3, 0.3, 0.1
	Necker Cube	B1	$\delta ,\theta$	0.98, 0.95	0.8, 0.6
		B4	$\alpha ,\beta ,\gamma$	0.96, 0.93, 0.91	0.4, 0.2, 0.2
Model 2	Mona Lisa	B2	$\theta$ *, $\alpha$ *, $\beta$	0.96, 0.8, 0.94	0.9, 0.6, 0.4
		B4	$\gamma$	0.83	0.8
		B5	$\delta$	-	-
	Necker Cube	B4	$\alpha ,\beta ,\gamma$	0.98, 0.93, 0.9	0.3, 0.2, 0.5
		B5	$\delta ,\theta$	-	-
Model 3	Mona Lisa	B1	$\theta$, $\alpha$ *, $\beta$ *	0.98, 0.96, 0.96	0.2, 0.6, 0.2
		B4	$\delta ,\gamma$	0.99, 0.88	0.7, 0.2
	Necker Cube	B1	$\delta ,\theta ,\alpha ,\beta$	0.98, 0.99, 0.95, 0.94	0.9, 0.8, 0.4, 0.4
		B4	$\gamma$	0.88	0.8
Model 4	Mona Lisa	B1	$\theta ,\alpha$	0.98, 0.97	0.9, 0.2
		B4	$\delta ,\beta ,\gamma$	1, 0.97, 0.89	0.5, 0.7, 0.6
	Necker Cube	B1	$\theta ,\alpha$	0.98, 0.96	0.8, 0.4
		B4	$\delta ,\beta ,\gamma$	0.96, 0.96, 0.91	0.7, 0.8, 0.2

* Refers to irregular behavior in accuracy vs. epoch graph shown in Figure 6.

shows the results of application of the models with AdaMax for classification. Comparing to all previous optimizers, AdaMax performs well for almost every frequency and model (except for model 2, which demonstrates B5 behavior for $\delta$ and $\theta$ frequencies). The accuracy for models 1, 3 and 4 are higher then 0.88. Moreover, model 1 for the Mona Lisa dataset demonstrates B1 behavior for all frequencies. The f1-score for model 4, the Mona Lisa dataset, and the $\delta$ frequency is 1.0. So, we need to introduce chaos [80] and dilution [81] to solve this problem for future work.

3.2.5. Adaptive Gradient Algorithm (adagrad)

Then, we use is adaptive gradient algorithm (adagrad) as optimizer.

Table 11. The results of using the models with adagrad optimizer for intensity classification.

Model	Dataset	Behavior	Frequency	f1-Score	Intensity
Model 1	Mona Lisa	B1	$\delta ,\theta ,\alpha ,\beta ,\gamma$	0.98, 0.97, 0.95, 0.91, 0.88	0.1, 0.9, 0.6, 1, 0.6
	Necker Cube			1, 0.96, 0.93, 0.85, 0.86	0.6, 0.5, 0.9, 0.4, 0.4
Model 2	Mona Lisa	B1	$\delta ,\theta ,\gamma$	0.98, 0.95, 0.83	0.2, 0.2, 0.9
		B4	$\alpha ,\beta$	0.94, 0.74	0.1, 0.7
	Necker Cube	B1	$\delta ,\theta ,\alpha$	0.96, 0.94, 0.92	0.2, 0.2, 0.1
		B4	$\beta ,\gamma$	0.75, 0.83	0.4, 0.2
Model 3	Mona Lisa	B1	$\delta ,\theta ,\alpha ,\beta ,\gamma$	0.98, 0.97, 0.96, 0.89, 0.86	0.5, 0.7, 0.3, 0.6, 0.9
	Necker Cube			1, 0.98, 0.95, 0.88, 0.88	0.5, 0.1, 1, 0.7, 1
Model 4	Mona Lisa	B1	$\delta ,\theta ,\alpha ,\beta ,\gamma$	0.98, 0.96, 0.96, 0.86, 0.86	0.5, 0.6, 0.6, 0.6, 0.3
	Necker Cube			0.99, 0.98, 0.95, 0.9, 0.84	0.3, 0.8, 0.1, 0.3, 0.9

shows the results of application the models with adagrad for classification. This algorithm does not demonstrate B5 behavior for any model. Moreover, all the models except model 2 demonstrate B1 behavior. One can see, that using lower frequency band gives higher f1-score. There are also 100% accuracy for $\delta$ frequency, Necker cube dataset, model 1 and 3. So, we need to introduce chaos [80] and dilution [81] to solve this problem for future work.

3.2.6. Extension of Adagrad (adadelta)

Next, we use the extension of Adagrad (Adadelta) as an optimizer.

Table 12. The results of using the models with the Adadelta optimizer for intensity classification.

Model	Dataset	Behavior	Frequency	f1-Score	Intensity
Model 1	Mona Lisa	B1	$\delta ,\theta ,\alpha ,\beta ,\gamma$	0.96, 0.4, 0.54, 0.2, 0.4	0.3, 0.5, 0.7, 0.8, 0.4
	Necker Cube			0.54, 0.51, 0.53, 0.3, 0.37	0.5, 0.4, 0.9, 0.3, 0.5
Model 2	Mona Lisa	B1	$\delta ,\theta ,\alpha ,\beta ,\gamma$	0.73, 0.59, 0.76, 0.66, 0.76	0.6, 0.4, 0.4, 0.4, 0.7
	Necker Cube			0.84, 0.68, 0.48, 0.53, 0.72	0.3, 0.1, 1, 0.8, 0.8
Model 3	Mona Lisa	B1	$\delta ,\theta ,\beta ,\gamma$	0.17, 0.25, 0.03, 0.15	0.8, 0.6, 0.3, 0.6
		B5	$\alpha$	-	-
	Necker Cube	B1	$\delta ,\theta ,\alpha ,\gamma$	0.07, 0.19, 0.15, 0.21	0.7, 1, 0.6, 1
		B5	$\beta$	-	-
Model 4	Mona Lisa	B1	$\delta ,\theta ,\alpha ,\beta ,\gamma$	0.14, 0.11, 0.1, 0.03, 0.18	0.6, 0.5, 0.9, 1, 0.1
	Necker Cube			0.32.0.11.0.14.0.03.0.18	0.7, 0.2, 0.9, 0.2, 0.1

shows the results of application of the models with Adagrad for classification. Although most models demonstrate B1 behavior, the f1-score is very low for many of them, and in some cases, it is close to zero. So, using the Adadelta optimizer does not perform good classification.

3.2.7. Stochastic Gradient Descent (SGD)

Then, we use stochastic gradient descent (SGD) as an optimizer.

Table 13. The results of using the models with an SGD optimizer for intensity classification.

Model	Dataset	Behavior	Frequency	f1-Score	Intensity
Model 1	Mona Lisa	B1	$\delta ,\theta ,\alpha ,\beta ,\gamma$	0.97, 0.99, 0.99, 0.97, 0.89	0.8, 0.8, 0.2, 0.1, 0.2
	Necker cube			0.98, 0.96, 0.96, 0.91, 0.88	0.2, 0.5, 0.4, 0.2, 0.8
Model 2	Mona Lisa	B1	$\delta$	0.73	0.7
		B2	$\gamma$	0.92	0.2
		B5	$\theta ,\alpha ,\beta$	-	-
	Necker cube	B2	$\beta ,\alpha$	0.19, 0.58	0.4, 0.1
		B5	$\delta ,\theta ,\gamma$	-	-
Model 3	Mona Lisa	B1	$\theta ,\alpha ,\beta ,\gamma$	0.96, 1, 0.96, 0.91	0.5, 0.9, 0.2, 1
		B2	$\delta$ *	0.25	0.1
	Necker cube	B1	$\theta ,\alpha ,\beta$	1, 0.99, 0.94	0.5, 0.4, 0.4
		B2	$\delta$	0.01	0.4
		B4	$\gamma$	0.94	0.6
Model 4	Mona Lisa	B1	$\delta ,\theta ,\alpha ,\beta ,\gamma$	0.98, 0.99, 0.98, 0.95, 0.9	0.7, 0.6, 0.5, 0.9, 1
	Necker cube			0.97, 0.98, 0.96, 0.93, 0.9	0.8, 0.1, 1, 0.7, 0.9

* refer as fault in epoch vs accuracy graph shown in Figure 6.

shows the results of application the models with SGD for classification. With this optimizer, models 1 and 4 (start and end with tanh) work perfect for both datasets because they demonstrate B1 behavior and high accuracy for all frequencies. Model 3 has good results for almost all frequencies except $\delta$ one. In contrast, model 2 demonstrates B5 behavior and does not perform well.

3.2.8. Follow The Regularized Leader (FTRL)

Finally, we use the follow the regularized leader (FTRL) algorithm for optimization. As one can see from

Table 14. The results of using the models with FTRL optimizer for intensity classification.

Model	Dataset	Behavior	Signal	f1-Score	Max Influence
Model 1, 3–4	Mona Lisa Necker Cube	B5	$\delta ,\theta ,\alpha ,\beta ,\gamma$	-	-
Model 2	Mona Lisa	B1	$\delta$	0.98	1
		B5	$\theta ,\alpha ,\beta ,\gamma$	-	-
	Necker cube	B1	$\delta$	0.98	0.7
		B5	$\theta ,\alpha ,\beta ,\gamma$	-	-

, it demonstrates B5 behavior for all the models and datasets. Only using the $\delta$ frequency band for model 2 gives good results.

4. Comparison of the Methods

4.1. Case Study I

As we know, right brain specialization is intuition and creativity while left brain specialization is language [82]. We found that using channels from different parts of the brain and filtering EEG signals in different frequency bands influences the classification process.

Table 15. The intensity which gives the maximal influence on the classification process for all participants in Case Study I.

Channel	Dataset	Frequency
		$\mathit{\delta }$	$\mathit{\theta }$	$\mathit{\alpha }$	$\mathit{\beta }$	$\mathit{\gamma }$
All EEG channels	Mona Lisa Necker cube	0.1	0.1	0.1	0.1	0.1
Channels from left brain hemisphere	Mona Lisa	0.7	0.7	0.9	0.9	0.8
	Necker cube	0.7	0.4	0.9	0.9|0.7	0.9
Channels from right brain hemisphere	Mona Lisa Necker cube	0.6	0.6	0.6	0.6	0.6

shows the intensity which gives the maximal influence on the classification process for all participants in Case Study I. As one can see, when we use all EEG channels or the ones from the right brain hemisphere, the same intensity gives the maximal influence for both datasets and all frequencies, while for left channels, the intensity differs for datasets and frequencies.

4.2. Case Study II

In Case Study II, we focus on difference keras [83] optimizers. As we mention in Table 6, we assign five behaviors. In this table, we omitted B5 behavior because it does not produce any output. The best behavior is B1, because both models have performed as a proper one. The next best one is the B3, because here, the validation performs good. The next best one will be B4, because its validation fails to perform well. So, the behavior activity from good to bad is $\mathrm{B}1>\mathrm{B}3>\mathrm{B}4>\mathrm{B}2$.

Table 16. Comparability within two datasets based on behavior and optimizer. Here, the cell value is the frequency bands whose f1-score is more than 70% and the number of the frequency corresponds to the behavior in brackets. B1/B2/B3/B4 are defined in Table 6.

Optimizer	Dataset	Behavior
		B1	B2	B3	B4
RMSprop	Mona Lisa	$\theta \left(1\right),\alpha \left(2\right)$	-	-	$\delta \left(2\right),\beta \left(2\right),\gamma \left(2\right),\alpha \left(1\right),\theta \left(1\right)$
	Necker cube	$\theta \left(2\right),\alpha \left(1\right)$	-	-	$\delta \left(2\right),\alpha \left(1\right),\beta \left(2\right),\gamma \left(1\right)$
Adam	Mona Lisa	$\theta \left(1\right),\alpha \left(1\right)$	$\beta \left(1\right)$	-	$\delta \left(2\right),\theta \left(1\right),\alpha \left(1\right),\beta \left(1\right),\gamma \left(2\right)$
	Necker cube	$\theta \left(2\right),\gamma \left(1\right)$	-	-	$\delta \left(1\right),\alpha \left(2\right),\beta \left(2\right),\gamma \left(1\right)$
Nadam	Mona Lisa	$\alpha \left(2\right),\delta \left(1\right)$	-		$\delta \left(1\right),\theta \left(2\right),\beta \left(2\right),\gamma \left(2\right)$
	Necker cube	$\delta \left(2\right),\alpha \left(1\right)$	-	-	$\delta \left(1\right),\theta \left(3\right),\alpha \left(1\right),\beta \left(2\right),\gamma \left(2\right)$
AdaMax	Mona Lisa	$\delta \left(1\right),\theta \left(3\right),\alpha \left(3\right),\beta \left(2\right),\gamma \left(1\right)$	$\theta \left(1\right),\alpha \left(1\right),\beta \left(1\right)$		$\delta \left(2\right),\beta \left(1\right),\gamma \left(3\right)$
	Necker cube	$\delta \left(2\right),\theta \left(3\right),\alpha \left(2\right),\beta \left(1\right)$	-	-	$\delta \left(1\right),\alpha \left(2\right),\beta \left(3\right),\gamma \left(4\right)$
Adagrad	Mona Lisa	$\delta \left(4\right),\theta \left(4\right),\alpha \left(3\right),\beta \left(3\right),\gamma \left(4\right)$	-	-	$\alpha \left(1\right),\beta \left(1\right)$
	Necker cube	$\delta \left(4\right),\theta \left(4\right),\alpha \left(4\right),\beta \left(3\right),\gamma \left(3\right)$	-	-	$\beta \left(1\right),\gamma \left(1\right)$
Adadelta	Mona Lisa	$\delta \left(2\right),\gamma \left(2\right)$	-	-	-
	Necker cube	$\delta \left(1\right),\gamma \left(1\right)$	-	-	-
SGD	Mona Lisa	$\delta \left(3\right),\theta \left(3\right),\alpha \left(3\right),\beta \left(3\right),\gamma \left(3\right)$	$\gamma \left(1\right)$	-	-
	Necker cube	$\delta \left(2\right),\theta \left(3\right),\alpha \left(3\right),\beta \left(3\right),\gamma \left(2\right)$	-	-	$\gamma \left(1\right)$
FTRL	Mona Lisa	$\delta \left(1\right)$	-	-	-
	Necker cube	$\delta \left(1\right)$	-	-	-

shows that the most frequent behavior of the RMSprop, Adam, Nadam, and AdaMax optimizers is B4, while for AdaMax, Adagrad, and SGD, it is B1. In our model, B1 is much better than B4. So, we can say that Adagrad is the best optimizer to use. Table 16 also gives information about which optimizer is suitable for which model. This table gives us a simple overview about linear and tangent models.

Another criterion we use for comparison of the model is the computational time (CT). Figure 7 illustrates the CT for all the models and optimizers which we have used in our research within Case Study II for $\delta$ frequency and the case of using all EEG channels. The lower the computational time of the machine learning method, the easier to use it. As one can see, the same optimizers marked by color require different times for the different models. The greatest difference was demonstrated by the Nadam optimization method (marked by brown color): more then 6000 s for models 1 and 3, 3500 s for model 2 and 2500 s for model 4; and RMSprop (pink): almost 4000 s for model 1 and around 2000 s for models 1, 3 and 4. The CT for Adadelta (blue), Adagrad (orange), AdaMax (red), and FTRL (violet) optimizers does not differ from one model to another. The CT for the Nadam optimizer is the highest for each model, so it is the most ineffective method among the others. In the opposite, the SGD optimization method (gray) requires the lowest CT for models 1–3, but for model 4, the most effective method is Adam (green).

To estimate the performance of the models, we measure precision, recall, f1-score, and specificity. Figure 8 illustrates the measures for all the models and optimizers. A good model should demonstrate precision, recall and f1-score values close to 1 and specificity value close to 0. As one can see, all optimizers except FTRL have 0 specificity and very high precision for model 1 and 3, and also all except FTRL and Adadelta have high recall and f1-score for those models. So, based on these measures, we can conclude that models 1 and 3 with Adagrad, Adam, AdaMax, Nadam, RMSprop, and SGD optimization functions are better than the others.

Table 17. Comparison with the existing studies.

Study	Year	Aim	Signals	Algorithm	Accuracy, %
El-Lone, R. et al. [84]	2015	Classification objects and animals by EEG	256 EEG channels	SVM	82.7
Spampinato, C. et al.	2017	Classification EEG data evoked by visual object stimuli	128 EEG channels	RNN-based model	84
Parekh, V. et al. [87]	2017	Image annotation system	14 EEG channels	EEG-Net	88
Zheng, X. et al. [86]	2020	Image classification by analyzing images and EEG signals	128 EEG channels + images	LSTMS-B	96.94
Kuc, A. et al. [51]	2021	Ambiguous stimuli classification	31 EEG channels	CNN	74
Current study	2022	Classification of the image’s intensity	31 EEG channels	MLP, adagrad optimizer	92.9

demonstrates a comparative study of methodologies and results between the current work and previous machine learning-based EEG studies aimed at classifying brain states In [84], the classical SVM classifier was used to classify EEG involving biological categories and non-biological categories, and it reached the classification accuracy of 82.7%. In [85], the RNN-based model was used for the classification of EEG data evoked by visual object stimuli, and 84% accuracy was achieved. In [86], the authors reached a great 96.94% f1-score, but we cannot it as the baseline since they used not only EEG signals for classification but also the original images. In our study, we use only EEG signals, and 31 channels is enough to achieve a 92.9% f1-score (average value for Model 1 over all frequencies) for the Adagrad optimizer. We should note the EEG-Net model presented in Ref. [87] where the authors used only 14 EEG channels to achieve 88% accuracy, which is a little bit less than our result. Overall, our model demonstrates a good ability to classify the EEG signals.

5. Conclusions

We have applied different machine-learning methods for the classification of brain states during visual perception. As data, we used 31 EEG channels filtered in $\delta$, $\theta$, $\alpha$, $\beta$, and $\gamma$ frequency bands corresponding to the perception of Mona Lisa and Necker cube images. In Case Study I, we used a deep learning model with eight layers, the tanh activation function, and the RMSprop optimizer. Using this model, we obtained the maximum accuracy which is 100% for almost every dataset that indicates an overfitting problem. To estimate the influence of different features on the classification process and make the method more interpretable, we use the SHAP’s library technique.

To avoid the overfitting problem of the first model, we introduced the second method (Case Study II), which has a data verification step. Here, we used four models with different combinations of two activation functions (tanh and relu) to cross-check which model works better (linear or tangent). We also used different optimizers to check how they perform for different models. We found that the best optimization method is Adagrad; it performs well for most of the frequency bands. In contrast, the FTRL method does not work at all. The list of optimizers from best to worst is: Adagrad > SGD > AdaMax > Adadelta > Nadam > RMSprop > Adam > FTRL. In addition, we found that only Adagrad works well for both linear and tangent models.