Data mining based approach to study the effect of consumption of caffeinated coffee on the generation of the steady-state visual evoked potential signals

Journal Pre-proof Data mining based approach to study the effect of consumption of caffeinated coffee on the generation of the steady-state visual evoked potential signals Kishore K. Tarafdar, Bikash K. Pradhan, Suraj K. Nayak, Anwesha Khasnobish, Sumit Chakravarty, Sirsendu S. Ray, Kunal Pal PII:

S0010-4825(19)30385-3

DOI:

https://doi.org/10.1016/j.compbiomed.2019.103526

Reference:

CBM 103526

To appear in:

Computers in Biology and Medicine

Received Date: 2 July 2019 Revised Date:

28 October 2019

Accepted Date: 28 October 2019

Please cite this article as: K.K. Tarafdar, B.K. Pradhan, S.K. Nayak, A. Khasnobish, S. Chakravarty, S.S. Ray, K. Pal, Data mining based approach to study the effect of consumption of caffeinated coffee on the generation of the steady-state visual evoked potential signals, Computers in Biology and Medicine (2019), doi: https://doi.org/10.1016/j.compbiomed.2019.103526. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Data mining based approach to study the effect of consumption of caffeinated coffee on the generation of the steady-state visual evoked potential signals Kishore K Tarafdar1, Bikash K Pradhan1, Suraj K Nayak1, Anwesha Khasnobish2, Sumit Chakravarty3, Sirsendu S. Ray1, and Kunal Pal1* 1

Department of Biotechnology and Medical Engineering, National Institute of Technology,

Rourkela, India-769008 2

TCS Research and Innovation, Kolkata, India- 700156

3

Department of Electrical Engineering, Kennesaw State University, Marietta, GA, USA- 30060

*

Correspondence to Dr. Kunal Pal (E-mail: [email protected]; Phone: +91-876-336-6085)

Abstract The steady-state visual evoked potentials (SSVEP), are elicited at the parieto-occipital region of the cortex when a light source (3.5-75 Hz) flickering at a constant frequency, stimulates the retinal cells. In the last few decades, researchers have reported that caffeine enhances the vigilance and the executive control of visual attention. However, no study has investigated the effect of caffeinated coffee on the SSVEP response, which is used for controlling the braincomputer interface (BCI) devices for rehabilitative applications. The current work proposes a data mining-based approach to gain insight into the alterations in the SSVEP signals after the consumption of caffeinated coffee. Recurrence quantification analysis (RQA) of the electroencephalogram (EEG) signals was employed for this purpose. The EEG signals were acquired at seven frequencies of photic stimuli. The stimuli frequencies were chosen such that they were distributed throughout the EEG frequency bands. The prominent SSVEP signals were identified using the Canonical Correlation Analysis (CCA) method. Several statistical features were extracted from the recurrence plot of the SSVEP signals. Statistical analyses using the t-test and decision tree-based methods helped to select the most relevant features, which were then

classified using Automated Neural Network (ANN). The relevant features could be classified with a maximum accuracy of 97%. This supports our hypothesis that the consumption of caffeinated coffee can alter the SSVEP response. In conclusion, utmost care should be taken in selecting the features for designing BCI devices. Keywords:

SSVEP,

EEG,

Caffeine,

Canonical

Correlation

Analysis,

Recurrence

Quantification Analysis, Multilayer Perceptron Network

Introduction Steady-state visually evoked potentials (SSVEPs) are the brain signals that are generated in response to the stimulation of the retina. The stimuli are generated by flashing a light source at a constant frequency. The frequency of the SSVEP signals has been reported to be more prominent in the frequency range of 3.5 Hz to 75 Hz [1]. From a practical standpoint, only the lower frequency stimulus (<30Hz) can effectively evoke a strong SSVEP response [2]. The SSVEP signals are observed as the electrical activities at the localized brain sites. Such electrical events can be repeatedly produced if the stimuli are provided under controlled conditions. However, the actual mechanism of SSVEP signal generation is yet to be deciphered. Staring at a flickering light that flashes at a constant frequency stimulates the human visual pathway. The flickering frequency is radiated throughout the brain. This stimulation produces electrical signals in the brain at the base frequency of the flashing light, as well as at its harmonics. For example, when staring at a flickering stimulus of 5 Hz, the brain will produce SSVEP signals at the base frequency of 5 Hz and its harmonics at “5n” Hz, where n ∈ N (a natural number) [3]. Practically, there is a marked reduction in the power of the SSVEP signals from the second harmonics

onwards. This has been attributed to the low signal-to-noise ratio of the SSVEP signals at high frequencies and can be accounted for the brain dynamics that act as a low pass filter [4-5]. It is important to note that coffee is a ubiquitous beverage that is extensively consumed across the globe. Coffee contains a significant amount of caffeine, which is a potent psychoactive compound [6]. It has been reported that caffeine increases cognitive performance, alertness, ability to concentrate and focus, short-term memory, and decreased fatigue [7]. The action of caffeine is mainly based on the antagonism of the adenosine receptor through the blockage of cyclic adenosine monophosphate (cAMP) formation [8]. Caffeine has also been reported to increase the formation of catecholamine, which in turn, plays an essential role in imparting stimulant properties to caffeine. The effectiveness of the stimulation activity of caffeine is dependent on the alteration in the local rate of cerebral energy metabolism. Researchers have studied the effect of caffeine on brain activity by analyzing EEG signals [9-10]. Literature suggests that caffeine is having a potential impact on the frontal, parietal, occipital, and central region of the scalp [8]. It has also been reported that SSVEP, due to its robustness and high signal to noise ratio, is most widely used in brain-computer interface (BCI) devices for rehabilitative applications [11]. Hence, it is quite feasible that the consumption of coffee may significantly affect the performance of the BCI devices by altering the SSVEP response. Taking inspiration from the above, in this study, we have attempted to understand the impact of caffeinated coffee on the lower frequency SSVEP signals (<30Hz). For this purpose, EEG signals were recorded in the presence of photic stimuli of seven frequencies, both before and after drinking the caffeinated coffee. Canonical correlation analysis (CCA) was applied to identify the SSVEP signals. Recurrence quantification analysis was carried out to extract eleven numbers of the statistical features of the SSVEP signals. Statistical analyses using t-test and

decision tree based methods were employed to find the important features, which were then further utilized as the input to the Automated Neural Network (ANN). Materials and methods A 32-channel portable EEG machine (MAXIMUS 32, Recorders and Medicare Solutions Pvt. Ltd., India), capable of recording 22-channel EEG signals, was used to acquire the EEG signals from the volunteers. The processing and the analysis of the EEG signals were carried out using MATLAB (R2017a, MathWorks Inc., USA) software installed in a 64-bit Linux computer (Dell Inc., USA) with Gentoo operating system (Gentoo Foundation Inc., USA) that had a clock speed of 2.26 GHz and 12 Gigabytes of RAM. Figure 1 depicts the overall process of the proposed method. The detailed description of each block is given below. Data acquisition The EEG signals were recorded from six healthy male volunteers (age group: 22-28 years) after obtaining their written consent to participate in the study. All the participants were students of the National Institute of Technology, Rourkela. The permission for the acquisition of the EEG signals was taken from the Institute Ethical Committee (IEC), National Institute of Technology, Rourkela. The recording was carried out using a wet electrode-based clinical EEG machine. An electrode-cap that deployed the electrodes as per the international 10-20 electrode positioning system was used for the acquisition of the signals [12]. Figure S1 shows the montage used for the experiment with A2 as the reference electrode. The recording was carried out at a sampling frequency (fs) of 256 Hz. The volunteers were advised for not consuming any food or beverage before 2 h of the experiment. Further, they were strictly forbidden from any caffeinated food items within 12 h before the recording of the EEG signals. The volunteers were made to sit comfortably on a chair in a dim-lit room. After the electrode-cap was mounted on the head

(Figure S2), the electrode impedance was adjusted to less than 20 kΩ (Figure S3) by applying a commercially available EEG conductive gel. A set of white LED [13] (Figure S4) was positioned at a distance of 60 cm in front of the eyes of the volunteer. The participants were asked to relax for 2 min with their eyes closed. After the relaxation period, the volunteers were asked to stare at the light source. The LEDs were programmed to provide photic stimuli for 20 s at the frequencies of 3 Hz, 5 Hz, 10 Hz, 15 Hz, 20 Hz, 25 Hz, and 30 Hz, respectively. The photic stimuli frequency values were chosen in such a way so as to understand the SSVEP responses in all the EEG bands, namely, delta (0-4 Hz), theta (4-8 Hz), alpha (8-13 Hz), beta (13-30 Hz) and lower gamma (30-100 Hz) [14]. The stimuli frequencies of 3 Hz, 5 Hz, 10 Hz, and 30 Hz were generated to capture the delta, theta, alpha, and lower gamma waves. The beta waves are split into three sections, low beta waves (12.5–16 Hz, "Beta 1 power"), beta waves (16.5–20 Hz, "Beta 2 power"), and high beta waves (20.5–28 Hz, "Beta 3 power") [15]. Thus, three stimuli frequencies (15 Hz, 20 Hz, and 25 Hz) were chosen so as to activate the specific sub-beta bands. The volunteers were advised to abstain from blinking during this 20 s. A time gap of 10 s was maintained between the two consecutive photic stimuli frequencies. During this period, the volunteers were allowed to blink. The EEG data were recorded continuously from the starting of the relaxation period until 1 min after the cessation of the 30 Hz stimulus. Then, the participants were served with a cup of caffeinated coffee. Each serving of the coffee contained 120 ml of hot water, 1.5 g coffee (Nescafé classic), 6.5 g powdered milk, and 7 g sugar. After 10 min, the EEG signal was re-recorded using the same protocol [16]. The recorded EEG signals were labeled into two categories, namely, 'N' and 'M'. The signals recorded before drinking the coffee were categorized as ‘N’ (non-caffeinated), while the signals recorded after drinking the coffee were regarded as “M” (caffeinated). After the

recording of the signals was complete, the signal portion during the relaxation period, intermittent 10 s periods between the photic stimuli sessions, and post- 30 Hz photic stimulation were discarded. Each extracted signal contained 5120 numbers of samples (signal length of 20 s). Structuring of the extracted EEG signal data The extracted EEG signals were arranged in a four-dimensional (4D) array. A separate 4D array for each category (non-caffeinated and caffeinated) was made. The dimension of the array was in the form of

×

×

× , where, C, N, V, and S are the number of channels (C=22), sample

points in each channel (N=5120), volunteers participated (V=6), and different frequencies of the photic stimuli (S=7), respectively. Hence, the number of EEG signals in each 4D array was equal to

×

× , i.e., 22 x 6 x 7 (=924). This means that the total number of the extracted

EEG signals used in this study was 2 ×

×

×

= 1848 . This included both non-

caffeinated and caffeinated signals. The 4D array was used for further processing. Figure 2 shows the schematic representation of the compiled dataset for both categories. Signal preprocessing All the extracted EEG signals were filtered using a 6th order band-pass Butterworth filter ( 1.5-80 Hz). The maximum range of stimulus frequency was 30 Hz, which had its first harmonic at 60 Hz. Hence, to avoid any distortion of the first harmonic during filtering, a safe band of extra 20Hz was used, and the higher cutoff frequency was set at 80 Hz. After filtering, the signals were passed through a

rectangular window-based moving average filter (length of the window, w= 40). The optimum window length for the signals was computed using the method of the “sum of absolute differences” (SAD) [17]. Then, the power spectral densities (PSD) of the individual extracted EEG channels were computed to check the presence of the photic stimuli frequencies. The rectangular window of size of 640 samples (i.e., 1/8 of N, the total length of the signal) was used.

During computation, an overlapping window width of 320 samples during the consecutive calculations was maintained [18]. Detection of SSVEP signals and frequencies The filtered signals were then processed to identify the SSVEP signals from the 22 channels. This was achieved by canonical correlation analysis (CCA) [3, 19]. The computation of canonical correlation requires two basis vectors. These vectors are regarded as canonical variates, namely, x and y. The x variate is the filtered signals (x[n]). The canonical variate, x[n], was formulated by selecting the filtered signals from the three channels that lay within the close proximity and formed a small triangular area over the scalp and the underlying cortex (Figure S5, Table S1). Hence, the x[n] canonical variate contains three 1D signals that are obtained from the three channels. These 1D signals are arranged in a row-wise manner within the variate. The order of the signals in the canonical variate does not affect the outcome of the study as we are interested only in the correlation coefficient. In our research, there were 22 such local patches. These local patches were able to cover the entire cortex (Figure S5). On the other hand, the y[n] variate is the particular photic stimulus signal and its harmonics (y[n]). For generating the canonical variate, y[n], a square signal having an amplitude of unity and a frequency that is equal to the flickering stimuli was created. The square waves were then decomposed, and the base frequency (y[n1]) and the first harmonic (y[n2]) signals were extracted. Thereafter, the canonical variate was created by storing the y[n1] and the y[n2] vectors in a column vector of size 2x1 (Equation 1). =

(1)

where y(n1) is the base frequency, and y(n2) is the first harmonic signal.

The higher-order harmonics were not considered due to the filtering effect exerted by our brain on the higher-order harmonics [3]. Since there were 7 photic stimuli frequencies, 7 such canonical variates were created. The calculation of the canonical correlation coefficients (CCC) at the localized cortex areas was carried out using a specific x[n] variate and all the y[n] variates. During the processing, the CCC values were calculated using the 1 s filtered signals from each of the channels. After that, the highest CCC values that are obtained from all the computations were compared. If the CCC value was found to be highest using ya[n] canonical variate for the “a” stimuli frequency, then the filtered signals were considered to be SSVEP signals of “a” frequency. This process was repeated sequentially to analyze the 20 s filtered signal. During this period, an accumulator kept track of the number of the recognized SSVEP signals. Herein, since an independent CCC computation was carried out for each photic stimuli frequency, the confusion between the principal component of one frequency with the harmonic of the other frequency could be avoided. The results were then combined to generate a 3D data matrix, whose x-axis is the EEG channels, the y-axis is the volunteers, and the z-axis is the stimuli frequency (Figure 3). The said matrix had a dimension of 22 x 6 x 7. Each of the boxes (coordinates in the 3D matrix) contained counts of the recognized SSVEP signals, which is specific to a particular channel, stimuli frequency, and volunteer. Thereafter, the accumulated accepted SSVEP signal count values for each of the volunteers at each specific stimuli frequency were added for all the volunteers. This resulted in the formation of a 2D matrix whose x-axis is the EEG channels, and the y-axis (z-axis in the 3D matrix) is the stimuli frequency (Figure 4). After that, all the counts for specific channels were added. This resulted in the loss of resolution in the stimuli frequency axis (y-axis), which generated a 1D plot where the independent variable was the EEG channel, and the dependent variable was the summation of counts across each of

the EEG channels. Then, the count values were normalized using the maximum count values (Equation 2). =

× 100

(2)

where N is normalized value, f is SSVEP count of an EEG channel, and f max is Maximum SSVEP count of an EEG channel.

Eight channels having the highest normalized values were selected as the critical channels. This helped in recognizing the cortex regions with high brain activity during the stimulation process in the non-caffeinated category [3]. The same channels were selected for the analysis of the caffeinated category. Further, another study was conducted to understand the effect of photic stimuli frequency on the activation of the cortex region. The method was similar to the channelselection process, with the exception that instead of adding the count values for a single channel, the count values of the selected EEG channels were summed for a particular stimuli frequency. Then, the count values were normalized. The SSVEP activation at each of the stimuli frequencies was compared in both the categories. Recurrence Quantification Analysis The EEG signals from the selected channels were further studied for their non-linear recurrence behavior. Each EEG signal (20 s) was broken down into 4 segments of 5 s. The extracted 5 s EEG segments were then used for further processing. Firstly, a d-dimensional phase space plot was reconstructed from each of the chosen EEG signals. The phase space plot was reconstructed by computing the embedding dimension, and time delay from the signals. The embedding dimension “d” was calculated using the False Nearest Neighbor (FNN) algorithm (Equation 3) [20].

 Rd2+1 ( t,τ ) − Rd2 ( t,τ )    Rd2 ( t,τ )  

1/2

=

x ( t +τ ) − x ( tr +τ ) Rd ( t,τ )

> R tol

(3)

where Rd is the measure of distance in phase space with embedding dimension d, and R

is the

tolerance threshold. The time delay (τ) was obtained by computing the average mutual information (Equation 4) [21]. ) ! ", $ = ∑.+, ∑*+, &'( /. 0. 123

456 78 98

(4)

45 78 46 :9; <

where &'( /. 0. is the joint probability distribution of the discrete-time variables X and Y at the instant x = /. , y = 0. .

&' /. is the probability distribution of X at the instant, x = /. . &( 0. is the probability distribution of Y at the instant, y = 0. . Subsequently, the phase space matrix was computed (Equation 5).

/, / H B /C G B / H−J A/ F A A D F = A / H − 2J . A . F A A . F A . @/) E @/ H − +1 J

/ H−J / H − 2J / H − 3J . . .

/ H − 2J / H − 3J / H − 4J . . .

. . . . . .

. . . . . .

/ H− K−1 J G . F . F . F . F / H− K− J E

(5)

where d is the embedded dimension, and J is the time delay. ⃗

/. ∈ O P ∀ R = 1,2, . . . . . ,

(6)

⃗

where /. represents the trajectory. The right-hand side matrix of the equation represents the phase space matrix, and the left-hand side vector indicates the trajectory of the phase space [22]. Subsequently, the recurrence plot and the recurrence quantification analysis (RQA) were carried out. The trajectory of this ddimensional phase space forms the line of identity (LOI) of the recurrence plot when the same is projected over a 2D map. The Euclidean distance between the R ST state of the trajectory, and

remaining all other states in the phase space was computed. The computed distances were mapped into the color information using the jet colormap function of the MATLAB software. Thus, a colored 2D recurrence plot was obtained for visual analysis of the recurrence behavior of the SSVEP signals [23]. The classical recurrence plot was also derived from the phase space plot for recurrence quantification analysis. It is assumed that a sphere of radius ‘r’ has its center on the trajectory of the reconstructed phase space. The phase space states (or points) that lie inside the sphere are labeled as “1”. On the other hand, the points that lie outside the sphere are marked as “0”. Traversing the sphere through the entire phase space trajectory resulted in the formation of a binary-valued 2D vector. This resulted in the generation of the recurrence plot of the signal [24]. This recurrence plot had multiple patterns of vertical, horizontal, and diagonal lines. The lines are formed due to the occurrence of the phase space states with a label of “1”. The line patterns provide qualitative information on the repetitive behavior of the SSVEP signals. The repetitive expression was quantified by the recurrence quantification analysis (RQA) [25]. The implementation of RQA resulted in the generation of 11 nonlinear statistical features from the recurrence plot. These features were: Recurrence rate (RR), Determinism (DET), Maximal line length (LMAX), Entropy (ENT), Laminarity (LAM), Trapping time (TT), Maximum length of the vertical lines (VMAX), Ratio of determinism and recurrence rate (RATIO), Average neighbours of each trajectory point (AN), Divergence (DIV), and Average length of all diagonals (AD) [23, 25-26]. Hence, a total of 56 (=8 channels x 7 stimuli) feature matrices were generated. Each of the feature matrices contained 11 features, and two labels “N” and “M”.

Feature reduction and signal classification The RQA features were grouped into a feature table with two class labels, i.e., “N” (noncaffeinated) and “M” (caffeinated). The critical features, which could better represent the non-

caffeinated and caffeinated SSVEP signals were identified using a t-test, and 3 non-linear decision tree-based methods, namely, Classification and regression trees (CART), Boosted trees (BT), and Random forest (RF) . The t-test is a linear method, which tries to identify the statistical significance of the mean values of the features [27]. It is used to reject the hypothesis that the group means are identical. If any difference is found between the group means, then it reflects the actual differences in the sample population of both the groups [28]. A p-value is used to calculate the statistical significance of a test within a specified interval by assuming that the group means are identical. p-value < 0.05 refers that the group means are different, and a value >0.05 indicates identical groups with equal mean [29]. CART is based on the methodology of binary recursive partitioning [30]. It classifies the given dataset using a decision tree. Since it does not provide the best results in classification, it is widely used as a preprocessing step of feature selection [31]. Further, the CART algorithm automatically produces variable importance (predictor importance) for each attribute. This ranking is based on the contribution of a feature in making the decision, and thus it identifies the critical features that were essential during the classification process [30]. The boosted tree method improves the performance of a single model by fitting a number of models and then merging the same for prediction [32]. In the initial stage, equal weights are assigned while training the decision tree. In the next step, the weights of the observation are optimized. Though each simple tree has a weak predictor accuracy, the accuracy can be improved when the trees are combined [33]. While selecting the splitting feature in the algorithm, a sorted list of relatively important features is automatically generated. The predictor importance of the boosted tree is the average of all the trees [34]. Random forest (RF) is a type of ensemble decision tree [35]. The construction of an RF is mainly based on two steps. First, dividing the training dataset into many subsets and using each subset to design a decision tree.

The splitting during the tree growth depends upon the number of features selected from a group of candidates that are randomly chosen from all the features. The importance of a feature can be calculated by permuting the values of one feature for all the samples and then optimizing the OOB (out of the bag) accuracy [36]. In each case, a predictor importance value signifies the more important feature for predicting the output. In the non-linear methods, a predictor importance value >0.95 was considered as statistically significant. The statistically significant features, which were identified using linear and non-linear statistical methods, were chosen as the input data for training the artificial neural networks. Various configurations of Multilayer Perceptron Neural Networks (MLPNN), a class of supervised feed-forward neural networks were trained by adjusting different training parameters of the neural networks [37]. The MLPNNs are widely used with various biomedical signals for solving pattern recognition and classification problems [38-39]. All the MLPNNs that were tested with our data had three-layer architecture, which consisted of an input layer, a hidden layer, and an output layer. The training algorithm, the activation function of the hidden layer, the activation function of the output layer, error function, and the number of neurons in the hidden layer were varied during the training of the MLPNNs. Levenberg-Marquardt [40], and BFGS Quasi-Newton Backpropagation training algorithms [4142] were used in our study. The hidden layer activation functions, used in this study, include hyperbolic tangent sigmoid transfer, log-sigmoid transfer, and exponential transfer functions [41, 43-44]. All these functions were non-linear functions. The output layer activation functions that were employed include Hyperbolic tangent sigmoid transfer (tansig), Log-sigmoid transfer (logsig), exponential transfer, Linear transfer (purelin), and Softmax transfer functions. Among them, the purelin and the Softmax functions were linear transfer functions [45]. The performance functions (sum of squared error (sse), and cross-entropy (ce)) were also varied during the

training process [46]. The number of Perceptrons in the hidden layer was varied from 3 to 25. 2000 numbers of MLPNNs were trained to find the optimum network for each of the 56 feature tables.

Results Signal preprocessing The progressive mapping of the brain in response to the seven photic stimuli was obtained from the software of the EEG machine (Figure 5). The brain mapping helps to identify the activity of the brain when a stimulus is given. In both the non-caffeinated, and the caffeinated categories, it was observed that during the initiation of a particular photic stimulus, the frontal cortex was activated. However, as time progressed, the activation was mainly found in the parieto-occipital region. The stimulation of the cortex in the parieto-occipital region was quite expected as this area processes the related visual signals [47-48]. Interestingly, the analysis of the brain mapping in the caffeinated category indicated that the frontal cortex was also activated during the stimulation process when the stimuli frequency was 5 Hz, 10 Hz, and 20 Hz. Further, at 15 and 25 Hz stimuli frequency, the activation of the frontal cortex was similar in both the noncaffeinated and caffeinated categories. The activation of the frontal cortex was not as strong as that of the cortex in the parieto-occipital region. The acquired EEG signals were applied to a bandpass filter to remove unwanted frequency components (Figure 6). Further, the bandlimited signals were subjected to moving average filtering. The taps (window size) of the moving average filter were optimized using the SAD algorithm. The SAD algorithm varies the number of taps of the moving average filter. It then correspondingly calculates the sum of the absolute differences between the original signal and

the processed signal for a particular number of the taps. Subsequently, the SAD is plotted against the number of taps. The plot initially increases in a linear fashion, which then attains a plateau phase. The number of taps where the plot reaches the plateau value is regarded as the optimized number of the taps. The SAD plot of the representative non-caffeinated and caffeinated categories has been shown in Figure S7. From the plots, it was observed that for all the signals, the SAD plots were flattened at the tap value of 40. Hence, the optimal width of the moving average filter for the pre-processing of the EEG signals was selected as 40. The use of a moving average filter helps to suppress the uncorrelated signal components with the consequent improvement in the signal-to-noise ratio (SNR) of the EEG signals [49-50]. In our case, we have also found a similar improvement in the EEG signal quality. It can be clearly seen that most of the perturbations in the signals were suppressed. The spectral analysis (PSD) of the filtered EEG channels revealed that the 3 Hz, 5 Hz, 10 Hz, 15 Hz, 20 Hz, and 30 Hz frequencies were significantly present in the PG1, PG2, FP1, FP2, FZ, CZ, PZ, P3, P4, O1, O2, OZ, T5, and T6. The PSD plot of Channel O2 at different photic stimuli for both the non-caffeinated and caffeinated states has been shown in Figure S6. The spectral analysis of the remaining EEG channels suggested that the response of the remaining channels were minimal to the visual stimuli frequencies.

SSVEP detection The normalized plot of the CCC count per channel has been plotted in Figure 7. From the normalized plot, eight numbers of the channels with the highest normalized counts were regarded as the important channels (represented in green color). These channels could be useful in understanding the effect of coffee on the cortex [51-52]. It was observed that the O1, O2, OZ,

PZ, P4, T4, T5, and T6 channels were mainly activated during the non-caffeinated state (Figure 7a). The results suggested that most of the occipital cortex region and the partial parts of the parietal and the temporal cortex regions were activated. This was quite expected because the photic stimuli are processed in the parietooccipital region of the brain [53-54]. However, previous literature suggests that the temporal region of the cortex is also activated during the photic stimulation [6, 55-56]. During the caffeinated state, the O1, O2, OZ, PZ, T4, T5, and T6 were mainly activated (Figure 7c). Interestingly, it was observed that all the activated regions during the non-caffeinated state were also activated during the caffeinated state. The only exception being the P4 channel, which was suppressed. After that, the percentage increase/ decrease in the SSVEP activity during the caffeinated state was calculated (Equation 7).

U=

VWX X

× 100

(7)

where C is the percent increase in SSVEP occurrence, A is the number of SSVEP occurrences after the stimulus, and B is the number of SSVEP occurrences before the stimulus. The SSVEP activity during the non-caffeinated state was taken as the control value. A plot of the analysis is provided in Figure 7c. From the investigation, it was found that there was a significant increase in the SSVEP activity of the frontal and the frontoparietal region during the caffeinated state. This result is in agreement with the progressive mapping that was obtained from the EEG machine. Subsequently, the effect of stimuli frequency on the cortex activity was also predicted. The eight important channels that were selected (O1, O2, OZ, PZ, P4, T4, T5, and T6) in the previous study were also used in this study. It was observed that the 10 Hz and 15 Hz stimuli frequencies were most active in the generation of the SSVEP signals in the non-caffeinated state (Figure 8a). This was followed by the activation of the 3 Hz stimulus. The activation of the SSVEP by the 5 Hz, and the 20 Hz stimuli frequencies was similar and was slightly lower than

the 3 Hz stimulus. The activation of the SSVEP by the 25 Hz stimulus frequency was lower than the 5 Hz, and the 20 Hz stimuli frequencies. The activation of the SSVEP by the 30 Hz stimulus frequency was lowest. A similar trend was followed during the caffeinated state also (Figure 8b). However, the activation of the 5 Hz stimulus frequency was slightly higher than the 20 Hz stimulus frequency. Previous literature suggests that there is a decreased SSVEP activity at higher frequencies [57-58]. In our study, we have also found that the SSVEP activity was comparatively lower in the higher stimuli frequencies (20 Hz, 25 Hz, and 30 Hz). The percentage increase/ decrease in the SSEVP activity as a function of the stimuli frequency was calculated (Figure 8c). It was observed that there was an increase in the SSVEP activity during the caffeinated state when the stimuli frequency was 3 Hz, 5 Hz, and 30 Hz. The increase in the SSVEP activity due to 30 Hz was found to be exceptionally high. There was no change in the SSVEP activity during the 10 Hz stimulus. The SSVEP activity during the 15 Hz, 20 Hz, and 25 Hz was suppressed.

Recurrence Quantification Analysis The fraction of the false nearest neighbors (i.e., the ratio of the number of false neighbors to the total pairs of points) is plotted as a function of the embedding dimension. The embedding dimension, for which the fraction of the false nearest neighbors goes below the threshold value, represents the optimal embedding dimension for the phase space plot. Figure 9a depicts a typical plot of the fraction of the false nearest neighbors vs. embedding dimension for the EEG signal of the O2 channel at 10 Hz photic stimulus for the caffeinated state. Then, the average mutual information (AMI) function of the given signal with its time-delayed version signal was calculated [26]. The optimal time delay is obtained from the first minimum of the AMI function

vs. time delay plot (Figure 9b). The phase space plot was generated (Figure 9c) using the embedding dimension and time delay information. This was further used for the generation of the recurrence plots of the extracted 5 s portions of the SSVEP signals (Figure 9d-e). The RQA of the recurrence plot of each of the 5 s extracted signals provided 11 statistical features. For each channel, an RQA feature matrix that contained 24 numbers of features for each category (non-caffeinated and caffeinated) was obtained. Hence, for the photic stimuli of 7 different frequencies, and 8 activated EEG channels, 56 numbers of the RQA feature matrices were obtained.

Feature reduction and signal classification Each of the 56 numbers of the RQA feature matrices was subjected to statistical tests, namely, ttest, CART, BT, and RF, for identifying the statistically significant RQA features. The statistically significant RQA features for 3Hz visual stimuli have been tabulated in Table 1. The detailed information about the important features has been provided in Table S2-Table S7. The representative tree architectures that were used for predicting the important features have been shown in Figure S8. Except for the SSVEP signal from the O1 channel, the t-test could not find any statistically significant RQA feature from any other EEG signals. However, the decision tree-based methods, i.e., CART, BT, and RF, were able to detect statistically significant RQA features from all the SSVEP signals. This suggested the presence of a nonlinear function that can segregate the RQA features from the two classes. After that, the statistically significant RQA features were used in various combinations as input for classification of the non-caffeinated and the caffeinated states using MLPNNs. The classification accuracies of the MLPNNs for all the SSVEP dominant EEG channels at all the stimuli frequencies have been summarized in Table 2. A classification efficiency of ≥95% was

considered to be significant in our study. It was observed that photic stimuli, which generated high SSVEP responses on the cortex (i.e., 10 Hz (alpha) and 15 Hz (low beta)), exhibited low classification accuracy with a maximum of 87 % (at EEG channel O2 and T4) for 10 Hz photic stimulus and 97 % (at EEG channel P4) for 15 Hz photic stimulus. The other 7 SSVEP predominant EEG channels showed relatively lower accuracy in classifying the caffeinated SSVEPs from non-caffeinated SSVEPs at 10 Hz and 15 Hz flickering visual stimuli. The average classification efficiencies at these stimuli frequencies were lowest. Conversely, the photic stimuli of frequencies >15 Hz (beta, high beta, low gamma) showed lower classification efficiency. The 3 Hz and 5 Hz photic stimuli generated high SSVEP responses, which were lower than the 10 Hz and 15 Hz photic stimuli. In these channels, the maximum classification accuracy was found to be 97% at the P4 channel and 95% at the O2 channel. Table S8 lists the detailed architecture of all the 56 MLPNNs.

Discussion This study was aimed at investigating the effect of caffeinated coffee on the activation of the SSVEP signals. The EEG signals were acquired from six volunteers during pre- and postconsumption of coffee at 7 different frequencies (3-30 Hz) of photic stimuli. The acquired EEG signals were analyzed for the identification of the SSVEP components. A progressive map of the brain for the stimulated regions showed that there was an activation of the frontal lobe at the onset of the stimuli. The activation of the frontal lobe may be attributed to the sudden stimulation of the rods and the cones at the onset of the stimuli. Thereafter, there was gradual activation of the parieto-occipital region over a 20 s time period. The acquired signals were then preprocessed to improve the signal-to-noise ratio. The identification of the SSVEP cortex areas was made by

employing the CCA method. Further, the SSVEP frequency component was also identified using the CCA method. Lin et al. (2006) had first proposed the method of CCA for the identification of the SSVEP signals. The group used the CCA method for locating both the activated cortex areas and the SSVEP frequencies. It was observed in their study that CCA provided a higher SSVEP signal recognition accuracy as compared to the conventional power spectral density analysis method [3]. In a similar study, conducted by Hakvoort et al., (2010) on 7 volunteers with 7 different stimuli, similar robustness in the detection of SSVEP signals using the CCA method was reported. This may be attributed to the inclusion of the harmonic frequencies in the CCA computation [59]. In our study, eight EEG channels (O1, O2, OZ, PZ, T4, T5, and T6) were found to contain SSVEP signals. Thus, a local patch of 8 EEG channels was formed on the cortex, which was considered as the area on the scalp with maximum SSVEPs. Zheng Yan et al. (2011) had analyzed functional connectivity in different brain regions during the SSVEP response [60]. Using PSDA (power spectral density analysis), the authors suggested a maximal SSVEP response in the lower frequency stimuli (alpha wave) and a lower SSVEP response in the beta wave. Similar results have been observed in the current study. The SSVEP count was comparatively higher in the frequency range of 3-15 Hz (maximum at 10Hz) than that of the higher band frequency range (25-30Hz), as shown in Figure 8. The enhanced SSVEP response in the alpha wave may be explained due to the increased information flow from the visual cortex to the occipital region after coffee ingestion. These findings are referring to direct connectivity among the band of flickering frequency, caffeine, and their response in the occipital region. However, more detailed analysis is required for its validation. Further, the results of the flow gain method suggested that the SSVEP response is not confined to the occipital region. It is also found in the parietal region. We found similar results in this study, i.e.,

there is an increased SSVEP response in the parietal region. It is well known that the brain's ability to attend visual stimuli based on their spatial location requires the parietal cortex [61]. It has been reported previously that the parietal cortex predicts the retinal consequences of the eye movements, updates the retinal coordinates of the remembered stimuli, and generates a continuously accurate representation of the visual space [61]. It also acts as a critical locus for the mental illustration of the visual world [62]. The authors have hypothesized that the parietal cortex may serve as a hub for information exchange between the visual and occipital cortexes. Ding et al. (2006), in a detailed study of the dependence of SSVEP power on the visual stimulus frequency between 2.5 Hz and 20 Hz, established that the SSVEP signals were strongly dependent on the flicker frequency. The group reported that the highest power was observed at 9.2 Hz (alpha) and 10 Hz (alpha) visual stimulus [63]. It was also reported that at higher flickering frequencies, the SSVEP signals were weak. Similar observations were made in our study where we have found that the SSVEP signals were very weak during the flickering frequencies of 20 Hz (beta), 25 Hz (high beta) and 30 Hz (low gamma). The weakest SSVEPs are generally found in the high-frequency range. However, in our case, 30 Hz flicker frequency significantly increased the SSVEP signals in the caffeinated states. Since the brain is often considered as a nonlinear and chaotic system [64-65], RQA was chosen for feature extraction of the signals from the SSVEP signals. Due to the limited size of the dataset, each signal was divided into 4 parts, each of 5 s duration, to construct a feature space that contained 48 instances [66]. The statistical methods, namely, t-test, CART, BT, and RF, were chosen for feature reduction to find statistically significant features. This helped to reduce any bias that could have been introduced in the analysis if only one method was considered. The feature reduction method helps to reduce the chances of over-fitting while training the MLPNNs.

It can be observed from the performance matrix of the MLPNNs (Table 2) that the SSVEP signals obtained due to 10 Hz (alpha) and 15 Hz (low beta) stimuli, showed lowest classification accuracy in each of the activated EEG channels. This may be accounted for the generation of strong SSVEP signals under both conditions. Due to this reason, the differences between the non-caffeinated and the caffeinated categories could not be segregated in the feature space. However, the P4 channel was an exception, where the classification efficiency was found to be 97%. It is important to note over here that during the identification of the SSVEP signals, it was found that during the caffeinated state, there was a suppression of the SSVEP signals in the P4 channel. This might have resulted in better classification accuracy. Also, it was observed from the CCA that a lesser amount of SSVEPs was produced by 25 Hz (high beta) stimulus. This was also observed during the PSD analysis as no significant peaks were observed in most of the EEG signals at 25 Hz (high beta) stimulus. Interestingly, the overall classification accuracy at 25 Hz stimulus was found to be relatively higher as compared to the other stimuli. The overall classification efficiency in the T4, T6, and P4 channels was significantly higher than the other channels, which might have resulted in the improved classification. This observation may be explained by the fact that the caffeine improved the SSVEP signals in the EEG channels, where the SSVEP signals were weak. A similar observation was also made for the 20 Hz (beta) stimulus. At 5 Hz (theta) stimulus, both the SSVEP detection and the caffeinated data classification exhibited good results. A somewhat similar result was also obtained for a 3 Hz (delta) stimulus. Caffeine is often found to increase energy in the delta band [67]. The highest classification efficiency for the 30 Hz stimuli frequency was 95%, with an average efficiency of 81.5%. Table 3 compares the proposed method with some state of the art SSVEP detection methods. Chen et al. (2015) have proposed a modified CCA for the detection of SSVEP signals

[68]. A filter bank canonical correlation analysis had been used to incorporate the fundamental and harmonic frequency components to enhance classification accuracy. The sub-bands of the SSVEP signals were extracted using the filter bank. In the next step, the CCA coefficients between these SSVEP signals and the reference sinusoidal signal was calculated and consequently used for target identification. The average classification accuracy achieved following this method was found to be 91.95%. In another study, a multiset CCA method was used for SSVEP detection [69]. The detection efficacy of the SSVEP signal was improved by optimizing the reference signal that was based on the common SSVEP features. The authors implemented the correlation maximization between the EEG signal and the optimized reference signal with sine-cosine waves at each stimulus frequency. The detection accuracy was calculated at different harmonics. The detection accuracy was >90% for a time window of four. Wang et al. (2016) have suggested a spatiotemporal method of feature extraction with multivariate linear regression (MLR) that improved the detection accuracy of the SSVEP signal [70]. A label matrix was constructed from the training data using principal component analysis (PCA) and MLR to find the optimal discriminative subspace. The testing data after reducing the dimensionality were projected to the label matrix. In the next stage, K-NN was applied to classify the subspace feature. It was observed that despite the time window used in most of the cases, the detection accuracy was ~95%. Yangsong Zhang et al. (2018) has explained a two-stage correlation component analysis (CORRCA) method to avoid the assumption used in the CCA method that canonical vectors must be orthogonal [71]. Spatial filters for each frequency was generated in the initial stage, and then reference signals were generated by group averaging across multiple blocks. Further, the correlation coefficients were computed using these two signals. The highest accuracy achieved in this case was ~80%. In comparison to different methods shown in Table 3.

The proposed method shows a good classification accuracy of 97%. However, it must be noted that the study protocol used in the current work was different from the above mentioned literature. Hence, providing a fair comparison may not be possible. The computational time was calculated using MATLAB 2017a environment on a 64 bit Linux computer with a Gentoo operating system (Gentoo Foundation Inc., USA) that had a clock speed of 2.26 GHz and 12 Gigabytes of RAM. The total execution time for SSVEP detection in a 22 channel EEG data was found to be 16.67 s, which is higher than the standard CCA method used in other studies [68-69]. The difference in the computational time may be due to the changes in the study protocol. In the current study, a single flickering frequency stimulus was considered at a time to detect its region of activation in the brain. Hence, to cover the whole spatiotemporal regions of the brain, 22 channel EEG signals were used, and the EEG signal from each channel was compared with the reference signal of each frequency stimuli. This may have added to its high computational efficiency. Further, the temporary variables and reference signals that are used for the calculation of the CCA coefficients are usually previously computed and stored. Hence, the computation time does not include these steps for calculating computational time during online processing [68]. However, the current study employs offline analysis, and no such prior calculation has been adopted. This may explain the reason for the increased computational time in our case.

Limitations and future work The participation of the lower number of volunteers is a limitation of the current study. Although we have achieved good performance accuracy, the study has to be repeated with a higher number of volunteers for proper validation of a larger population. Inter-participant differences may exist among the volunteers due to their different habituation of coffee ingestion. However, the current

study focuses on average group-level analysis. It is worthy of exploring the individual differences in response to SSVEP stimulus after coffee ingestion using the proposed method in the future. In this study, a total of 22 channels was used for collecting the EEG data. The smaller electrode systems lead to a reduced spatial resolution. On the contrary, a larger electrode system reduces comfort and extended application time [72]. A number of studies have used the 22electrode EEG system for BCI applications using motor imaginary, and source localization [7375]. A smaller electrode system was used in our experiment considering the above disadvantages of the larger electrode system. But, in the future, we would like to continue this study with the larger electrode system to have an in-depth insight into the spatiotemporal response of the brain signals. As described in the previous literature, the effect of SSVEP is more pronounced in the frequency range below 30 Hz [2]. Hence, the maximum frequency range considered in this study was 30 Hz. However, different research can be conducted to find the effect of coffee on higher stimuli frequency. The effect of caffeinated coffee and other drugs on different stimuli such as ERP, P300, and hybrid systems (ERP-SSVEP, P300-SSVEP, etc.) can also be studied in the future to understand the difference in the brain response to different drugs. Different feature reduction and classification methods may also be explored to improve the efficiency further.

Conclusion The EEG signals were acquired from 6 volunteers at the photic stimuli of 3 Hz, 5 Hz, 10 Hz, 15 Hz, 20 Hz, 25 Hz, and 30 Hz. These photic stimuli were selected such that they appear in the delta, theta, alpha, beta, and low gamma bands of the EEG signals, respectively. The acquisition of the EEG signal was carried out both during the pre- and post-consumption of coffee. The progressive mapping of the brain was also obtained from the EEG recording machine. Visual observation of the brain maps suggested the initial activation of the frontal lobe with the

exposure of the photic stimuli, which gradually moved to the parietooccipital region as the time progressed. The preprocessing of the acquired EEG signals was performed using band-pass, and moving average filtering. The preprocessed EEG signals were subjected to CCA for the identification of the SSVEP signals. The results suggested that 8 channels in the parieto-occipital region of the cortex can provide information about the SSVEP signals. The effects of caffeine on the SSVEP signals were studied by the data mining-based approach. Since the brain is considered to be a chaotic nonlinear system, 11 nonlinear statistical features were calculated for the SSVEP signals using RQA. The t-test, along with 3 nonlinear tests (CART, BT, and RF), was performed on each of these feature tables to identify the statistically significant features that could well classify the non-caffeinated and the caffeinated states. The significant features from each feature tables were then used to train MLPNNs to differentiate between the non-caffeinated and the caffeinated states of the SSVEP signals it was observed that channel O2 showed a classification accuracy of ≥95% at 5 Hz (theta), 20 Hz (beta) and 30 Hz (low gamma) stimuli frequency. The channel P4 showed a classification accuracy of 97% at 3 Hz (delta) and 15 Hz (low beta) stimuli frequency. The O1 channel had a classification accuracy of 95% at 20 Hz (beta) visual stimulus. In gist, it can be concluded that the consumption of caffeinated coffee, in general, does not significantly affect the generation of the SSVEP signals in the parietooccipital channels. However, the SSVEP signals in the P4 channel, present in the parietal region, are markedly suppressed. Further, it was found that the consumption of coffee significantly increases the SSVEP signal content in the frontal lobe channels. Hence, the selection of the EEG channels and the SSVEP signal feature should be carefully carried out so that the performances of the BCI devices are not compromised.

REFERENCES: [1] F. Beverina, G. Palmas, S. Silvoni, F. Piccione, S. Giove, User adaptive BCIs: SSVEP and P300 based interfaces, PsychNology Journal, 1 (2003) 331-354. [2] P. Xu, C. Tian, Y. Zhang, W. Jing, Z. Wang, T. Liu, J. Hu, Y. Tian, Y. Xia, D. Yao, Cortical network properties revealed by SSVEP in anesthetized rats, Scientific reports, 3 (2013) 2496. [3] Z. Lin, C. Zhang, W. Wu, X. Gao, Frequency recognition based on canonical correlation analysis for SSVEP-based BCIs, IEEE transactions on biomedical engineering, 53 (2006) 2610-2614. [4] Q. Wei, M. Xiao, Z. Lu, A comparative study of canonical correlation analysis and power spectral density analysis for SSVEP detection, 2011 Third International Conference on Intelligent HumanMachine Systems and Cybernetics, IEEE, 2011, pp. 7-10. [5] D. Lesenfants, D. Habbal, Z. Lugo, M. Lebeau, P. Horki, E. Amico, C. Pokorny, F. Gomez, A. Soddu, G. Müller-Putz, An independent SSVEP-based brain–computer interface in locked-in syndrome, Journal of neural engineering, 11 (2014) 035002. [6] F. Burdan, Caffeine in Coffee, Coffee in health and disease prevention, Elsevier2015, pp. 201-207. [7] M.M. Lorist, M. Tops, Caffeine, fatigue, and cognition, Brain and cognition, 53 (2003) 82-94. [8] J.V. Higdon, B. Frei, Coffee and health: a review of recent human research, Critical reviews in food science and nutrition, 46 (2006) 101-123. [9] J. Snel, M.M. Lorist, Effects of caffeine on sleep and cognition, Progress in brain research, Elsevier2011, pp. 105-117. [10] M. Siepmann, W. Kirch, Effects of caffeine on topographic quantitative EEG, Neuropsychobiology, 45 (2002) 161-166. [11] Z. İşcan, V.V. Nikulin, Steady state visual evoked potential (SSVEP) based brain-computer interface (BCI) performance under different perturbations, PloS one, 13 (2018) e0191673. [12] G.H. Klem, H.O. Lüders, H. Jasper, C. Elger, The ten-twenty electrode system of the International Federation, Electroencephalogr Clin Neurophysiol, 52 (1999) 3-6. [13] C.S. Herrmann, Human EEG responses to 1–100 Hz flicker: resonance phenomena in visual cortex and their potential correlation to cognitive phenomena, Experimental brain research, 137 (2001) 346-353. [14] J.S. Kumar, P. Bhuvaneswari, Analysis of Electroencephalography (EEG) signals and its categorization–a study, Procedia engineering, 38 (2012) 2525-2536. [15] M. Rangaswamy, B. Porjesz, D.B. Chorlian, K. Wang, K.A. Jones, L.O. Bauer, J. Rohrbaugh, S.J. O’connor, S. Kuperman, T. Reich, Beta power in the EEG of alcoholics, Biological psychiatry, 52 (2002) 831-842. [16] M. Saifudinova, M. Bachmann, J. Lass, H. Hinrikus, Effect of Coffee on EEG Spectral Assymmetry, World Congress on Medical Physics and Biomedical Engineering, June 7-12, 2015, Toronto, Canada, Springer, 2015, pp. 1030-1033. [17] J. Fowers, G. Brown, P. Cooke, G. Stitt, A performance and energy comparison of FPGAs, GPUs, and multicores for sliding-window applications, Proceedings of the ACM/SIGDA international symposium on Field Programmable Gate Arrays, ACM, 2012, pp. 47-56. [18] N. Srinivasan, Cognitive neuroscience of creativity: EEG based approaches, Methods, 42 (2007) 109-116. [19] D.R. Hardoon, S. Szedmak, J. Shawe-Taylor, Canonical correlation analysis: An overview with application to learning methods, Neural computation, 16 (2004) 2639-2664. [20] H.D. Abarbanel, M.B. Kennel, Local false nearest neighbors and dynamical dimensions from observed chaotic data, Physical Review E, 47 (1993) 3057. [21] J. Martinerie, A.M. Albano, A. Mees, P. Rapp, Mutual information, strange attractors, and the optimal estimation of dimension, Physical Review A, 45 (1992) 7058. [22] N. Marwan, M.C. Romano, M. Thiel, J. Kurths, Recurrence plots for the analysis of complex systems, Physics reports, 438 (2007) 237-329. [23] Y. Chen, H. Yang, Multiscale recurrence analysis of long-term nonlinear and nonstationary time series, Chaos, Solitons & Fractals, 45 (2012) 978-987.

[24] H. Yang, Multiscale recurrence quantification analysis of spatial cardiac vectorcardiogram signals, IEEE transactions on biomedical engineering, 58 (2010) 339-347. [25] C.L. Webber Jr, N. Marwan, Recurrence quantification analysis, Theory and Best Practices, (2015). [26] S.K. Nayak, A. Bit, A. Dey, B. Mohapatra, K. Pal, A review on the nonlinear dynamical system analysis of electrocardiogram signal, Journal of healthcare engineering, 2018 (2018). [27] T.C. Urdan, Statistics in plain English, Routledge2011. [28] T.K.P. Shri, N. Sriraam, Comparison of t-test ranking with PCA and SEPCOR feature selection for wake and stage 1 sleep pattern recognition in multichannel electroencephalograms, Biomedical Signal Processing and Control, 31 (2017) 499-512. [29] I. Kamkar, S.K. Gupta, D. Phung, S. Venkatesh, Stable feature selection for clinical prediction: Exploiting ICD tree structure using Tree-Lasso, Journal of biomedical informatics, 53 (2015) 277-290. [30] S. Chebrolu, A. Abraham, J.P. Thomas, Feature deduction and ensemble design of intrusion detection systems, Computers & security, 24 (2005) 295-307. [31] L. Gómez-Chova, J. Calpe, E. Soria, G. Camps-Valls, J.D. Martin, J. Moreno, CART-based feature selection of hyperspectral images for crop cover classification, IEEE, pp. III-589. [32] A.M. Youssef, H.R. Pourghasemi, Z.S. Pourtaghi, M.M. Al-Katheeri, Erratum to: Landslide susceptibility mapping using random forest, boosted regression tree, classification and regression tree, and general linear models and comparison of their performance at Wadi Tayyah Basin, Asir Region, Saudi Arabia, Landslides, 13 (2016) 1315-1318. [33] L. Zaorálek, J. Platoš, V. Snášel, Patient-adapted and inter-patient ecg classification using neural network and gradient boosting, (2018). [34] F. Pan, T. Converse, D. Ahn, F. Salvetti, G. Donato, Feature selection for ranking using boosted trees, ACM, pp. 2025-2028. [35] R.G. Kumar, Y.S. Kumaraswamy, Investigating cardiac arrhythmia in ECG using random forest classification, International Journal of Computer Applications, 37 (2012) 31-34. [36] Q. Zhou, H. Zhou, T. Li, Cost-sensitive feature selection using random forest: Selecting low-cost subsets of informative features, Knowledge-based systems, 95 (2016) 1-11. [37] H. Ramchoun, M.A.J. Idrissi, Y. Ghanou, M. Ettaouil, Multilayer Perceptron: Architecture Optimization and Training, IJIMAI, 4 (2016) 26-30. [38] A. Subasi, Automatic recognition of alertness level from EEG by using neural network and wavelet coefficients, Expert systems with applications, 28 (2005) 701-711. [39] A. Güven, S. Kara, Classification of electro-oculogram signals using artificial neural network, Expert systems with applications, 31 (2006) 199-205. [40] J.J. Moré, The Levenberg-Marquardt algorithm: implementation and theory, Numerical analysis, Springer1978, pp. 105-116. [41] B. Karlik, A.V. Olgac, Performance analysis of various activation functions in generalized MLP architectures of neural networks, International Journal of Artificial Intelligence and Expert Systems, 1 (2011) 111-122. [42] R. Setiono, L.C.K. Hui, Use of a quasi-Newton method in a feedforward neural network construction algorithm, IEEE Transactions on Neural Networks, 6 (1995) 273-277. [43] S. Gomar, M. Mirhassani, M. Ahmadi, Precise digital implementations of hyperbolic tanh and sigmoid function, 2016 50th Asilomar Conference on Signals, Systems and Computers, IEEE, 2016, pp. 1586-1589. [44] P.d.B. Harrington, Sigmoid transfer functions in backpropagation neural networks, Analytical Chemistry, 65 (1993) 2167-2168. [45] A.M. Vukicevic, G.R. Jovicic, M.M. Stojadinovic, R.I. Prelevic, N.D. Filipovic, Evolutionary assembled neural networks for making medical decisions with minimal regret: Application for predicting advanced bladder cancer outcome, Expert systems with applications, 41 (2014) 8092-8100. [46] P. Sadowski, Notes on backpropagation, homepage: https://www. ics. uci. edu/pjsadows/notes. pdf (online), (2016).

[47] B.T. Volpe, J.E. Ledoux, M.S. Gazzaniga, Information processing of visual stimuli in an ‘extinguished’field, Nature, 282 (1979) 722. [48] J.J. Foxe, G.V. Simpson, S.P. Ahlfors, Parieto-occipital∼ 1 0Hz activity reflects anticipatory state of visual attention mechanisms, Neuroreport, 9 (1998) 3929-3933. [49] K. Hasenstab, C.A. Sugar, D. Telesca, K. McEvoy, S. Jeste, D. Şentürk, Identifying longitudinal trends within EEG experiments, Biometrics, 71 (2015) 1090-1100. [50] R. Schneider, S. Lau, L. Kuhlmann, S.J. Vogrin, M. Gratkowski, M.J. Cook, J. Haueisen, Matching pursuit based removal of cardiac pulse-related artifacts in EEG/fMRI, Universitätsbibliothek Ilmenau2016. [51] A. Nehlig, J.-P. Armspach, I.J. Namer, SPECT assessment of brain activation induced by caffeine: no effect on areas involved in dependence, Dialogues in clinical neuroscience, 12 (2010) 255. [52] C.-A. Park, C.-K. Kang, Y.-D. Son, E.-J. Choi, S.-H. Kim, S.-T. Oh, Y.-B. Kim, C.-W. Park, Z.-H. Cho, The effects of caffeine ingestion on cortical areas: functional imaging study, Magnetic resonance imaging, 32 (2014) 366-371. [53] B. Wittevrongel, E. Khachatryan, M.F. Hnazaee, E. Carrette, L. De Taeye, A. Meurs, P. Boon, D. Van Roost, M.M. Van Hulle, Representation of steady-state visual evoked potentials elicited by luminance flicker in human occipital cortex: An electrocorticography study, Neuroimage, 175 (2018) 315-326. [54] A. Beydoun, S.H. Schechter, W. Nasreddine, I. Drury, Responses to photic stimulation in patients with occipital spikes, Electroencephalography and clinical neurophysiology, 107 (1998) 13-17. [55] P. Kidmose, D. Looney, M. Ungstrup, M.L. Rank, D.P. Mandic, A study of evoked potentials from ear-EEG, IEEE Transactions on Biomedical Engineering, 60 (2013) 2824-2830. [56] H.A. Lamti, M.M.B. Khelifa, A.M. Alimi, P. Gorce, Influence of mental fatigue on P300 and SSVEP during virtual wheelchair navigation, 2014 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, IEEE, 2014, pp. 1255-1258. [57] W. Yijun, W. Ruiping, G. Xiaorong, G. Shangkai, Brain-computer interface based on the highfrequency steady-state visual evoked potential, IEEE, pp. 37-39. [58] P.F. Diez, V.A. Mut, E.M.A. Perona, E.L. Leber, Asynchronous BCI control using high-frequency SSVEP, Journal of neuroengineering and rehabilitation, 8 (2011) 39. [59] G. Hakvoort, B. Reuderink, M. Obbink, Comparison of PSDA and CCA detection methods in a SSVEP-based BCI-system, Centre for Telematics & Information Technology University of Twente, (2011). [60] Z. Yan, X. Gao, Functional connectivity analysis of steady-state visual evoked potentials, Neuroscience letters, 499 (2011) 199-203. [61] C.L. Colby, M.E. Goldberg, The updating of the representation of visual space in parietal cortex by intended eye movements, Science, 255 (1992) 90-92. [62] J.J. Todd, R. Marois, Capacity limit of visual short-term memory in human posterior parietal cortex, Nature, 428 (2004) 751. [63] J. Ding, G. Sperling, R. Srinivasan, Attentional modulation of SSVEP power depends on the network tagged by the flicker frequency, Cerebral cortex, 16 (2005) 1016-1029. [64] J.P. Pijn, J. Van Neerven, A. Noest, F.H.L. da Silva, Chaos or noise in EEG signals; dependence on state and brain site, Electroencephalography and clinical Neurophysiology, 79 (1991) 371-381. [65] A. Babloyantz, C. Lourenço, Brain chaos and computation, International Journal of Neural Systems, 7 (1996) 461-471. [66] A. Khasnobish, S. Bhattacharyya, A. Konar, D. Tibarewala, A.K. Nagar, A Two-fold classification for composite decision about localized arm movement from EEG by SVM and QDA techniques, The 2011 International Joint Conference on Neural Networks, IEEE, 2011, pp. 1344-1351. [67] M. Clubley, C. Bye, T. Henson, A. Peck, C. Riddington, Effects of caffeine and cyclizine alone and in combination on human performance, subjective effects and EEG activity, British Journal of Clinical Pharmacology, 7 (1979) U57-U63.

[68] X. Chen, Y. Wang, S. Gao, T.-P. Jung, X. Gao, Filter bank canonical correlation analysis for implementing a high-speed SSVEP-based brain–computer interface, Journal of neural engineering, 12 (2015) 046008. [69] Y.U. Zhang, G. Zhou, J. Jin, X. Wang, A. Cichocki, Frequency recognition in SSVEP-based BCI using multiset canonical correlation analysis, International journal of neural systems, 24 (2014) 1450013. [70] H. Wang, Y. Zhang, N.R. Waytowich, D.J. Krusienski, G. Zhou, J. Jin, X. Wang, A. Cichocki, Discriminative feature extraction via multivariate linear regression for SSVEP-based BCI, IEEE Transactions on Neural Systems and Rehabilitation Engineering, 24 (2016) 532-541. [71] Y. Zhang, E. Yin, F. Li, Y. Zhang, T. Tanaka, Q. Zhao, Y. Cui, P. Xu, D. Yao, D. Guo, Two-stage frequency recognition method based on correlated component analysis for SSVEP-based BCI, IEEE Transactions on Neural Systems and Rehabilitation Engineering, 26 (2018) 1314-1323. [72] I. Rejer, Ł. Cieszyński, Independent component analysis for a low-channel SSVEP-BCI, Pattern Analysis and Applications, 22 (2019) 47-62. [73] D.H. Krishna, I.A. Pasha, T.S. Savithri, Classification of EEG motor imagery multi class signals based on cross correlation, Procedia Computer Science, 85 (2016) 490-495. [74] M.E. Sharbaf, A. Fallah, S. Rashidi, EEG-based multi-class motor imagery classification using variable sized filter bank and enhanced One Versus One classifier, IEEE, pp. 135-140. [75] A. Zaitcev, G. Cook, W. Liu, M. Paley, E. Milne, Source localization for brain-computer interfaces, Brain-Computer Interfaces, Springer2015, pp. 125-153.

List of Figures and Tables List of Figures: Figure 1. Block diagram of the proposed method. Figure 2. Schematic of the computation process of structuring of 22 channel raw EEG signals of 6 volunteers into 4-dimensional (C×N×V×S) arrays, one each for non-caffeinated (N) and caffeinated (M) signals, where, C is the number of EEG channels (C = 22), N is the number of sample points in each EEG signal (N = 5120), V is the number of volunteers (V = 6) and S is the number of visual stimuli (S = 7). Figure 3. Schematic representation of the computation of SSVEP signals across 22 EEG channels Figure 4. Schematic representation of the computation of SSVEP signals for 7 photic stimuli

Figure 5. Progressive mapping of the brain. (a) Non-caffeinated, and (b) caffeinated EEG signals with 7 photic stimuli. Figure 6. Pre-processing of the EEG signals. (a) Raw EEG signal, (b) Band-limited EEG signal and (c) Moving mean filtered EEG signal from channel OZ of a volunteer for at 5 Hz (theta) visual stimulus Figure 7. (a) SSVEP activity before stimulus (non-caffeinated) on EEG channels, (b) SSVEP activity mapping on cortex (non-caffeinated), (c) SSVEP activity after stimulus (caffeinated) on EEG channels, (d) SSVEP activity mapping on cortex (caffeinated), (e) Percentage change observed in SSVEP activity before and after stimulus Figure 8. (a) SSVEP counts vs. stimuli frequency (non-caffeinated), (b) SSVEP counts vs stimuli frequency (caffeinated), (c) Changes in SSVEP activity observed for each of the 7 photic stimuli. Figure 9. (a) Mutual information test, (b) False nearest neighbor test, (c) Phase space plot, (d) Recurrence plot (color), (e) Recurrence plot (b/w) List of Tables: Table 1. Statistically predominant RQA features of all channels at 10 Hz visual stimuli Table 2. Classification accuracy of the 8 SSVEP predominant channel with respect to all the 7 visual stimuli Table 3. Comparison of the proposed method with some state of the art methods

EEG Signal Acquisition

Preprocessing

Detection of SSVEP Signals (CCA)

Feature Extraction (RQA)

Feature Reduction (t-test, CART, BT, RF)

Classification (MLP)

Figure 1. Block diagram of the proposed method

Figure 2. Schematic of the computation process of structuring of 22 channel raw EEG signals of 6 volunteers into 4-dimensional (C×N×V×S) arrays, one each for non-caffeinated (N) and caffeinated (M) signals, where, C is the number of EEG channels (C = 22), N is the number of sample points in each EEG signal (N = 5120), V is the number of volunteers (V = 6) and S is the number of visual stimuli (S = 7).

Figure 3. Schematic representation of the computation of SSVEP signals across 22 EEG channels

Figure 4. Schematic representation of the computation of SSVEP signals for 7 photic stimuli

Figure 5. Progressive mapping of the brain. (a) Non-caffeinated, and (b) caffeinated EEG signals with 7 photic stimuli.

Figure 6. Pre-processing of the EEG signals. (a) Raw EEG signal, (b) Band-limited EEG signal and (c) Moving mean filtered EEG signal from channel OZ of a volunteer for at 5 Hz (theta) visual stimulus

Figure 7. (a) SSVEP activity before stimulus (non-caffeinated) on EEG channels, (b) SSVEP activity mapping on cortex (non-caffeinated), (c) SSVEP activity after stimulus (caffeinated) on EEG channels, (d) SSVEP activity mapping on cortex (caffeinated), (e) Percentage change observed in SSVEP activity before and after stimulus

Figure 8. (a) SSVEP counts vs. stimuli frequency (non-caffeinated), (b) SSVEP counts vs. stimuli frequency (caffeinated), (c) Changes in SSVEP activity observed for each of the 7 photic stimuli.

Figure 9. (a) Mutual information test, (b) False nearest neighbor test, (c) Phase space plot, (d) Recurrence plot (color), (e) Recurrence plot (b/w)

Table 1. Statistically predominant RQA features of all channels at 3 Hz visual stimuli EEG channel

Methods

RQA features

Mean ± SD Without stimulus (B)

With Stimulus (A)

Predictor importance

p-value

T5

t-test CART BT RF

-AN RR RR

-1249.521 ± 88.591 5.199 ± 6.340 5.199 ± 6.340

-1245.313 ± 58.838 4.464 ± 5.336 4.464 ± 5.336

-1 1 1

-----

O1

t-test CART BT

-ENT DIV LMAX ENT

-2.299 ± 0.342 0.014± 0.015 368.333 ± 479.470 2.299 ± 0.342

-2.267 ± 0.322 0.013± 0.010 340.958± 479.268 2.267 ± 0.322

-1 1 1 1

------

BT RF

-ENT RR TT DIV

-2.378 ± 0.366 4.009 ± 4.543 4.125 ± 0.789 0.005 ± 0.005

-2.326 ± 0.288 4.685 ± 5.320 4.171± 1.000 0.007 ± 0.007

-1 1 1 1

------

OZ

t-test CART BT RF

-AN DET TT

-1264.348 ± 93.932 47.308 ± 18.271 4.047± 1.031

-1259.315 ± 66.506 44.787 ± 17.091 3.866± 0.814

-1 1 1

-----

T4

t-test CART BT RF

-RATIO DET TT

-22.634 ± 21.354 52.416 ± 23.264 4.477 ± 1.709

-23.705 ± 26.099 50.514 ± 23.743 4.362 ± 1.530

-1 1 1

-----

T6

t-test CART BT RF

-TT LAM AD

-3.967 ± 0.959 39.605 ± 31.801 5.348 ± 0.641

-3.990 ± 1.019 38.664 ± 34.084 5.136 ± 0.525

-1 1 1

-----

--

--

RF PZ

P4

t-test CART

t-test CART BT

O2

RF

DET LMAX DIV RR

49.812 ± 17.875 412.583 ± 496.047 0.013 ± 0.014 4.225 ± 6.049

44.605 ± 17.033 401.375 ± 484.625 0.009 ± 0.007 2.579 ± 2.553

1 1 1 1

t-test CART BT RF

-ENT LAM DET

-2.294 ± 0.380 50.036 ± 27.144 49.854 ± 19.302

-2.353 ± 0.234 46.641 ± 26.773 45.970 ± 20.673

-1 1 0.99

-----

Table 2. Classification accuracy of the 8 SSVEP predominant channel with respect to all the 7 visual stimuli (*color representation: : 85% > accuracy ≥ 80%,

: accuracy ≥ 90%,

: 90%> accuracy ≥85% ,

: 80% > accuracy ≥ 70%,

: 70% > accuracy ≥ 60%)

Table 3. Comparison of the proposed method with some state of the art methods SSVEP signal detection method

Accuracy (%)

Reference

FBCCA MsetCCA PCA+MLR Correlated component analysis CCA+RQA

91.95 90.00 95.00 80.00 97.00

Chen et al. [68] Zhang et al. [69] Wang et al. [70] Yangsong Zhang et al. [71] Proposed method

Highlights: • • • • •

The effect of caffeinated coffee on SSVEP signals was studied. Activated regions of the cortex were identified using canonical correlation analysis. Significant RQA features were identified and subsequently classified. SSVEP signal generation was not hampered in the parieto-occipital region. Interestingly, frontal lobe activity of the brain was significantly increased.