Higuchi fractal dimension is a unique indicator of working memory content represented in spiking activity of visual neurons in extrastriate cortex

1.   Introduction

The primary function of working memory brain areas is to help store information in one's mind while doing different and complicated tasks in parallel ^[1]. In fact, in cognitive science, through storing information, working memory is assumed to be the basis of high-level activities such as decision-making and recognition, learning, overt and covert attention, social interactions and oral communication ^[2,3,4]. In psychology and cognitive science, working memory and attention are considered to be heavily interconnected as their malfunction has been widely reported in many cognitive disorders ^[5,6], and also their functional interdependence would lead to appropriate psychophysical outcomes ^[7].

Considering the significant role of working memory in the attention process, many studies have been performed in order to understand the underlying neural mechanisms of working memory. For instance, memory-guided saccade (MGS) task, which is able to activate the working memory process in the brain, has been widely used in many and electrophysiological and fMRI studies ^[8,9,10,11]. As reported in ^[8,9,12], working memory can be represented as an increase in spike rates recorded in brain areas, including the medial superior temporal, lateral prefrontal cortex, frontal eye field, and lateral intraparietal area. In the occipital cortex, studies found that memorizing a visual object (not its spatial information) increases the firing rate of V4 neurons ^[13]. However, there is no sign of spatial working memory in the neurons' spiking activities in the early and mid-level visual areas, such as middle temporal (MT) cortex and area V4 ^[14].

In signal processing, different studies have been conducted to identify features characteristic of fundamental processes or diseases ^[15,16,17]. Fractal-based and entropy-based features are important nonlinear features that measure the complexity of data ^[18] and are widely used in signal processing for different purposes.

Although the spiking activity of MT neurons shows no increase in their average firing rate during working memory ^[14], a recent study claimed that the dimensionality of neurons' spiking activity increases during the deployment of spatial working memory ^[19]. This paper is devoted to analyzing the ability of nonlinear and classical features to detect memory content. The following parts of the paper are organized as follows: Section 2 describes the data used in this study and its acquisition process. Section 3 introduces the nonlinear (including fractal- and entropy-based) and classical (including transform- and statistical-based) features in brief. Section 4 presents the study's results in detail and discusses the most important findings. Finally, Section 5 concludes the paper.

2.   Data

The data used in this study were gathered from the spiking activity of 131 neurons within area MT. To record the firing rate signals, a five- and a seven-year-old male macaque monkeys, which had been already trained to perform the MGS task, were employed. The signals were recorded using an electrode with 16 arrays (recoding sites) over 11 recording sessions. The recording chambers were mounted on the skull of two monkeys' over area MT during the anesthetized surgical procedures. Before data collection, the site of brain areas within the recording chambers was confirmed through single-electrode recordings. As the initial settings before starting the task and data acquisition, the monkeys were seated in a custom-built chair 28.5 centimeters away from a monitor of 24 inches and 144 Hz refresh rate. While their heads were restrained, a syringe pump was set in their mouth to deliver juice as a reward (see ^[8,10] for more details).

Note that the experimental procedures were performed in accordance with 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. Figure 1 demonstrates the four main stages of the MGS task:

Figure 1.  a) four phases of the MGS task: the appearance of FP for 1000 ms (fixation period), wherein the monkeys should gaze at the FP; the presented visual cue for 1000 ms (visual period), wherein the monkey should keep looking at the FP; the disappearance of the visual cue for 1000 ms (memory period), wherein the monkeys remained fixated at the FP; the disappearance of the FP (saccade period), wherein the monkeys should make a saccade toward the memorized location. b) In the visual period, the cue can occur in one of the six positions, including five IN conditions (red dots) and one OUT condition (green dot). The dashed line indicates the schematic receptive field of a sample neuron.

DownLoad: Full-Size Img PowerPoint

- Fixation period: This phase starts as soon as a fixation point (FP) appears in the middle of the monitor. The monkeys should fixate on the FP on the monitor for one second.

- Visual period: In this phase, a cue appears either in the same visual hemifield (one of the red dots in Figure 1b; IN conditions) or in the opposite one (green dot in Figure 1b; OUT condition), while the monkeys are needed to keep gazing at the FP. This phase takes a second, and then the cue disappears.

- Memory period: After the cue goes off, the monkeys are still required to hold their gaze on the FP for another one second. Since the process of memorizing the cue's location occurs during this phase, it is called the memory period.

- Saccade period: In this phase, the FP goes off, and the monkeys perform a saccade toward the recalled cue's locations within 500 ms.

During the task, the actual time of the visual cue's occurrence in the monitor was recorded using a photodiode, and the signals were digitalized at a 32 kHz sampling frequency. Moreover, the presentation of the cue on the monitor and reward delivery was monitored using the MonkeyLogic toolbox in MATLAB software. Electrode arrays were employed to record the spiking activity of 131 MT neurons while the animals performed the task. Note that the recorded data were digitalized and stored at 32 kHz for further analysis.

Figure 2a demonstrates the averaged response of a sample neuron's spiking activity during the MGS task, and Figure 2b shows the normalized response of 131 neurons within the area MT. The purple signal refers to the IN condition, and the blue one corresponds to the OUT condition. As the cue appeared in the visual period, the firing rate of neurons increased transiently. Similarly, another transient increase in the spiking activity of MT neurons was observed in response to the disappearance of the visual stimulus. However, no such effect was seen for the OUT condition as the visual stimulus was presented outside the neuron's receptive field. The same scenario is seen in the average normalized response of 131 MT neurons.

Figure 2.  a) The average response of a sample neuron during the MGS task and b) the normalized response averaged over 131 neurons. In both sections, the purple plots are the averaged response during IN conditions, and the blue plots are the OUT condition's response.

DownLoad: Full-Size Img PowerPoint

3.   Methods

In this paper, we examined the ability of nonlinear and classical features to detect the content of working memory from the neural spiking activity of the extrastriate cortex. Therefore, two main categories of features, including nonlinear and classical features, were studied.

3.1.   Nonlinear features

Fractal dimension (FD) and entropy are two main categories of nonlinear features that quantitatively measure the complexity of a time series ^[20,21]. Table 1 demonstrates the definition of the six most used algorithms, namely Higuchi FD (HFD), Katz FD, Generalized Hurst exponent (GHE), Margaos-Sun FD (MSFD), Leibovich-Toth FD (LTFD), and fractal volatility (FV), to calculate the FD of a given time series.

Table 1.  The definition of fractal-based features.

Feature	Definition
HFD [22]	${\left\langle{\frac{N-1}{{k}^{2}\left[\frac{N-m}{k}\;\right]}\sum \limits_{i=1}^{\left[\frac{N-m}{k}\;\right]}\left|x\left(m+ik\right)-x\left(m+\left(i-1\right)k\right)\right|}\right\rangle}_{k}\propto {k}^{-D}$ $for\; m={\mathrm{1, 2}}, \dots k\; and\; k={\mathrm{1, 2}}, \dots, {k}_{max}\left({k}_{max}=30\right)$	(1)
KFD [23]	$D=\frac{{\mathit{log}}\left(N\right)}{{\mathit{log}}\left(N\right)+{\mathit{log}}\left(\frac{{\mathit{max}}\left(dist\left(1, j\right)\right)}{\sum _{i=1}^{N}dist\left(i-1, i\right)}\;\right)}$ $for\; j={\mathrm{1, 2}}, \dots, N$	(2)
GHE [24]	$\frac{\left\langle{{\left|x\left(t\right)-x\left(t-\tau \right)\right|}^{q}}\right\rangle}{\left\langle{{\left|x\left(t\right)\right|}^{q}}\right\rangle}\sim{\left(\frac{\tau }{\upsilon }\right)}^{qD}$ $for\; 1\le \tau \le {\tau }_{max}\; and\; q=1$	(3)
MSFD [25]	${\mathit{ln}}\left(\frac{\sum _{n=1}^{N}\left[\left(x\left(t\right){\oplus}_{s}{b}_{t}^{\oplusϵ}\right)-\left(x\left(t\right){\ominus }_{s}{b}_{t}^{\oplusϵ}\right)\right]}{{\left(\frac{2k}{N}\right)}^{2}}\right)\propto D{\mathit{ln}}\left(\frac{1}{\frac{2k}{N}}\right)$ $for\; S={\mathrm{0, 1}}, 2, \dots, N, 0\le t\le T, ϵ={\mathrm{1, 2}}, \dots, {ϵ}_{max}, \; and\; k={\mathrm{1, 2}}, \dots, {k}_{m}$	(4)
LTFD [26]	$D=\underset{ϵ\to 0}{{\mathit{lim}}}\frac{{\mathit{log}}\left({N}_{bc}\left(ϵ\right)\right)}{{\mathit{log}}\left(\frac{1}{ϵ}\right)}$	(5)
FV [27]	$D=\underset{ϵ\to 0}{{\mathit{lim}}}\frac{{\mathit{log}}\left({N}_{bc}\left(ϵ\right)\right)}{{\mathit{log}}\left(\frac{1}{ϵ}\right)}$	(6)

---

 | Show Table

DownLoad: CSV

In Table 1, $x\left(t\right)$ is a given time series, $N$ is the number of data samples, $m$ is the initial time, $k$ is the time interval, $T$ is the period of the time series $x\left(t\right)$, $ϵ$ is the considered scale, ${k}_{m}$ is set based on the rules in ^[25], $b$ is the structuring element, and ${N}_{bc}\left(ϵ\right)$ is the number of pieces that covers the graph of $x\left(t\right)$. Also, $\left\langle{\dots }\right\rangle$ denotes averaging over time, $\oplus$ and $\ominus$ respectively refer to the erosion and dilation operations. Note that the FV algorithm is a box-counting method, including random-walk processes.

Similarly, Table 2 contains the definitions of six well-known methods of measuring the entropy of a time series, including Approximate Entropy (ApEn), Sample Entropy (SampEn), Shannon Entropy (ShanEn), Permutation Entropy (PermEn), Fuzzy Entropy (FuzEn), and corrected Conditional Entropy (CondEn).

Table 2.  The definition of entropy-based features.

Feature	Definition
ApEn [28]	$E=\underset{N\to \infty }{{\mathit{lim}}}\left(\frac{\sum _{i=1}^{N-m+1}\; {\mathit{log}}\; {C}_{i}^{m}\left(r\right)}{N-m+1}-\frac{\sum _{i=1}^{N-m+1}\; {\mathit{log}}\; {C}_{i}^{m+1}\; \left(r\right)}{N-m+1}\right)$	(7)
SampEn [29]	$E=-{\mathrm{ln}}\frac{\sum _{i=1}^{N-m}\; {\mathit{log}}\; {A}_{i}^{m}\left(r\right)}{\sum _{i=1}^{N-m}\; {\mathit{log}}\; {B}_{i}^{m}\left(r\right)}$	(8)
ShanEn [30]	$E=-\sum \limits_{L}{p}_{L}\; {\mathit{log}}\; {p}_{L}$	(9)
PermEn [31]	$E=-\sum\limits _{L}{\pi }_{L}\; {\mathit{log}}\; {\pi }_{L}$	(10)
FuzEn [32]	$E={\mathrm{ln}}\left(\frac{1}{N-m}\sum \limits_{i=1}^{N-m}\left(\frac{{\sum }_{j=1, j\ne i}^{N-m}\; {D}_{ij}^{m}}{N-m-1}\right)\right)-{\mathrm{ln}}\left(\frac{1}{N-m}\sum \limits_{i=1}^{N-m}\left(\frac{{\sum }_{j=1, j\ne i}^{N-m}\; {D}_{ij}^{m+1}}{N-m-1}\right)\right)$	(11)
CondEn [30]	$E=\sum\limits _{L}{p}_{L}\; {\mathit{log}}\; {p}_{L}-\sum \limits_{L-1}{p}_{L-1}\; {\mathit{log}}\; {p}_{L-1}+{E}_{c}\left(L\right)$	(12)

---

 | Show Table

DownLoad: CSV

In Table 2, ${C}_{i}^{m}\left(r\right)$ is the number of blocks of successive values (with length $m$) that are similar to a certain block within the resolution $r$, ${A}_{i}^{m}\left(r\right)$ is the model vector, ${B}_{i}^{m}\left(r\right)$ is the template vectors, ${D}_{ij}^{m}$ is the similarity degree of two time series ${x}_{i}\left(t\right)$ and ${x}_{j}\left(t\right)$ (of length $m$), ${p}_{L}$ is the joint probability, ${\pi }_{L}$ is the relative frequency, $L$ is the dimension of the phase space, and ${E}_{c}\left(L\right)$ is the corrective term.

3.2.   Classical features

Transform-based feature extraction is an important technique in signal processing. Therefore, such transform-based features, as well as the statistical ones, were considered as the two most used groups of classical features ^[33,34]. Among the famous transform-based features, Discrete Wavelet Transform (DWT), Discrete Fourier Transform (DFT), Short-Time Fourier Transform (STFT), Discrete Cosine Transform (DCT), Hilbert Transform (HT), and Stockwell Transform (ST) are formulated in Table 3.

Table 3.  The definition of Transform-based features.

Feature	Definition
DWT [35]	$W\left(k, l\right)=\frac{1}{s}\sum\limits _{n}\sum \limits_{m}x\left(m.n\right)\psi \left(\frac{m-k}{s}\right)\psi \left(\frac{n-l}{s}\right)$	(13)
DFT [36]	$F\left(k\right)=\sum \limits_{n=0}^{N-1}x\left(n\right){e}^{-\frac{j2\pi }{N}kn}\; for\; k=0, \dots, N-1$	(14)
STFT [37]	$SF\left(m, f\right)=\sum \limits_{n=0}^{N-1}x\left(n\right)g\left(n-{w}_{0}\right){e}^{-j2\pi fn}$	(15)
DCT [38]	$C\left(k\right)=\sqrt{\frac{2}{N}}\sum \limits_{n=1}^{N}\frac{x\left(n\right)}{\sqrt{1+{\delta }_{k1}}}{\mathrm{cos}}\left(\frac{\pi }{2N}\left(2n-1\right)\left(k-1\right)\right)$	(16)
HT [39]	$H\left(t\right)=\frac{1}{\pi }{\int }_{-\infty }^{+\infty }\frac{x\left(\tau \right)}{t-\tau }d\tau$	(17)
ST [40]	$S\left(i{T}_{s}, \frac{n}{N{T}_{s}}\right)=\sum \limits_{m=0}^{N-1}F\left(\frac{m+n}{N{T}_{s}}\right){e}^{-\frac{2{\pi }^{2}{m}^{2}}{{n}^{2}}}{e}^{\frac{j2\pi mi}{N}}$	(18)

---

 | Show Table

DownLoad: CSV

In Table 3, $\psi$ is the mother Wavelet function, $k$ and $l$ are the positions of the wavelet $\psi$ with the scale of $s$, $g\left(M\right)$ is the window function of length $M$ and the center of ${w}_{0}$, $f$ is the frequency, ${\delta }_{k1}$ is Kronecker delta, and ${T}_{s}$ is the time sampling interval.

Table 4 describes six statistical measures as the classical and primary features of a given time series, i.e., mean, variance, skewness, kurtosis, median, and mode. In Table 4, ${f}_{i}$ is the frequency of data value, and $N$ is the total of data points.

Table 4.  The definition of statistical-based features.

Feature	Definition
Mean	${\mu }_{1}=\frac{1}{N}\sum\limits _{n=1}^{N}x\left(n\right)$	(19)
Variance	${\mu }_{2}=\frac{1}{N}\sum\limits _{n=1}^{N}{\left(x\left(n\right)-{\mu }_{1}\right)}^{2}$	(20)
Skewness	${\mu }_{3}=\frac{1}{N.{\mu }_{2}^{\frac{3}{2}}}\sum\limits _{n=1}^{N}{\left(x\left(n\right)-{\mu }_{1}\right)}^{3}$	(21)
Kurtosis	${\mu }_{4}=\frac{1}{N.{\mu }_{2}^{2}}\sum\limits _{n=1}^{N}{\left(x\left(n\right)-{\mu }_{1}\right)}^{4}.$	(22)
Median	$med=\left\{\begin{array}{l}x\left(\frac{n}{2}\right)\; \; \; \; \; for\; odd\; n\\ \frac{1}{2}\left(x\left(\frac{n-1}{2}\right)+x\left(\frac{n+1}{2}\right)\right)\; \; \; for\; even\; n\end{array}\right.$	(23)
Mode	$mod={\mathrm{max}}\left({f}_{i}\right)$ $for\; i=1, \dots, N$	(24)

---

 | Show Table

DownLoad: CSV

4.   Results and discussion

Before applying the methods, briefly described in Tables 1–4, the spiking activity of each MT neuron was averaged across trials for IN and OUT conditions. Then the averaged spiking activity in IN conditions as well as the OUT condition were used for feature extraction. Note that the first 400 ms of the memory period were excluded from any further analysis to prevent any data contamination caused by the disappearance of the visual target (see Figure 2a). In the next step, a total of 24 features, including fractal- and entropy-based features as the nonlinear features and transform- and statistical-based features, were extracted from the IN and OUT conditions in both fixation and memory periods by applying the algorithms summarized in Tables 1–4. An indicator of the memory content is expected to make a significant difference between the IN and OUT conditions during the memory period while showing no difference between the conditions during the fixation period. Therefore, using the Wilcoxon signed-rank test, the features of IN conditions were compared to the corresponding values of the OUT condition in both memory and fixation periods. Tables 5–8 express the results (p-values) of the examination for each class of nonlinear and classical features.

Table 5.  The comparison of the fractal-based features to detect the memory content from the average spiking activity of MT neurons in fixation and memory periods (between IN and OUT conditions).

Index	p-value
	Fixation Period	Memory Period
HFD	0.4637	< 0.001
KFD	$< 0.001$	$< 0.001$
GHE	$< 0.001$	$< 0.001$
MSFD	0.1295	< 0.001
LTFD	$< 0.001$	$< 0.001$
FV	$< 0.001$	$< 0.001$

---

 | Show Table

DownLoad: CSV

Table 6.  The comparison of the entropy-based features to detect the memory content from the average spiking activity of MT neurons in fixation and memory periods (between IN and OUT conditions).

Index	p-value
	Fixation Period	Memory Period
ApEn	$< 0.001$	$0.1827$
SampEn	$< 0.001$	$0.0114$
ShanEn	0.5488	< 0.001
PermEn	$< 0.001$	$< 0.001$
FuzEn	$< 0.001$	$0.0022$
CondEn	0.2927	< 0.001

---

 | Show Table

DownLoad: CSV

Table 7.  The comparison of the transform-based features to detect the memory content from the average spiking activity of MT neurons in fixation and memory periods (between IN and OUT conditions).

Index	p-value
	Fixation Period	Memory Period
DWT	$0.0449$	$0.6004$
DFT	$< 0.001$	$< 0.001$
STFT	$< 0.001$	$< 0.001$
DCT	$< 0.001$	$< 0.001$
HT	$0.0069$	$0.0006$
ST	$< 0.001$	$< 0.001$

---

 | Show Table

DownLoad: CSV

Table 8.  The comparison of the statistical-based features to detect the memory content from the average spiking activity of MT neurons in fixation and memory periods (between IN and OUT conditions).

Index	p-value
	Fixation Period	Memory Period
Mean	$0.0301$	$0.6165$
Variance	$< 0.001$	$< 0.001$
Skewness	0.0706	< 0.001
Kurtosis	$0.2252$	$0.2377$
Median	$0.0007$	$0.2856$
Mode	$< 0.001$	$0.1814$

---

 | Show Table

DownLoad: CSV

Looking more closely at Table 5, it can be noticed that the HFD (Figure 3a, b) and MSFD (Figure 4a, b) may show a difference in memory period (${p}_{HFD}\; and\; {p}_{MSFD} < 0.001$) while no difference was detected in the fixation period (${p}_{HFD} = 0.4637$ and ${p}_{MSFD} = 0.1295$). Based on Table 6, the same effect was also observed for the ShanEn (Figure 5a, b; fixation period: ${p}_{ShanEn} = 0.5488$; memory period: ${p}_{ShanEn} < 0.001$) and CondEn (Figure 6a, b; fixation period: ${p}_{CondEn} = 0.2974$; memory period: ${p}_{CondEn} < 0.001$). Although four nonlinear features showed an effect (in memory period) compared to their baseline (fixation period), none of the transform-based features depicted such an effect (see Table 7). On the other hand, one of the statistical-based features, skewness in Table 8, showed a slight difference between IN and OUT conditions in the memory period (${p}_{skewness} = 0.0706$) compared to the fixation period (${p}_{skewness} < 0.001$).

Figure 3.  First row: The HFD values of IN versus OUT conditions during a) fixation and b) memory periods. Second row: The HFD values during fixation versus memory periods of c) IN and d) OUT conditions. Three asterisks show $p < 0.001$ and $n.s.$ indicates $p > 0.05$.

DownLoad: Full-Size Img PowerPoint

Figure 4.  First row: The MSFD values of IN versus OUT conditions during a) fixation and b) memory periods. Second row: The MSFD values during fixation versus memory periods of c) IN and d) OUT conditions. Three asterisks show $p < 0.001$ and $n.s.$ indicates $p > 0.05$.

DownLoad: Full-Size Img PowerPoint

Figure 5.  First row: The ShanEn values of IN versus OUT conditions during a) fixation and b) memory periods. Second row: The ShanEn values during fixation versus memory periods of c) IN and d) OUT conditions. Three asterisks show $p < 0.001$ and $n.s.$ indicates $p > 0.05$.

DownLoad: Full-Size Img PowerPoint

Figure 6.  First row: The CondEn values of IN versus OUT conditions during a) fixation and b) memory periods. Second row: The CondEn values during fixation versus memory periods of c) IN and d) OUT conditions. Three asterisks show $p < 0.001$ and $n.s.$ indicates $p > 0.05$.

DownLoad: Full-Size Img PowerPoint

The above analysis suggests that the HFD, MSFD, ShanEn, CondEn, and skewness could detect the presence of the working memory since they could perfectly distinguish IN condition from the OUT condition in the memory period, while no difference was detected in the fixation period. However, the desired indicator of the working memory should be sensitive to the presence of memory modulation as well. Therefore, we performed further analysis to identify whether these features can detect the presence of working memory or not. In other words, the same effect mentioned above should be noticed between the fixation (when no memory is present) and memory periods of IN condition. The comparison of HFD values of memory IN versus fixation IN and memory OUT versus fixation OUT (Figure 3c, d) revealed that HFD could capture the presence of working memory since it made a difference between memory IN and fixation IN (${p}_{HFD} < 0.001$), while it showed no difference between memory OUT and fixation OUT (${p}_{HFD} = 0.4735$). Furthermore, it can be seen that the value of HFD increases in the presence of working memory ($\Delta HF{D}_{memory\; IN-memory\; OUT} = 0.0708$ and $\Delta HF{D}_{memory\; IN-fixation\; IN} = 0.0710$), indicating that the spiking activity of MT neurons becomes more complex as the monkey maintains the spatial information. The same analysis was performed on MSFD (Figure 4c, d), ShanEn (Figure 5c, d), CondEn (Figure 6c, d), and skewness (Figure 7c, d). However, a significant difference was observed in the comparison of memory IN versus fixation IN ($p < 0.001$) and memory OUT versus fixation OUT ($p < 0.001$). Thus, it can be concluded that the MSFD, ShanEn, CondEn, and skewness perhaps are also sensitive to other cognitive factors such as vigilance, arousal and awareness as well as working memory. In any case, they cannot be considered as indicators unique to working memory.

Figure 7.  First row: The Skewness values of IN versus OUT conditions during a) fixation and b) memory periods. Second row: The Skewness values during fixation versus memory periods of c) IN and d) OUT conditions. Three asterisks show $p < 0.001$ and $n.s.$ indicates $p > 0.05$.

DownLoad: Full-Size Img PowerPoint

5.   Conclusions

In this paper, a comprehensive analysis was performed to examine the ability of the nonlinear features, including fractal- and entropy-based features, and classical ones, such as transform- and statistical-based features. Six fractal-based (HFD, KFDm GHE, MSFD, LTFD, and Fv), six entropy-based (ApEn, SampEn, ShanEn, PermEn, FuzEn, and CondEn), six transform-based (the average value of DWT, DFT, STFT, DCT, HT, and ST) and six statistical-based (mean, variance, skewness, kurtosis, median, and mode) features were extracted from the spiking activity of MT neurons. Using the Wilcoxon signed-rank test, values of each feature in IN and OUT conditions were compared during the memory and fixation period. As a result, HFD and MSFD from the fractal-based features, ShanEn and CondEn from the entropy-based features, and skewness from the statistical-based features were found to be capable of capturing the memory content. However, more investigations revealed that among the five mentioned features, HFD was the only indicator that was unique to the working memory since it shows a significate difference in the comparison of memory IN versus memory OUT and memory IN versus fixation IN while showing no difference between fixation IN and fixation OUT as well as memory OUT and fixation OUT. Such effect was not observed for any of the other four mentioned features. In addition, our results were consistent with the ones declared in ^[19]. Although the rest of mentioned features, in fact, captured the working memory, they could also detect other cognitive factors (such as vigilance and arousal), which make them not unique to working memory. As four out of five features were nonlinear, it can be concluded that the working memory alters the nonlinearity of visual responses in area MT that cannot be detected using simple statistical measures or powerful transform-based features.

Since the brain is a complex structure, it exhibits complex behavior reflected in biological signals such as neural activities. Thus, it is not far-fetched that this complexity can be captured by the nonlinear features measuring the complexity of nonlinear and nonstationary signals, such as the fractal dimension. Previously, we showed that although no signs of spatial working memory could be observed directly in the average spiking activity of neurons in the area MT, the nonlinear properties of spiking patterns were significantly changed under the influence of working memory captured by HFD. In this study, we examined the hypothesis that the effect of working memory can be noticed in the other feature domains, such as time, frequency, time-frequency, and entropy. However, the results revealed that linear features are not qualified to capture the content of working memory in the activity of individual neurons in the area MT. On the other hand, we found other criteria (except for the HFD) sensitive to the visual information in the memory period, whilst they were not unique to the working memory. Moreover, our results confirm that maintaining visual information increases the complexity or the dimensionality of the neural spiking activity. Overall, our results shed more light on the underlying neural mechanisms of spatial working memory and may yield insights regarding the neural substrates of other cognitive functions, such as attention in future studies.

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.