The autoregressive (AR) Burg approach is applied to estimate the absolute power spectrum density (PSD) for each channel (Faust et al., 2008). To reach the objective of power comparisons in four frequency bands and different experimental states, the EEG data are segmented into 1-s epochs with 0.5-s overlap. A power spectrum analysis is performed for each time window. At a certain subband, the final power value at a certain time is obtained by averaging all power values of the frequency band on the current time window interval. The theoretical basis of the AR method is that the given signal is the output of a linear system whose input is white noise, which is (1)r(n)=-∑i=1kair(n-i)+u(n) where a_i represents the AR coefficients, k is the order of the AR model, and u(n) is the white noise with a variance of σ². The order k can be described by AR parameters {a1,a2,…,an,σ2}. The PSD is defined as (2)P^AR(f)= σ2|1+∑1paie-j2πfi|2

Principal component analysis is the common linear dimensionality reduction technique (Duffy et al., 1992). In this research, PCA is applied to the EEG of all subjects. To avert the case where channels with excessive fluctuations impact decoding, the EEG data are smoothed through a Gaussian kernel before employing dimensionality reduction techniques. In this research, PCA is performed to the EEG data of a single-trial level rather than trial averaged from all subjects to avoid the inaccuracy caused by averaging the responses across entire trials. The purpose of PCA is to extract principal components (PCs) from smoothed EEG data. Meanwhile, all PCs are arranged in descending order according to the magnitude of explained variances. The first 2 PCs are selected to decode visual stimuli.

Compared with PCA, t-SNE is a non-linear dimensionality reduction technique. t-SNE calculates the probability similarity for points using a normalized Gaussian kernel in a high-dimensional space (Van der Maaten and Hinton, 2008). Correspondingly, in low-dimensional space, the similarity is calculated through t-distribution. The similarity between sample points in high-dimensional space can be expressed by the asymmetric point probability distribution p_ij. For the t-SNE algorithm, the p_ij is defined as follows: (3)pij=exp(-||ri-rj||22σi2)∑k∑l≠kexp(-||rk-rl||22σi2),∀i,j and i≠j (4)pij=pj|i+pi|j2N where σ_i is the variance of a Gaussian centered on r_i. Due to the different distribution density of sample points in the dataset, σ_i corresponding to different sample points is also different. The denser the distribution of data point is, the smaller σ_i is. The sparser the distribution of data point is, the larger σ_i is. The value of σ_i is calculated by binary search.

Takagi-Sugeno-Kang is a fuzzy model (Takagi and Sugeno, 1985). For an input dataset U={u1,u2,…,un}∈Rd and the corresponding class label O = {o₁, o₂, …, o_n} (oi=[oi,1,oi,2,…oi,j]T). When the ith sample data belong to jth class, o_i,j = 1; Otherwise, o_i,j = 0. The kth fuzzy inference rules are usually described as (5)Rk:IF u1 is A1k∧u2 is A2k∧…∧ud is Adk, THEN fk(u)=β0k+β1ku1+…+βdkud,k=1,…,K where the input vector of each fuzzy rule is represented as u=[u1,u2,…ud]T, Aik are Gaussian antecedent fuzzy sets subscribed by the input variable u_i of Rule k, the number of fuzzy rules is K, ∧ is a fuzzy conjunction operator, the kth fuzzy inference rule f_k(u) can be expressed by a linear combination of the input vector u=[u1,u2,…ud]T, and βik are linear coefficients.

If the input vector of each rule is u, the output of TSK is calculated as (6)õ=∑k=1Kμk(u)fk(u)∑k ′Kμk ′(u)=∑k=1Kμ~k(u)fk(u) where μk(u)=∏i=1dμAik(ui) is the fuzzy membership function and the normalized fuzzy membership function of the kth fuzzy rule is denoted as μ~k(u)=μk(u)∑k ′Kμk′(u). μAik(ui) is Gaussian membership function for fuzzy set Aik, which is defined as (7)μAik(ui)=exp[-(ui-cik)2δik] where the element cik is kth cluster center parameters that can be obtained through the classical fuzzy c-means (FCM) (Bezdek et al., 1984) clustering algorithm: (8)cik=∑j=1Nxjkuji∑j=1Nxjk where the width parameter δik can be estimated by: (9)δik=l·∑j=1Nujk(xjk-cik)2∑j=1Nujk where x_jk ∈ [0, 1] is the fuzzy membership between jth input vector and kth cluster (k = 1, 2, …, k), and the constant l is the scale parameter.

Given u_n as the input sample, we can obtain (10)un,e=(1,unT)T (11)ũnk=μ~k(un)ue (12)ρ(un)=[(ũn1)T,(ũn2)T,…,(ũnK)T]∈RK(d+1) (13)βk=(β0k,β1k,…,βdk,)T (14)βg=[(β1)T,(β2)T,…,(βK)T] Thereby, the output result õ_n of a TSK fuzzy decoder for the input vector u_n can be described as (15)õn=βgTρ(un)

The NB decoder is a type of Bayesian decoder that determines the probabilities of different outcomes, which can predict the most probable outcome. For the input dataset U={u1,u2,…,un}∈Rd with m classes and a sample C = {c₁,c₂, …, c_n}, Bayes' theorem can be defined as (16)P(ui|C)=P(C|ui)P(ci)P(C) Because P(C) is constant for all classes, the probability that C belongs to each class can be estimated by (17)P(ui|C)∝P(ci)∏k=1nP(ck|ui)=∏k=1nP(ck|ui)[P(ui)]n-1 According to this formula, the class with the highest probability is taken as the class of C.

An LSTM-RNN based on an RNN architecture has been well befitting in decoding brain activities (Tsiouris et al., 2018). The LSTM-RNN is constituted by memory cell c, forget gate f, input gate i, and output gate o. LSTM-RNN can accommodate the information flow through a cell via these gates, and each gate can be expressed by (18)ft=σsigmoid(Wf,zut+Wf,1ht-1+bf) (19)it=σsigmoid(Wi,zut+Wi,1ht-1+bi) (20)ot=σsigmoid(Wo,zut+Wo,1ht-1+bo) where b is the vector of biases, h denotes the vector of the hidden layer, and W represents the weight matrix of recurrent connection from the input gate to the output gate. The input vector u_t is the latent factors at time t and σ_sigmoid(·) represents the sigmoidal activation function. The subscripts of weight matrix W indicate the corresponding gates that include forgetting gate f, input gate i, and output gate o. The information flow of the cell memory can be updated by (21)cu= σtanh(Wc,xut+Wc,1ht-1+bc) (22)ct=ft⊗ct-1+it⊗cu (23)ht=ot⊗σttanh where σ_tanh and ⊗ represent the hyperbolic tangent function and the element-wise product, respectively.

Linear Kalman filter is a commonly generative decoding method based on the linear dynamical. The hidden state y_t is a linear combination of the hidden state y_t−1 (Wu et al., 2003). The system model can be described as (24)yt=Ayt-1+Qt-1 (25)lt=Hyt+Vt where A denotes the system model parameter matrix, H is the observation model parameter matrix, and Q and V are the process noise and observation noise that obeys a Gaussian distribution. The observation vector l_t is the vector of latent factors. After converging the Kalman gain to its stable state early, we initialized y₀ to the average of the observation vector at the beginning of every trial to predict y_t.

We first investigate whether there exist specific frequencies that are dominant in the digital cognitive process. EEG recordings of the awake volunteers are performed to explore the changes in EEG evoked by digital cognition (Figure 2A). The EEG signals are divided into four major frequency bands involved in sensory consciousness, information coding, cognitive memory, selective attention: delta (0.5–4 Hz), theta (4–8 Hz), alpha (8–15 Hz), beta (15–30 Hz) (Figure 2B). In a dark behavioral chamber, presentations of visual stimuli (30 s) induce alpha frequency oscillations, but the oscillations persist only ~10 s during the digital stimuli period (Figures 2C,G). The frequency spectra show that the power of the alpha band increases during the stimuli period compared with that during the prestimuli and poststimuli periods (Figure 2D). Furthermore, the spectral power across four frequency bands is compared in Figures 2G,H. It is found that alpha band power increases during the stimuli period with p < 0.01 whereas delta band power decreases with p < 0.05 (one-way ANOVA followed by Bonferroni's post-hoc test).

The previous study has shown that specific brain areas may participate in specific cognitive functions (Gatti et al., 2020). The spectral power of the three periods across four brain lobes is calculated (Figures 2E,F). Alpha band power increases with the varying level in different brain lobes during the stimuli period, especially in the frontal and the occipital lobes. The variation of power of each EEG channel data in alpha is further analyzed and it is obtained that channels with relatively higher power spectrum values are mainly concentrated in the frontal and the occipital lobes (Figure 2I). In the following dimensionality reduction and decoding analysis, we are mainly focusing on the alpha frequency band and the first 10-s EEG signals during visual tasks. These results demonstrate that the alpha frequency band is the dominant frequency. Meanwhile, the frontal and the occipital lobes might be closely related to digital cognitive function.

To identify the dynamics of brain states of different cognitive processes, dimensionally reduction techniques that include PCA and t-SNE are further leveraged to estimate latent factors that reflect the subject's brain. Linear dimensionality reduction based on PCA is first performed. It is the ineffective separation of brain dynamics at different cognitive task epochs by PCA (Figure 3A). Three states are not possible to be distinguished by chosen viewing angle. It describes the cumulative explained variance of brain states as a function of the number of PCs in Figure 2B. It is obtained that the brain states are high-dimensional since explaining 80% of the total variance requires 9 dimensions at least. It is difficult to capture key differences during different task epochs using PCA, and thus, dimensionality reduction fails to achieve effectively in a fashion. To overcome this issue, a non-linear dimensionality reduction technique, t-distributed stochastic neighbor embedding (t-SNE), is applied to the time series of brain activities.

The non-linear dimensionality reduction technique t-SNE is used on the single trial, and the distributions of latent factors are shown in Figure 4A. Since the t-SNE is based on manifold learning, the latent factors after dimensionality reduction are divided into different manifolds. Except for a very small number of latent factors, the latent factors of the same task can be gathered in a region and the t-SNE clusters coincide well with distinct tasks. There are obvious classification boundaries among different tasks, which provides good conditions for the following decoding. More importantly, clusters in the rest states including prestimuli and poststimuli states are located in the same area. Latent factors from 10 experiments are further analyzed in the same way. The similarity of differences is measured via the pairwise Euclidean distances between brain states in the t-SNE space. The similarity between the clusters of two different tasks (Pre-Num1, Num1-Post, Num1–Num5) is lower than the similarity between the same tasks (pre–post) in Figure 4B. Different visual stimuli drive different transitions among brain states resulting in a lower similarity of the brain states between two different visual stimuli (Num1–Num5).

To investigate the transition process of brain states, the consecutive brain states are temporally connected in the t-SNE map on the trial-by-trial level, which is shown in Figure 4C. It illustrates that the brain states dynamically evolve into cyclic manifolds which are distinct with different tasks over time, whereas the evolutive paths of brain states are highly similar during the same task. Afterward, to further analyze the difference among temporal evolutions of brain states during different cognitive tasks, the partial correlation coefficient is calculated between each brain state in a trial and a reference brain state determined by averaging all brain states of five tasks. The brain states are projected into 2D space spanned by reference brain states for the number of 1 and the number of 5. Strikingly, the bifurcation toward different cognitive tasks could be observed well at the beginning of visual stimuli at the single-trial level (Figure 4D), which further suggests that different visual cognitive tasks drive different transitions among brain states.

Given our analyses at the level of brain states and the latent factors under visual stimuli, it can be concluded that decoding brain activities by the non-linear dimensionality reduction technique is reliable. After dimension reduction, multiple machine learning decoders including linear decoder Kalman filter (KF), non-linear decoder long–short-term memory in recurrent neural networks (LSTM), fuzzy learning decoder Takagi-Sugeno-Kang fuzzy system (TSK), and naïve Bayes decoder (NB) are employed to decode brain activities. The decoding performance is evaluated for each combination of two-dimensionality reduction methods and four decoders (Figure 5). For the PCA method, the frequency band energy of the first two principal components is extracted, and the values x₁, x₂ are used as input vector x of the decoders. For the t-SNE method, the first two principal components x₁, x₂ are retained directly to obtain the neural manifold, which is used as the input vector x to the decoders. The decoding is implemented through 10-fold crossvalidation, which means that 9 out of 10 samples are selected as a training set, and the remaining 1 out of 10 samples are selected as the test set in each classification process. The classification process is repeated ten times, and the average value of the 10 times is taken as the final result. TSK decoder displays the encouraging consequence with 0.8198 accuracies with t-SNE and 0.6971 accuracies with PCA. The features extracted by the non-linear manifold learning method have better decoding performance than that of PCA. The performance of TSK is compared with three common methods (NB, LSTM, LKF) for EEG decoding. The accuracy data are normally distributed as the results of the Q-Q plot are approximately on a straight line (Figure 6). The classification results are shown in Table 1, Figure 7. Table 1 shows the classification accuracy of different decoders with two whole-brain features extracted by PCA and t-SNE. Meanwhile, TSK can achieve high performance in decoding EEG. The comparison of the decoding accuracy for visual stimuli across different latent factors and decoders is depicted in Figure 7. Friedman's test with Bonferroni correction is used to analyze the difference between the decoding accuracy of four decoders. The Friedman's test reveals the main effect of decoders on the decoder accuracy when using t-SNE (χ² = 20.97, p < 0.01) or when using PCA (χ² = 25.51, p < 0.01). When using t-SNE, the linear decoder KF shows lower accuracy than other decoders (p < 0.01) and TSK shows higher accuracy than NB (p < 0.05), whereas there is no difference between TSK and NB, LSTM. When using PCA, the TSK shows higher accuracy than other decoders (p < 0.01), and NB and KF also show higher accuracy than LSTM (p < 0.01). These indicate that TSK has a good decoding effect in both linear and non-linear dimensionality reduction techniques. Meanwhile, for all decoders, the t-SNE dimensionality reduction technique has a better decoding effect than PCA.

The brain activities are first decoded using the TSK decoder. The example of TSK with five rules is illustrated in Figure 8. It illustrates the corresponding membership functions of each fuzzy set obtained from five fuzzy rules in Figure 8A. Each membership function corresponds to a fuzzy linguistic description which is very high, high, medium, low, and very low. To provide further explanation, the antecedent parameters that include centers and standard variance are calculated, which are presented in Figure 8A. By the permutation of five linguistic descriptions, there can be twenty-five rules. Take five of them as an example. Based on the linguistic expressions and the corresponding linear function, the fuzzy rules can be shown as follows:

Rule 1: IF x₁ is very high AND x₂ is very high, (26)THENf1=[1.8268-0.0113x1-0.0042x2,0.2707-0.0144x1-0.0119x2,-0.2911-0.0038x1-0.0052x2,0.0370+7.2537e-4x1+9.4636e-4x2,-0.0289-8.1321e-4x1-0.0010x2]

Rule 2: IF x₁ is very low AND x₂ is medium, (27)THENf2=[-0.0954-4.2010e-4x1-0.0013x2,-0.0146-6.0436e-4x1-8.2373e-5x2,0.1950+0.0020x1-0.0026x2,1.5912-0.0135x1+5.2341e-4x2,-0.8781-0.0111x1-0.0034x2]

Rule 3: IF x₁ is high AND x₂ is very low, (28)THENf3=[-0.5874-0.0093x1-0.0024x2,1.7819-0.0032x1-0.0021x2,-1.0105-0.0120x1+8.3280e-4x2,0.0375-8.7808e-4x1-2.2777e-5x2,0.0814-0.0023x1+0.0027x2]

Rule 4: IF x₁ is low AND x₂ is low, (29)THENf4=[-0.1187+0.0017x1+5.0977e-4x2,-1.1958+0.0088x1+0.0153x2,3.1966+0.0236x1-0.0070x2,0.0791+0.0020x1+0.0027x2,-0.6716+0.0071x1-0.0145x2]

Rule 5: IF x₁ is medium AND x₂ is high, (30)THENf5=[0.1576+0.0030x1-0.0014x2,-1.0899-0.0099x1+0.0139x2,-0.7449+0.0117x1-0.0041x2,0.0255+6.8180e-4x1-1.8482e-5x2,2.4971+0.0055x1+0.0163x2] Based on these five fuzzy rules above, the mechanism of the decoding process is depicted in Figure 8B. The low-dimensional characteristics of different tasks are taken as the inputs of the decoder, and the label vector is predicted by the TSK system's decision output. The sum of the decision outputs f0i(i=1,2,…,5) based on the five rules is indicated by f⁰ and the decision outputs are handled by setting the maximal element in f⁰ to 1 and others to 0. In this example, the decision output is number 1.

The brain activities are decoded with the low-dimensional features extracted from a single brain lobe, two brain lobes, and three brain lobes (Figure 9). The result shows that the frontal and the occipital lobes are the top two highest accuracies when the features are extracted from a single brain lobe. The combination of the frontal and the occipital lobes is the highest accuracy in the two lobes combination. Additionally, in the three lobes combination, the frontal and occipital and parietal lobes are the highest effective. We hence infer that the frontal and the occipital lobes with high accuracy may be the key areas of the brain's visual cognition process. It can be concluded that finding the dominant brain regions can improve decoding efficiency and accuracy, which is essential for BCI.

In this work, we propose a method to decode cortical activity under visual stimuli from a single-trial level using the neural manifold. Through power spectrum analysis, it is found that visual cognition induces significant changes in the alpha band (Hogendoorn and Burkitt, 2018). Different visual stimuli lead to different spiking of neurons, which represent different information encoded by the human brain. The information is recorded in EEG and is exposed as the bifurcated manifold by dimension reduction methods. By further analyzing the latent factors extracted from the alpha band and using it to decode the brain activities, it is found that five major clusters that are associated with the visual stimuli distinctively occupy separate subregions.

The dimensionality reduction techniques are compared in terms of decoding performance. As excepted, it is confirmed that the latent factors obtained by t-SNE well-decoded brain activities. Decoding visual information using latent factors estimated by the t-SNE results in a superior performance to another case using PCA-estimated latent factors regardless of the decoder types, that is, TSK, NB, LSTM, and LKF. The research on the neurodynamic of the cerebral cortex has consistently uncovered low-dimensional manifolds that capture a crucial portion of neural variability (Arieli et al., 1996). PCA is widely used in motor cortex dynamics research. Gallego et al. (2018) confirmed the existence of a neural manifold in the motor cortex of monkeys by PCA. However, for behaviors whose dynamics explore non-motor brain areas, estimating neural manifolds using linear techniques requires latent factors of higher dimensions than the true dimensions in the data. EEG is commonly considered to have significant chaotic characteristics, linear methods provided poor estimates of the neural manifold (Figure 3). When brain activity is reduced to the same dimensions with xx, we can obtain a better recognition performance on the latent factors decoded by t-SNE during different cognitive tasks (Figure 4). The improvement in decoding brain activities using t-SNE suggests that the brain activities lie on a low-dimensional, non-linear manifold in the high-dimensional space and t-SNE can truly model the hidden state space of the brain in the cognitive process.

This may be due to TSK's good non-linear approximation ability as the performance of decoders is significantly more improved using latent factors estimated by t-SNE than PCA. Additionally, TSK could transform the decision output into a value between 0 and 1 through fuzzy rules, which is closer to human thoughts than other dichotomies of either 0 or 1 (Yu et al., 2019a). These findings suggest that we can design a simple fuzzy learning decoder (TSK) with t-SNE while achieving well performance. Multiple supervised machine learning methods that include TSK, NB, LSTM, KF are employed to classify different cognitive tasks based on two-dimensionality reduction techniques. Using latent factors estimated by t-SNE to decode visual information has higher accuracy than those estimated by PCA. Among these decoders, TSK with a promising classification result has been proved to be useful in the decoding analysis of neuroimaging data in many applications (Yu et al., 2019b).

According to the classification results, it is obtained that the number of brain lobes is not proportional to the classification performance (Figure 7). The frontal and the occipital lobes are found to be the most effective areas to classify different cognitive tasks. Compared with traditional methods, the visual decoding approach proposed has lower implementation difficulty and higher performance. Future work will consider finding key channels and frequencies to further reduce the amount of data analysis and improve decoding accuracy. The dominant brain lobes are found based on the decoding accuracy, but how to select features more efficiently and accurately is necessary to be further considered. In future works, it will be combined common spatial patterns (CSPs) with deep learning theory, further ensuring optimal channels and latent factors that are closely related to cognitive processes to achieve higher decoding accuracy.