Applying machine learning techniques to detect the deployment of spatial working memory from the spiking activity of MT neurons

1.   Introduction

Machine-learning algorithms are well-known approaches with a wide range of applications in biomedical signal processing, such as classification ^[1], regression ^[2], and optimization ^[3]. In medical applications, diagnosis, recognition, or prediction are the most important tasks that can be handled using classification techniques ^[4]. In literature, a wide range of applications can be found based on the classification techniques such as emotion recognition ^[5], text classification ^[6], activity recognition ^[7], and epileptic seizure classification ^[8]. In most classification techniques, the classification performance strongly depends on the input features of the classification algorithms. Hence, feature extraction, the process of summarizing the data into some indexes, and feature selection, the process of finding the optimum combination of the extracted features, play important roles in classification approaches ^[9].

Among the diverse features that can be extracted from any data, fractal-based features are known as powerful nonlinear tools for measuring data complexity ^[10]. For example, it was shown that an individual's eye movements are closely related to fractal patterns ^[11]. Also, it was revealed that the fractality of the electroencephalography (EEG) signals is significantly reduced in schizophrenia ^[12]. Furthermore, it found that the memory content could increase the fractality of the EEG signals ^[13]. Another class of powerful features is the features derived from transformed data. For instance, the features based on the time-frequency transforms strongly help to predict atrial fibrillation ^[14], diagnosis of heart failure ^[15], and detection of mitral valve prolapse ^[16]. Such transforms are sometimes not directly used for feature extraction; however, they can be beneficial in the data preparation processes prior to feature extraction ^[17].

The study of the brain's cognitive functions, such as working memory and over and covert visual attention, are interesting areas of research. Thus, many studies have been conducted to investigate how these top-down functions influence neural spiking activity in different areas of the brain ^[18,19,20]. Nonetheless, different studies claimed that no memory-related modulation could be found in the neural spiking activity of the middle temporal (MT) cortex ^[21,22,23]. On the other hand, a recent study showed that working memory increased the fractionality of the firing rate signals ^[24]. Many studies have shown the neural correlate of working memory as an increase in the firing rate of neurons in the prefrontal ^[19,25,26], parietal ^[27], and visual ^[25] cortices; however, no such memory-related change in the average spiking activity of neurons has been reported in extrastriate cortex including V4 ^[19] and the middle temporal (MT) cortex ^[19,25] even after applying machine learning techniques ^[25]. In this paper, we will focus on the maintenance of spatial information, which has been shown to be sent directly from frontal eye field to the extrastriate cortex through feedback connections in form of persistent spiking activity ^[19]. Although spatial working memory does not increase the average firing rate of extrastriate neurons, it significantly enhances the sensitivity of individual neurons to incoming visual stimuli ^[26], alters the correlated activity of the population of neurons ^[18], increases the power of the local field potential (LFP) in the frequency band of alpha-beta, enhances the spike-phase coherency of MT responses in the same frequency range ^[28] and increases the fractionality of MT spiking activity ^[24].

As mentioned above, a recent study revealed that when the spiking activity of neurons is mapped to the fractal dimension feature space, the content of working memory can be captured ^[24]. In this study, we examined whether a set of linear and/or nonlinear features could reveal the deployment of spatial working memory from the spiking activity of neurons in the area MT. In this regard, we used two different learners, three feature selection algorithms, and two cross-validation methods to show the robustness of our results. The remaining parts of the paper are arranged as follows: Section 2 describes the studied data, extracted features, selection methods, classification algorithms, and classification assessment criteria. Section 3 presents the results, and Section 4 discusses the results and concludes the paper.

2.   Materials and methods

2.1.   Data description

All experimental procedures were performed under the National Institutes of Health Guide for the Care and Use of Laboratory Animals, the Society for Neuroscience Guidelines and Policies. The protocols for all experimental, surgical, and behavioral procedures were approved by the Montana State University Institutional Animal Care and Use Committee.

In this study, the spiking activity of 131 neurons (stored at 32 kHz), recorded in 11 sessions using electrode arrays, was used. These signals were recorded from the area MT (Figure 1a) of the two male macaque monkeys' brains (five and seven years old). The monkeys, already acquainted with carrying out the memory-guided saccade (MGS) task, were positioned on a customized chair in front of a monitor (24 inches with 144 Hz refresh rate) at a distance of 28.5 cm from their eyes. Initially, during a surgery in which the monkeys were anesthetized, the recording chambers were mounted on the monkeys' skulls in the MT area. The recordings of single electrodes were used to confirm that the chambers are mounted in a desired (MT) area of the monkeys' brains, During the task, the monkeys' heads were restrained, and they received juice as a reward through a syringe pump. This reward delivery, as well as the visual stimulus presentation procedures, were controlled using the MonkeyLogic toolbox in MATLAB software. Moreover, a photodiode was used to record the actual time of the visual stimulus incidence. Then the recorded data were digitalized with a sampling frequency of 32 kHz and stored.

Figure 1.  a) The area MT of the brain wherein the firing rate signals were recorded. b) Potential positions for the visual stimulus during the visual period of the MGS task. The red dots indicate the locations in the same hemifield (IN conditions), and the green dot indicates the location in the opposite hemifield (OUT condition) as the receptive field of the recorded neuron. The dashed line shows an imaginary receptive field of a sample neuron. c) The schematic representation of the MGS task, including the fixation (appearance of the FP), visual (appearance of visual stimulus), memory (disappearance of the visual stimulus), and saccade (after the disappearance of the FP) periods, as well as the corresponding response of a sample neuron. The average neural response during IN and OUT conditions are shown in green and pink, respectively.

DownLoad: Full-Size Img PowerPoint

The MGS task commenced with the appearance of a fixation point (FP) in the center of the monitor. In this period, called the fixation period, the monkeys are required to fixate at the FP for 1000 ms. After that, while the monkeys were fixating on the FP, a visual stimulus appeared in one of the positions, whether IN (same visual hemifield as the neuron's receptive field; displayed with red dots) or OUT (opposite hemifield relative to the RF of the recorded neuron; displayed with a green dot) conditions shown in Figure 1b. The stimulus remained for 1000 ms and then disappeared (visual period). As soon as the visual cue disappeared, the memory period started and lasted for 1000 ms. During the memory period, the monkeys had to keep gazing at the FP and memorize the location of the disappeared cue. Finally, the monkeys were obliged to make a saccade to the remembered location after the disappearance of the FP to receive a reward. The MGS phases are simply portrayed in Figure 1c, along with the recorded spiking activity for a sample MT neuron during the MGS task (more details can be found in see ^[18,19]).

2.2.   Data preparation

Data preparation plays a vital role in obtaining the best results for machine-learning approaches. Here, in the first step, the average spiking activity of individual neurons across trials was obtained. In the next step, the average firing rate signals of each neuron in IN conditions and OUT were considered for extracting the features. In order to eliminate the spiking activity related to the disappearance of the visual stimulus in the memory period, the first 400 ms of this period were trimmed for the feature extraction step. In the final step, the smoothed signals were considered for the following processes.

2.3.   Feature extraction

Extracting distinguishable features is an essential step in detecting the presence of the working memory. This section introduces the six most used fractal features and some statistical measures. Also, four frequently used transforms in signal processing are described.

2.3.1.   Features based on fractal dimension

Fractal dimension (FD) is an index of complexity that refers to an object's non-integer dimension geometrically ^[29]. FD can be obtained using different algorithms; however, in general, it can be obtained based on the number of blocks forming a pattern or covering the data graph. Here, six famous algorithms for calculating the FD are described.

Higuchi fractal dimension (HFD): HFD is an accurate box-counting method that is the main algorithm for calculating the FD of a graph ^[30]. Any time series ${x}_{t}:x\left(1\right), x\left(2\right), \dots , x\left(N\right)$ with a finite number of samples ($N$) can be expressed as the $k$ sets of ${x}_{m}^{k}$ where

Here, $m = \mathrm{1, 2}, \dots , k$ is the initial time, $k = \mathrm{1, 2}, \dots , {k}_{m}$ is the delay, and $R = \left[\frac{N-m}{k}\right]$. Note that $\left[\dots \right]$ denotes the integer part of the internal number. Accordingly, the curve length of each subset can be defined as

Letting ${\left\langle{{L}^{k}}\right\rangle}_{m}$ defines the average value of $m$ curve lengths, it can be written that

where $D$ is FD obtained using the Higuchi algorithm. In this paper, ${k}_{m} = 30$ was determined based on trial and error.

Katz fractal dimension (KFD): KFD is a distance-based technique for obtaining the FD value that uses the averaged Euclidean distance of two consecutive points in a time series ^[31]. If $L$ defines the sum of Euclidean distance between every two successive samples ($L = \sum \nolimits_{i = 1}^{N-1}dis{t}_{Euclidean}\left({n}_{i}, {n}_{i+1}\right)$) and ${L}_{m}$ denotes the maximum value of the Euclidean distance between the first and $j$th sample (${L}_{m} = \mathit{max}\left(dis{t}_{Euclidean}\left({n}_{1}, {n}_{j}\right)\right)\;for\;j = 2, \dots , N$), the FD can be obtained as

In the above equation, $N$ is the number of samples in a given time series, and $D$ is the FD value with the Katz technique.

Generalized Hurst exponent (GHE): GHE is known as a fractal index indicating the long-rage dependence of time series that uses the $q$th-order moment of the distribution ^[32]. Assuming $x\left(t\right)$ as the studied time series with $N$ samples, the the $q$th-order moment of the distribution can be defined as

where $\tau$ is the delay and ${\left\langle{\dots }\right\rangle}_{t}$ indicates the average value of the internal value over the total duration. Then, the fractal dimension of the time series can be obtained through

where ${T}_{s}$ is the sampling time, and $H$ refers to the GHE of the time series. In this study, the first-order moment was considered.

Margaos and Sun fractal dimension (MSFD): MSFD is a morphological method for calculating FD since it tries to cover the data graph employing the morphological operators, i.e., erosion and dilation ^[33]. The support-limited erosions and dilations with the support set $s$ ($s = \mathrm{1, 2}, \dots , N$) and using the structuring element of $b$ can be formulated as

where $k = \mathrm{1, 2}, \dots , {k}_{m}$. Consequently, the morphological cover that surrounds the data graph can be obtained as

Here ${k}_{m}$ is set according to the rule mentioned in ^[34]. Finally, the morphological FD of the time series can be acquired as the angular coefficient of linear regression of $\mathit{ln}\left(\frac{{C}_{k}}{{\left(\frac{2k}{N}\right)}^{2}}\right)$ vs. $\mathit{ln}\left(\frac{1}{\frac{2k}{N}}\right)$.

Leibovich and Toth fractal dimension (LTD): LTD is the fast implementation of the box-counting algorithm ^[35]. The box-counting algorithm is based on the number of blocks to which the data graph can be split. So, the FD can be found as

where $d$ is the size of the blocks and $D$ is the FD base on the original box-counting method. The LTD algorithm implements the box-counting method by withdrawing too small and too large blocks; therefore, it is faster than the original algorithm ^[36].

Fractal Volatility (FV): FV computes the FD of a time series using the box-counting method with the rand-walk process ^[37]. In other words, it splits the data into blocks of size $d$ by performing the random-walk process.

2.3.2.   Features based on transforms

Discrete Wavelet Transform (DWT): DWT describes any data by a set of weighted orthonormal wavelets called mother wavelets ^[38] according to the following relation

where $N$ is the total number of samples, ${X}_{m, n}$ is the transformed data, $\psi$ is the orthonormal wavelet. Moreover, $m$ and $n$ are two control parameters regarding the translation and dilation operation of the discrete wavelet and the data. DWT includes temporal information as well as frequency information in different scales. In this paper, the db4 was selected as the orthonormal wavelet for further analysis.

Discrete Fourier Transform (DFT): DFT is the most fundamental transformation that describes the data based on some weighted sinusoidal functions. DFT can be defined

where $k = \mathrm{1, 2}, \dots , N$ and ${X}_{k}$ is the transformes data for a specific frequency.

Discrete Short-Time Fourier Transform (DSTFT): DSTFT is the time-frequency version of the DFT transform that contains only the frequency information. Thus, DSTFT can be helpful when DFT becomes insufficient for analysis, such as nonstationary signals ^[39]. DSTFT is defined as

where $\omega \left(L\right)$ is the temporal window of size $L$.

Discrete Stockwell Transform (DST): DST is the extension of the DWT while having a close relationship with DSTFT. Also, it provides a frequency-dependent resolution since the sinusoidal functions are fixed in time, and a scalable Gaussian window operates the dilation and translation ^[40]. DST can be described as

where ${T}_{s}$ is the inverse sampling frequency, $i{T}_{s}$ defines the window $\tau$, $\frac{n}{N{T}_{s}}$ is the frequency domain, and $H\left(n\right)$ is the DFT of the input data.

2.3.3.   Moments of distribution

Moments of distribution are basically statistical measures, the primary features that can be simply obtained from any data. The first-, second-, third-, and forth-order moments are called mean, variance, skewness, and kurtosis, which are described as follows

It should be noted that $N$ is the number of data samples.

Median is another statistical measure that can be helpful whenever the mean of that data is not a good measure of the distribution. The median can be defined as follows

where $k = \mathrm{1, 2}, \dots , N$. Maximum and minimum values of the samples within a specific period are the other two famous statistical measures used in this paper.

2.4.   Feature selection

In signal processing, feature selection is an optimization method leading to the optimum features for classification. Therefore, feature selection was used as a preprocessing step for machine-learning problems and is of particular importance when the data or extracted features are of high dimension ^[41]. In this subsection, three popular algorithms, namely, Genetic Algorithm (GA), Particle Swarm Optimization (PSO), and Ant Colony Optimization (ACO) algorithms, are briefly described.

2.4.1.   Genetic Optimization

GA is an evolutionary algorithm that can be applied for search and optimization problems inspired by species' natural selection and evolution. In this algorithm, a population of chromosomes, which are, in fact, the solution to the problem, are directly selected for the next generation. Based on the selected species, the new generation, including a few numbers of new chromosomes, is also created using crossover and mutation techniques. This process is heuristically repeated until the solution converges to the optimal value ^[42]. In this article, five primary chromosomes and 100 repetitions with a crossover rate of 0.6 and a mutation rate of 0.001 were selected as the initial parameters of the GA algorithm for selecting the optimum features.

2.4.2.   Particle Swarm Optimization

PSO is an evolutionary computational algorithm that tries to find an optimized solution inspired by natural social behaviors, such as the group behavior observed in schools of fish and flocks of birds. According to this idea, some initial particles are localized in the search space, each having its own position and velocity. These particles move within the search space to find the best solution. In general, particles make their next moves based on their best-experienced position as well as the best-experienced position of the whole population of particles called a swarm. So, the movement of the particles in the search space is not dependent on the gradient, and therefore, it can be applied for differentiable and non-differentiable optimization problems ^[43]. In this paper, we set cognitive and social factors ${C}_{1} = {C}_{2} = 2$ and inertia weight $w = 1$ with five initial particles and 100 maximum iterations as the assumed parameters for the PSO algorithm used for obtaining the optimum features.

2.4.3.   Ant Colony Optimization

ACO is a probabilistic and graph-based evolutionary algorithm that was first proposed based on the natural behavior of ants in hunting their food. In real life, ants leave pheromones on the pass to guide others to resources as they explore their surroundings, and thus, they can find the shortest pass to the food. Inspired by this cooperative-based technique, optimization problems can be solved and handled. First, some initial artificial ants are positioned in the parameter space that move to a solution stochastically. The pheromone trails, which specify the edges of the graph in the ACO algorithm, are obtained for each ant, and the best solution is selected. For the following steps, the edges of the graph become updated and guide the artificial ants toward the solution. This process is repeated iteratively until the solution converges to an optimal value ^[44]. In this paper, five initial artificial ants, 100 allowed iterations, $\alpha = \tau = \eta = 1$, $\rho = 0.2$, and $\beta = 0.1$ are selected as the parameters needed for applying the ACO algorithm to find optimum features.

2.5.   Feature classification

Feature classification is the process of assigning or predicting the label of new data based on the trained model or information gained from the observed data ^[45]. The Support Vector Machine (SVM) and the K-Nearest Neighbor (KNN) classifiers are two supervised and most-used machine-learning algorithms for classification, which are briefly described below.

2.5.1.   Support Vector Machine classifier

The original SVM algorithm can be used for classifying two classes using a linear boundary; however, it has been extended for the classification of multi-class data. The SVM classifier builds an optimum hyperplane with the largest margin—with the maximum distance from the nearest data to the decision boundary—that can distinguish the data with the highest accuracy. If the data is not linearly distinguishable, using nonlinear kernels, the data will be implicitly mapped into the higher dimension wherein the data is distinguishable with a linear boundary ^[44]. This paper uses the SVM classifier with a three-order polynomial kernel function to classify the neuronal spiking activity in fixation and memory periods.

2.5.2.   K-Nearest Neighbor classifier

The KNN algorithm classifies each new data based on voting between the class of k-closest observed data to the new sample. Therefore, to classify any new data, the KNN classifier needs to find the k-nearest training data by computing the distances of the new sample from all other samples in the parameter space. Therefore, although the KNN algorithm is simple to implement, it may be time-consuming due to its computational costs. This paper employs the KNN classifier with $k = 3$ and standardized Euclidean distance to detect the presence of working memory using the firing rate data.

2.5.3.   Classification evaluation

To assess the classification performance, different criteria are introduced in the literature ^[16]. Accuracy, the most well-known evaluation criterion for classification, is defined based on the number of samples correctly labeled by the classifier versus the number of samples (see Eq (19)). Sensitivity and specificity are more specific assessment criteria since they respectively show the performance of the classifier in detecting and not detecting the target class. Therefore, sensitivity is defined as the number of samples correctly labeled as the target class by the classifiers versus the actual number of samples in the target class (see Eq (20)). On the other hand, specificity is defined as the number of samples correctly labeled as the non-target class by the classifiers versus the actual number of samples in the non-target class (see Eq (21)).

Here, $tp$, $tn$, $fp$, and $fn$ refer to the true positive, true negative, false positive, and false negative indexes.

3.   Results

Figure 2 shows the average normalized response of 131 MT neurons during the MGS task (see Methods). According to this figure, on average, no change in spiking activity is observed when comparing the response of neurons before and after the visual stimulus (i.e., fixation and memory periods, respectively). The inset bar graph in Figure 2 reveals no significant difference between the average spiking activity of MT neurons in memory vs. fixation periods in both IN and OUT memory conditions (${p}_{fixation\;IN\;vs, memory\;IN} = 0.385$ and ${p}_{fixation\;OUT\;vs, memory\;OUT} = 0.385$).

Figure 2.  The normalized response of 131 MT neurons during the MGS task, including the fixation, visual, memory, and saccade periods. The average response in IN and OUT conditions are shown in green and pink, respectively. The inset bar graphs show the average firing rate of MT neurons in IN and OUT conditions during the fixation and memory periods.

DownLoad: Full-Size Img PowerPoint

To detect the presence of the memory, we used the neural spiking activity in IN conditions during the fixation (when no working memory is involved) and memory (when deployment of top-down working memory signals is present) periods for feature extraction, selection, and classification. In the feature extraction step, 41 features are extracted from the IN conditions in fixation and memory periods. Therefore, the feature vector can be described as:

- Index 1–6 are the fractal-based features, including HFD, KFD, GHE, MSFD, LTFD, and FV.

- Index 7–34 are the transform-based features, including mean (index 7–10), variance (index 11–14), kurtosis (index 15–18), skewness (index 19–22), median (index 23–26), maximum (index 27–30), and minimum (index 31–34) of the components of DWT, DFT, DSTFT, and DST, respectively.

- Index 35–41 are statistical features, including mean, variance, kurtosis, skewness, median, maximum, and minimum values of the firing rate signals within the fixation and memory periods.

Finally, the selected features from the fixation IN and memory IN were classified using the SVM and KNN classifiers. Based on the selection method, classification was performed in four cases. To estimate the performance of the classifiers on the data, the k-fold cross-validation method (with $k = 10$) and its balanced version called the A-test algorithm (with $k = 10$ and ten iterations) were employed as the procedure for performing the classification. Unlike the K-fold cross-validation method, the A-test algorithm ensures that the same number of samples from each class is included in all training and testing folds. Therefore, the A-test method might be more reliable, especially when classes contain unequal data samples.

3.1.   All extracted features combined

In the first case, all extracted features, including 41 linear and nonlinear features, were used for the classification step. The results of the classification performance are shown in Figure 3. The details of the results can be found in Table 1.

Figure 3.  The average performance of the classifiers SVM and KNN using the K-fold (with $K = 10$) and A-test (with $K = 10$ and 10 iterations) cross-validation methods. In this case, the classification was performed using all extracted features.

DownLoad: Full-Size Img PowerPoint

Table 1.  Classification performance results ($mean\%\pm STD\%$) of fixation IN versus memory IN using all extracted features.

Classifier	Assessment criterion	K-fold ($\boldsymbol{K}={\bf{10}}$)	A-test ($\boldsymbol{K}={\bf{10}}, \boldsymbol{i}\boldsymbol{t}\boldsymbol{e}\boldsymbol{r}={\bf{10}}$)
SVM	Accuracy	${\bf{98.85\pm 2.60}}$	$98.77\pm 0.16$
	sensitivity	$97.69\pm 5.19$	$97.54\pm 0.32$
	specificity	$100\pm 0$	$100\pm 0$
KNN	Accuracy	${\bf{98.46\pm 1.99 }}$	$97.92\pm 0.41$
	sensitivity	$97.69\pm 3.72$	$97\pm 0.67$
	specificity	$99.23\pm 2.43$	$98.85\pm 0.54$

---

 | Show Table

DownLoad: CSV

Figure 3 shows that the SVM and KNN classifiers SVM and KNN classifiers performed closely ($Accurac{y}_{SVM} = 98.85\pm 2.6$; $Accurac{y}_{KNN} = 98.46\pm 1.99$) in distinguishing fixation IN and memory IN data. However, SVM performed slightly better than KNN ($Accurac{y}_{SVM}-Accurac{y}_{KNN} < 1\mathrm{\%}$) in both cross-validation approaches. Table 1 reveals that the average accuracy is higher in the K-fold approach; however, the standard deviation mentioned in the A-test method is considerably lower. This shows that in different iterations, the average accuracy is not changed remarkably, and thus, the results are valid.

3.2.   GA selecting method

Here, 20 features were selected using the GA feature selection method, including three fractal-based features (including the HFD, the KFD, and the FV), 14 transform-based features (including the mean of DFT, the mean of DST, the variance of DFT, the variance of DSTFT, the variance of DST, the kurtosis of DWT, the kurtosis of DSTFT, the skewness of DFT, the skewness of DSTFT, the median of DWT, the median of DFT, the median of DST, the maximum of DWT, and the minimum DWT), and three statistical features (including the skewness, the median, and the minimum values). Figure 4 and Table 2 demonstrate the classification results using the 20 GA-based selected features.

Figure 4.  The average performance of the classifiers SVM and KNN using the K-fold (with $K = 10$) and A-test (with $K = 10$ and 10 iterations) cross-validation methods. In this case, the classification was performed using the feature selected by GA.

DownLoad: Full-Size Img PowerPoint

Table 2.  Classification performance results ($mean\%\pm STD\%$) of fixation IN versus memory IN using the GA as the feature selection method.

Classifier	Assessment criterion	K-fold ($\boldsymbol{K}={\bf{10}}$)	A-test ($\boldsymbol{K}={\bf{10}}, \boldsymbol{i}\boldsymbol{t}\boldsymbol{e}\boldsymbol{r}={\bf{10}}$)
SVM	Accuracy	$98.46\pm 2.69$	${\bf{99.04\pm 0.33}}$
	sensitivity	$96.92\pm 5.38$	$98.08\pm 0.65$
	specificity	$100\pm 0$	$100\pm 0$
KNN	Accuracy	${\bf{99.23\pm 1.62}}$	$99.19\pm 0.12$
	sensitivity	$98.46\pm 3.24$	$98.38\pm 0.24$
	specificity	$100\pm 0$	$100\pm 0$

---

 | Show Table

DownLoad: CSV

According to Figure 4, unlike Figure 3, the KNN classifier performed slightly more effectively than the SVM classifier ($Accurac{y}_{KNN}-Accurac{y}_{SVM} < 1\mathrm{\%}$). Table 2 shows that the GA-selected feature improved the classification performance, especially for the KNN classifier ($Accurac{y}_{KNN} = 99.23\pm 1.62$). In addition, the KNN classifier has more reliable results since it generally has a lower standard deviation.

3.3.   PSO selecting method

Employing the PSO feature selection method, 17 features were selected, including three fractal-based features (including the HFD, the KFD, and the FV), 12 transform-based features (including the mean of DFT, the mean of DST, the variance of DFT, the variance of DSTFT, the kurtosis of DWT, the kurtosis of DSTFT, the skewness of DFT, the median of DWT, the median of DFT, the median of DST, the maximum of DWT, the minimum of DWT), and two statistical features (including the variance, and the maximum values). Similar to the previous subsections, Figure 5 and Table 3 illustrate the classification performance considering the PSO-based selected features.

Figure 5.  The average performance of the classifiers SVM and KNN using the K-fold (with $K = 10$) and A-test (with $K = 10$ and 10 iterations) cross-validation methods. In this case, the classification was performed using the feature selected by PSO.

DownLoad: Full-Size Img PowerPoint

Table 3.  Classification performance results ($mean\%\pm STD\%$) of fixation IN versus memory IN using the PSO as the feature selection method.

Classifier	Assessment criterion	K-fold ($\boldsymbol{K}={\bf{10}}$)	A-test ($\boldsymbol{K}={\bf{10}}, \boldsymbol{i}\boldsymbol{t}\boldsymbol{e}\boldsymbol{r}={\bf{10}}$)
SVM	Accuracy	$99.23\pm 2.43$	${\bf{99.50\pm 0.26}}$
	sensitivity	$98.46\pm 4.87$	$99.00\pm 0.52$
	specificity	$100\pm 0$	$100\pm 0$
KNN	Accuracy	${\bf{98.85\pm 1.86}}$	$98.77\pm 0.24$
	sensitivity	$97.69\pm 3.72$	$97.77\pm 0.24$
	specificity	$100\pm 0$	$99.77\pm 0.37$

---

 | Show Table

DownLoad: CSV

In contrast to Figure 3, in this case, Figure 5 shows that the PSO-selected features helped the SVM classifier improve its performance ($Accurac{y}_{SVM} = 99.50\pm 0.26$) more considerably than the KNN classifier ($Accuracy = 98.85\pm 1.86$). Besides, Table 3 reveals that the KNN classifier has more reliable results due to the lower standard deviation on average.

3.4.   ACO selecting method

Using the ACO algorithm for selecting the optimum features, a total of 20 features were selected, including three fractal-based features (including the HFD, the GHE, and the FV), 14 transform-based features (including the skewness of DSFTF, the variance of DWT, the kurtosis of DSTFT, the skewness of DWT, the kurtosis of DFT, the mean of DST, the skewness of DFT, the min of DFT, the mean of DWT, the variance of DSTFT, the variance of DST, the maximum of DFT), and five statistical features (including the skewness, the maximum, the variance, the kurtosis, and the minimum values). Figure 6, as well as Table 4, contains the results of performing the classification using the ACO-based selected features.

Figure 6.  The average performance of the classifiers SVM and KNN using the K-fold (with $K = 10$) and A-test (with $K = 10$ and 10 iterations) cross-validation methods. In this case, the classification was performed using the feature selected by ACO.

DownLoad: Full-Size Img PowerPoint

Table 4.  Classification performance results ($mean\%\pm STD\%$) of fixation IN versus memory IN using the ACO as the feature selection method.

Classifier	Assessment criterion	K-fold ($\boldsymbol{K}={\bf{10}}$)	A-test ($\boldsymbol{K}={\bf{10}}, \boldsymbol{i}\boldsymbol{t}\boldsymbol{e}\boldsymbol{r}={\bf{10}}$)
SVM	Accuracy	$99.23\pm 1.62$	${\bf{99.42\pm 0.20}}$
	sensitivity	$98.46\pm 3.24$	$98.85\pm 0.41$
	specificity	$100\pm 0$	$100\pm 0$
KNN	Accuracy	$99.62\pm 1.22$	${\bf{99.65\pm 0.12}}$
	sensitivity	$99.23\pm 2.43$	$99.31\pm 0.24$
	specificity	$100\pm 0$	$100\pm 0$

---

 | Show Table

DownLoad: CSV

Figure 6 illustrates that using the ACO-selected features, the average performance of both classifiers in detecting the presence of working memory is significantly enhanced compared to Figure 3, wherein all features were involved in the classification procedure ($Accurac{y}_{SVM} = 99.42\pm 0.20$; $Accurac{y}_{KNN} = 99.65\pm 0.12$). From Table 4, it can be seen that the KNN classifier not only has slightly better performance ($Accurac{y}_{KNN}-Accurac{y}_{SVM} < 0.5\mathrm{\%}$) but also, due to the lower standard deviation, has a more reliable performance.

4.   Discussion

The brain is the most complex system in the human body. This complexity is reflected in the signals recorded from the brain. Thus, brain-associated data such as EEG or the spiking activity of neurons predominantly have nonlinear properties. This nonlinearity can be captured by FD, which is an index of complexity. To obtain FD of a time series, different algorithms have been proposed, such as HFD ^[30], KFD ^[31], GHE ^[32], MSFD ^[33], LTD ^[35], and FV ^[37]. However, transform-based features and statistical indexes are popular in signal processing. The main objective of this paper was to examine the ability of machine-learning methods to detect the presence of working memory using various linear and nonlinear features. Therefore, we used several algorithms to obtain the FD value (including HFD, KFD, GHE, MSFD, LTFD, and FV) and different transforms to obtain the frequency and/or time-frequency components (including statistical measures of DWT, DFT, DSTFT, and DST) of the average spiking activity of MT neurons. Also, we included some of the important statistical indexes (including mean, variance, kurtosis, skewness, median, maximum, and minimum values) in the feature set.

Selecting the optimum features can be an optimization problem that, in most cases, results in improving the classifiers' performance. Such methods mainly focus on finding the best set of features that can lead to the best classification result. In this way, it can also help reduce the feature space's dimensionality and simplify the classification problem. Hence, after performing the classification with all extracted features, we examined whether the feature selection method could enhance the classification performance. Accordingly, we used GA, PSO, and ACO algorithms to select the optimum features for detecting the presence of working memory. In the classification step, the ability of two machine-learning algorithms, namely SVM and KNN classifiers, were employed to detect the presence of memory. KNN is typically considered a nonlinear classifier as it has a nonlinear decision boundary whilst SVM can be a nonlinear classifier if it uses a nonlinear kernel function. In general, the nonlinearity of the decision boundary enables a classifier to learn and distinguish the data classes more precisely and define the membership's probability to each data class for new data. For this reason, we used SVM with a three-order polynomial kernel function and KNN with three nearest neighbors. It should be noted that cross-validation methods can determine the validity of classification results, particularly when the number of samples is not too high. Therefore, K-fold (with 10 folds) and A-test (with 10 folds and 10 iterations) cross-validation methods were performed to show the results' robustness.

The best classification performance of the SVM and KNN classifiers is summarized in Figure 7. According to this figure, when no selection algorithms were employed, the SVM classifiers led to the best performance ($Accurac{y}_{SVM} = 98.85\pm 2.60$). The same result can be seen in the case where the PSO algorithm was used to select the features ($Accurac{y}_{SVM} = 99.50\pm 0.26$). In contrast, when the GA- and ACO-selected features were used as the input of the classification algorithms, the KNN classifier reached higher average accuracy ($Accurac{y}_{KNN} = 99.23\pm 1.62$ and $99.65\pm 0.12$, respectively). Moreover, Figure 7 shows that employing the selection method can improve the classification performance since, in all cases, the average accuracy was grown compared to the case where no selecting method was employed.

Figure 7.  The best performance of KNN and SVM classifiers employing no selection algorithm as well as the GA, PSO, and ACO selection methods.

DownLoad: Full-Size Img PowerPoint

In total, Figure 7 reveals that when the features were selected using the ACO algorithm, the KNN classifier could detect the presence of working memory with the accuracy of $99.65$% and the standard deviation of $0.12$ (using the A-test cross-validation method), which is the highest obtained average accuracy among the studied cases. In this case, three out of six FD-base features (including the HFD, the GHE, and the FV), 14 out of 28 transform-based features (including the skewness of DSFTF, the variance of DWT, the kurtosis of DSTFT, the skewness of DWT, the kurtosis of DFT, the mean of DST, the skewness of DFT, the min of DFT, the mean of DWT, the variance of DSTFT, the variance of DST, the maximum of DFT), and five out of seven statistical-based features (including the skewness, the maximum, the variance, the kurtosis, and the minimum values) were involved in the classification procedures.

To obtain the neural code for spatial working memory represented in the firing rate of visual neurons, we compared the neural responses of MT neurons during the memory period (where the monkey is actively memorizing a location) with the neural activity during the fixation period (when no working memory is present). This approach could be questioned as one can argue that any differences between the neural responses during memory and fixation periods could occur due to other cognitive signals such as arousal or expectation than spatial working memory. Here we review a series of neurophysiological evidence revealing the dependence of these response changes (i.e., the differences of neural responses between memory and fixation period) on the content of working memory:

- By measuring the neural responses of extrastriate neurons to visual probes presented during fixation and memory periods, a recent study showed a strong modulation of RF profile in V4 and MT neurons that were dependent on the content of working memory ^[19]. It was shown that during the maintenance of spatial information, only the neurons whose RFs during the fixation period were close to the remembered location, expanded and shifted their RFs towards that location during the memory period.

- It was also found that the encoding of the visual probe's location by the population of MT neurons was enhanced during the memory period compared to the fixation period. In detail, this was measured by the ability of MT individual neurons' firing activity to discriminate two different visual probes (two-point discriminability). In fact, this memory-related enhancement in two-point discriminability only occurred to those visual probes that were presented near the locus of working memory.

- At the level of LFP, it has been shown that the amount of information regarding the visual input conveyed by the alpha-beta phase of spike times increases during the memory period compared to the fixation period. This phenomenon only occurred to those visual probes which were presented near the remembered location.

- Furthermore, it was shown that the discrimination between visual probes, based on the phase of each spike in the alpha-beta frequency range, is enhanced during the memory period compared to fixation. This discrimination enhancement was observed for the visual probes presented near the locus of working memory.

As the differences between the neural responses during memory and fixation periods occur in a spatial-specific manner (i.e., near the locus of working memory), it would be very unlikely to relate these response changes to any other cognitive signals, such as arousal, than spatial working memory.

Acknowledgments

Yaser Merrikhi collected this dataset in the laboratory of Dr. Behrad Noudoost at Montana State University, Bozeman, MT, USA. The lab's experiments were supported by MSU start-up fund, Whitehall 2014-5-18, NIH R01EY026924, and NSF143221 and 1632738 grants to Dr. Noudoost. We would like to thank Dr. Noudoost for sharing this dataset. This work is funded by the Centre for Nonlinear Systems, Chennai Institute of Technology, India, vide funding number CIT/CNS/2022/RP-006.

Conflict of interest

The authors declare there is no conflict of interest.