## Abstract

With the development of neural recording technology, it has become possible to collect activities from hundreds or even thousands of neurons simultaneously. Visualization of neural population dynamics can help neuroscientists analyze large-scale neural activities efficiently. In this letter, Laplacian eigenmaps is applied to this task for the first time, and the experimental results show that the proposed method significantly outperforms the commonly used methods. This finding was confirmed by the systematic evaluation using nonhuman primate data, which contained the complex dynamics well suited for testing. According to our results, Laplacian eigenmaps is better than the other methods in various ways and can clearly visualize interesting biological phenomena related to neural dynamics.

## 1  Introduction

The remarkable computational capability of the brain resides in its complex neural networks, which tend to encode a piece of information through the coordinated activities of a large population of neurons. For instance, a single motor or cognitive behavior may involve hundreds or thousands of neurons (Sanes & Donoghue, 2000; Churchland, Cunningham, Kaufman, Ryu, & Shenoy, 2010; Rigotti et al., 2013). The need to study neuron populations has led to the rapid development of large-scale recording technologies (Alivisatos et al., 2013; Cunningham & Yu, 2014; Breakspear, 2017), as well as numerous studies about the simultaneous recording of a large number of neurons (Kipke et al., 2010; Spira & Hai, 2013; Viventi et al., 2011; Kerr & Denk, 2008; Mank et al., 2008; Cui et al., 2013; Brunton, Johnson, Ojemann, & Kutz, 2016; Kozai, Eles, Vazquez, & Cui, 2016; Chen, Canales, & Anikeeva, 2017). However, the constantly increasing size and dimensionality of neural activity data are posing new challenges to data analysis. The major challenge is how to convert such high-dimensional dynamic data into an intuitive form suitable for human interpretation. To achieve this, the data dimensionality needs to be reduced to the level on which the entire population dynamics in an experiment can be visualized in a single image, with the constraint of preserving as much useful information as possible to produce informative visualization.

Linear dimensionality reduction methods such as principal component analysis (PCA) and multidimensional scaling have been widely used in neuroscience research (Churchland, Yu, Sahani, & Shenoy, 2007; Churchland et al., 2010, 2012; Mante et al., 2013; Young & Yamane, 1992; Rolls & Tovee, 1995; Op de Beeck, Wagemans, & Vogels, 2001; Kayaert, Biederman, & Vogels, 2005; Kiani, Esteky, Mirpour, & Tanaka, 2007; Lehky & Sereno, 2007). Although these methods can visualize neural activity changes at different epochs of an experiment, they may fail to capture a nonlinear relationship among the spike signals (Richardson, 2008). Therefore, nonlinear dimensionality-reduction methods, such as isometric feature mapping (Isomap) and $t$-distributed stochastic neighbor embedding ($t$-SNE), have been applied to neural population analysis in recent studies (Adamos, Laskaris, Kosmidis, & Theophilidis, 2010, 2012; Song, Wang, & Viventi, 2017; Dimitriadis, Neto, & Kampff, 2018; Vargas-Irwin, Brandman, Zimmermann, Donoghue, & Black, 2015). However, applying these methods to large-scale data is challenging because of their high computational complexity. For instance, the $t$-SNE method needs to model the distribution of neighbors of each data point and minimize the difference between high-dimensional and low-dimensional distribution. Therefore, when there are many sampling points in the neural data, the $t$-SNE method takes a long time to perform the optimization steps, causing difficulties in converging to an optimal solution with limited computational resources. Therefore, in large-scale neural population data analysis, it is essential to find a method that can respond to external condition changes fast and effectively.

In this letter, we propose to use the Laplacian eigenmaps (LapEig) method to reduce the dimension of a spike count matrix for visualization. The method represents an efficient computational approach for nonlinear dimensionality reduction that has locality-preserving properties and a natural connection to clustering. In the proposed method, the essential structural features of neural population data are preserved. The proposed method is tested and validated experimentally using neuron population activity data obtained from the dorsal premotor area (PMd) of monkeys. The experimental results show that neural activity patterns can be determined based on the monkeys' behavior. Moreover, different epochs of the movement can be visually identified. Compared to the commonly used linear dimensionality-reduction method (the PCA method) and nonlinear dimensionality-reduction methods (the Isomap and $t$-SNE methods), we found that LapEig achieved higher classification accuracy and lower within-class distance. Nevertheless, LapEig had much lower computational complexity than the two nonlinear methods, so it could be used as an effective tool for exploring dynamic processes in large-scale neural population data.

## 2  Laplacian Eigenmaps

LapEig represents a graph-based dimensionality-reduction method, proposed by Belkin and Niyogi (2003) that is formulated in order to preserve the distance graph of the input data points. In other words, the LapEig method tries to keep neighboring data points in the same space as the neighbors in a low-dimensional space. Given a data set of $m$ points ${x1,…,xm⊂Rl$, the LapEig method first computes an adjacency map by treating each data point as a graph node and creating edges from each node to its $k$ nearest neighbors. The edge weight $Wi,j$ between two nodes $(xi,xj)$ is defined by
$Wi,j=e-∥xi-xj∥22σiandjareconnected,0otherwise,$
(2.1)
where $σ$ denotes the scale factor. The weight function ensures that points that are closer to each other have larger edge weights. Therefore, finding a $d$-dimensional ($d representation $yi$ of $xi$ by minimizing $∑i=1m∑j=1m∥yi-yj∥2Wij$ leads to a solution that preserves neighboring relationships. Using the weight matrix $W$ and the degree matrix $D$, a diagonal matrix $Dii=∑j=1mWij$ can be obtained and the optimization problem can be further deduced to the following form:
$minYtr(YTLY)s.t.YTDY=I,$
(2.2)
where $L$ is called the Laplacian matrix, and $L=D-W$, and the $i$th row of $Y$ is $yiT$. The constraint defines the scale of the solution and eliminates trivial solutions such as all zero points. The deduced form can be efficiently solved by computing the first $d$ eigenvectors of $L$ with a time complexity of $O((d+k)m2)$ because $L$ is usually sparse ($k≪m)$ (Suzuki, 2015).

The key step in applying the LapEig method to neural activity data analysis is to measure the distance between two neural activity patterns. The raw form of neural activity data obtained through electrophysiological recording denotes the time series of voltage measurement, which can be sorted into spikes of one or more neurons. After the spike sorting, the activity of each neuron is represented by a vector consisting of binary values in which 0 or 1, respectively, means there is one or no spike at the corresponding time point. The spike count array can be further compressed by counting the spikes at $n$ evenly distributed time points $(t1,t2,…,tm)$, resulting in a new vector $Si=(Si1,Si2,…,Sim)T$, where $Sij$ denotes the number of spikes of the $i$th neuron within the time window [$tj-tbintj]$. Given $l$ neurons, the LapEig method takes the matrix $S=(S1,S2,…,Sl)T$ as its input for dimensionality reduction by treating each column vector as a data point.

## 3  LapEig Method Validation Using Primate Cortical Ensemble Activity Data

We used the neural population activity data of the PMd of monkeys to validate the LapEig method. The neural signals were collected by a 96-channel Blackrock microelectrode array. The care of the monkeys and the experimental protocols complied with the Guide for the Care and Use of Laboratory Animals (China Ministry of Health) and were approved by the Animal Care Committee at Zhejiang University, China. To test our method, we conducted two types of experiments using different data sets. In the first experiment, the single-object task (SO task), the monkeys were trained to grasp a single object by using two types of gestures, and in the second experiment, the tuning table task (TT task), monkeys were trained to grasp four different objects by using four different types of gestures. The experimental data and actions of the monkeys corresponding to different epochs were collected by the sensors and then recorded.

In this experiment, two macaque monkeys were trained to perform a delayed grasping task. The experimental flowchart and gesture types are shown in Figure 1. Each monkey was trained to grasp a handle on a vertical plexiglass board with its left hand while sitting in front of the board. The monkey needed to grasp the handle with either a power or hook gesture depending on the instructions given. A trial started after the monkey placed the acting hand on the resting button below the handle. Once the button was pressed, the lower LED turned red (fix signal), indicating the start of the 700 ms fix epoch, which was followed by lighting up one of the two upper LEDs which instructed the monkey how to grasp the handle in the trial (gesture cue). The red LED indicated the power gesture, and the green LED indicated the hook gesture. After seeing the gesture cue, the monkey had to keep the button pressed for a random period ranging from 600 ms to 1000 ms, which denoted a delayed plan epoch, until the go signal was turned on, instructing the monkey to release the button, reach out its hand, and grasp the handle. Monkeys were requested to hold the appropriated grip for 300 ms. The monkeys received a liquid reward when they had performed a correct grip. All trials were interleaved randomly. For the purpose of a subsequent analysis, we subdivide a trial into six epochs: (1) fix: From when fix signal was turned on to gesture cue was turned on (700 ms); (2) cue: from when gesture cue was turned on to gesture cue was turned off (800 ms); (3) plan: from when gesture cue was turned off to go signal was turned on (600–1000 ms); (4) move: from when go signal was turned on to the beginning of object touching; (5) hold: from object touching to go signal was turned off (500 ms); and (6) reward: after go signal was turned off.

Figure 1:

The experimental flowchart and gesture types of the single-object task. (Top) Gesture cue signals that instruct the power and hook gestures. (Middle) Flowchart and six different epochs of a trial. (Bottom) The different states (on and off) of the buttons and LEDs.

Figure 1:

The experimental flowchart and gesture types of the single-object task. (Top) Gesture cue signals that instruct the power and hook gestures. (Middle) Flowchart and six different epochs of a trial. (Bottom) The different states (on and off) of the buttons and LEDs.

Similar to the single-object task, two macaque monkeys were trained to perform a delayed turning table task. In this task, monkeys were seated in front of the turning table and were trained to grasp one of four different objects—handle, flatbed, ring, and small column—corresponding to the four actions of grasping, gripping, hooking, and pinching, respectively. A trial started after the monkeys pressed the button and the two-color LED light turned to red. Monkeys needed to keep pressing the button for 700 ms, which is represented by the fix epoch in Figure 2. After the fix epoch, the light turned on, and the target object was shown to the monkey. After 800 ms, the light turned off, and the monkeys had to keep holding the button for the next 400 to 800 ms. After the LED light turned green, the monkeys had to release the button and start to grasp the object with the corresponding gesture. When the performed grip was correct and the object was held for 500 ms, the trial was considered successful. Then the LED light turned off, and the monkeys received water as a reward. After that, the turning table randomly chose the next object and waited for the next trial to begin.

Figure 2:

The experimental flowchart and gesture types of the turning table task. (Top) The four objects and the corresponding gestures (grasp, grip, hook, and pinch). (Middle) Flowchart and the six epochs of a trial. (Bottom) Different states (on and off) of the buttons and LEDs.

Figure 2:

The experimental flowchart and gesture types of the turning table task. (Top) The four objects and the corresponding gestures (grasp, grip, hook, and pinch). (Middle) Flowchart and the six epochs of a trial. (Bottom) Different states (on and off) of the buttons and LEDs.

The interval between every two trials was 1.5 seconds. As shown in Figure 2, we subdivided a trial into six epochs: (1) fix: From the beginning of the experiment to the moment the light was on; (2) light: from the moment the light was on to the moment the light was off; (3) plan: from the moment the light was off to the moment the go signal was turned on; (4) move: from the moment the go-signal was turned on to object touching (the moment when the object was touched); (5) hold: from the moment the object was touched to the moment the go signal turned off; and (6) reward: after the moment the monkey was given the water reward.

### 3.2  Electrophysiological Recording

After the training (the performance was higher than 85% for training lasting five to seven months), a 96-channel microelectrode array, 4.0 mm $×$ 4.0 mm (Blackrock Microsystems, USA), was implanted in the right PMd area. The neural data were collected by the neural processing system (Blackrock Microsystems, Salt Lake City, UT) with the sampling frequency of 30 kHz and the bandpass filters with the bandwidth of 250 Hz to 7.5 kHz. We detected the spike signals by using the threshold detection method where the threshold was set to $-$5.5 times the root mean square (RMS) of the baseline signal level. The single-unit activities were further isolated by the Offline Sorter (Plexon, Dallas, TX).

The data collection process for each of the two tasks lasted 15 days. To facilitate the evaluation of the LapEig method performance, we extracted one-day data from each of the two data sets to be used in the following analysis. After sorting, 215 units were isolated in the single-object task. During the experiment, which lasted 26 minutes, the monkeys achieved a success rate of 94.2%. In the turning table task, 240 units were isolated; the experiment lasted 33 minutes, and the monkeys achieved a success rate of 97.1%. Four example units for each task are shown in Figure 3, where it can be seen that different units had different firing patterns for different gestures. Then we used the LapEig method to perform dimensionality-reduction analysis on these two data sets to test the effect of the classification and visualization of the proposed method.

Figure 3:

Raster plot of four example units of each task. Spikes are aligned to the start time of the go signal (at 0 s). The color points in the raster displays denote the spike firing time in trials corresponding to different gestures. The black vertical ticks in the raster displays denote the experimental events. (a) Example units of the single-object task. (b) Example units of the turning table task.

Figure 3:

Raster plot of four example units of each task. Spikes are aligned to the start time of the go signal (at 0 s). The color points in the raster displays denote the spike firing time in trials corresponding to different gestures. The black vertical ticks in the raster displays denote the experimental events. (a) Example units of the single-object task. (b) Example units of the turning table task.

### 3.3  Dimensionality Reduction

For both tasks, spikes were counted at every 100 ms starting from a chosen time point using the time window with the length $tbin=500ms$. The overlap between two adjacent time windows introduced a smoothing effect to stabilize the distance measurement. The two free parameters $k$ and $σ$ of the LapEig method were generally determined empirically. In this study, we used $k=12$ and $σ=1$, based on spike sorting studies that show that setting $k$ to a value between 5 and 21 and $σ$ to be bigger than 0.6 could yield good results (Ghanbari, Papamichalis, & Spence, 2010; Chah et al., 2011). After being reduced to three dimensions by the LapEig method, the neuron dynamics data were plotted on a computer screen using 3D data visualization.

In the single-object task, the average time of a trial was about 5 seconds, and we selected the data of 80 continuous trials for subsequent analysis and processing (with an overall duration of 400 seconds). The monkey made the correct action in all the selected trials. In the selected data set, there were 40 trials for each of the two gestures. Overall, there were 215 neuron samples, each represented as a vector of 4000 spike counts (four initial vectors were calculated using data from $-$400 ms to 0 ms).

Figure 4a shows a line chart of 4000 time points after dimensionality reduction, connected in chronological order. We found that the dimensionality reduction resulted in two rings in 3D space that partially overlapped. The connection of the lines represented the change in neural population activity with time. Figure 4a illustrates that most of the trials changed in one direction without oscillating back and forth. These results suggested that the neural population activity gradually deviated from the baseline as the experiment was progressing and then returned to the baseline when the experiment was completed. In order to demonstrate the relationship between movement epoch and population activity, we used different colors to denote different experimental epochs. At the beginning of the experiment, the population was in the baseline state. After the experiment started, the monkeys concentrated on preparing for the experiment before gesture cue was turned on and the firing patterns of the neural population began to deviate from the baseline level, but the monkeys did not know which gesture to make, so the neural population activity did not show any selectivity. When gesture cue was turned on, the monkeys started to prepare themselves for the next epoch. The firing of the corresponding neural populations began to separate different gestures. After the monkeys made the corresponding gestures, more motor-related neurons showed selectivity, and the differences between the neural population activities became more obvious. Then the neural population activity began to decrease, and it returned to the baseline level when the trials came to an end and the monkeys were rewarded.

Figure 4:

Dimensionality-reduction result of the LapEig method. (a) The result of the SO task. The trials related to the hook gesture are represented by black, blue, and green lines, and the trials related to the power gesture are presented by magenta, red, and yellow lines. Each color represents a different epoch of a trial; black and magenta illustrate the period from the trial start to gesture cue was turned on; blue and red illustrate the period from when gesture cue was turned on to go signal was turned on; and green and yellow illustrate the period from when go signal was turned on to the trial end. The circle denotes the time gesture cue was turned on in each trial, and the asterisk denotes the time go signal was turned on in each trial. (b) The result of the TT task. The trials related to grip, grasp, hook, and pinch are shown in blue, black, green, and red, respectively. The dot denotes the time light was turned on in each trial, and the asterisk denotes the time go signal was turned on in each trial.

Figure 4:

Dimensionality-reduction result of the LapEig method. (a) The result of the SO task. The trials related to the hook gesture are represented by black, blue, and green lines, and the trials related to the power gesture are presented by magenta, red, and yellow lines. Each color represents a different epoch of a trial; black and magenta illustrate the period from the trial start to gesture cue was turned on; blue and red illustrate the period from when gesture cue was turned on to go signal was turned on; and green and yellow illustrate the period from when go signal was turned on to the trial end. The circle denotes the time gesture cue was turned on in each trial, and the asterisk denotes the time go signal was turned on in each trial. (b) The result of the TT task. The trials related to grip, grasp, hook, and pinch are shown in blue, black, green, and red, respectively. The dot denotes the time light was turned on in each trial, and the asterisk denotes the time go signal was turned on in each trial.

Based on the experimental results, the LapEig method could clearly distinguish different motion and action epochs. However, only two types of gestures were used in the experiment, which might not be sufficient to test the LapEig method performance well. Therefore, to further verify the LapEig method, we performed an additional experiment involving more gestures. In addition, since there were more gestures in the turning table task than in the single-object task, we took the data of 120 continuous trials for analysis and verification (with the overall duration of 600 seconds). The four gestures were respectively prompted 30 times by a pseudo-random sequence, which denoted 120 trials in total. In all 120 trials, the monkeys made the correct gestures. After sorting, we had 240 isolated units. The generated spike count matrix $S$ had the dimensionality of $240×6000$ (four initial vectors were calculated using data from $-$400 ms to 0 ms). Finally, matrix $S$ was imported into the LapEig algorithm for dimensionality reduction. The visualization results are shown in Figure 4b.

Figure 4b shows 6000 time points of matrix $S$ after dimensionality reduction. Similar to the previous experiment, the dimensionality-reduction results of this experiment were clearly clustered into four rings, corresponding to the four gestures performed in the experiment. At the same time, when the population activity started to show a separation characteristic, the neural population firing was first separated into two clusters, and then, after a period of time, the neural population firing was separated into four clusters. In addition, in the most obvious clustering stage of the population activity, three clusters were close to each other, while one was farther away.

According to the experimental design, the following conclusion can be drawn. Before light was turned on, monkeys did not know which gesture to make, so the neural population activity did not show any selectivity. After light was turned on, the monkeys started to prepare themselves for the corresponding gestures, but some of the gestures were similar so the distinction was not obvious. When the go signal was turned on, the motor-related neurons began to participate in the process of gesture execution, and the distinctions among different gestures became more and more obvious. When the gesture was finally completed, the distinction between the clusters became obvious. At that point, the within-class distance of similar gestures was relatively small, and the between-class distance of highly different gestures was relatively large. The change of the grip gesture was particularly obvious, and the between-class distance of the other three gestures experienced a small increase.

### 3.4  Quantitative Analysis of Dimensionality-Reduction Results

In order to quantify the clustering effect of different epochs in a trial, we used a $k$-nearest-neighbors ($k$NN) classifier implemented by four-fold cross-validation to evaluate the accuracy of clustering under the two paradigms. The clustering results are shown in Figure 5. Figure 5a shows the clustering performance of the single-object task and Figure 5b the clustering performance of the turning table task. As shown in Figure 5a, the clustering performance was in accordance with that in Figure 3. In the fix epoch, the clustering accuracy was near chance level. In the cue epoch, clustering accuracy started to rise and reached about 95% in the plan epoch. In the move epoch, the clustering accuracy was the highest, reaching almost 100% until the end of the hold epoch. After the monkeys got the reward, the clustering accuracy dropped rapidly to chance level at the end of the trial. The results in Figure 5b are similar to those in Figure 5a. In the fix epoch, clustering accuracy was near chance level. In the light epoch, clustering accuracy started to rise and reached nearly 65% in the plan epoch. In the move epoch, clustering accuracy continued to increase and reached a maximum of about 95% in the hold epoch. After the monkeys got the reward, clustering accuracy began to decline and slowly returned to chance level at the end of the experiment.

Figure 5:

Clustering performance of the two tasks. The $k$-nearest-neighbors (kNN) classifier was implemented by leave-one-out cross-validation to evaluate clustering performance. The plots were aligned to the go signal as indicated by the vertical lines at 0 s. The black vertical lines separate different epochs of a trial. (a) Results of the SO task. (b) Results of the TT task.

Figure 5:

Clustering performance of the two tasks. The $k$-nearest-neighbors (kNN) classifier was implemented by leave-one-out cross-validation to evaluate clustering performance. The plots were aligned to the go signal as indicated by the vertical lines at 0 s. The black vertical lines separate different epochs of a trial. (a) Results of the SO task. (b) Results of the TT task.

The results shown in Figure 5 are in accordance with those in Figure 3. Since there were only two gestures in the single-object task, when the cue was given in the cue epoch, different types of activities were clearly separated in different clusters until the end of the trial (reward epoch). In the turning table task, the activities of the neural population did not show selectivity in the fix epoch before the light was turned on but the monkeys began to prepare themselves for different gestures in the light and plan epochs after the light was turned on. However, there was no movement, so the classification effect was not obvious. In the move epoch after the go signal was turned on, the motor-related neurons started to participate (Raos, Umiltá, Gallese, & Fogassi, 2004) and the distinctions between different gestures increased accordingly. In the hold epoch, the gestures were fully formed, and the distinctions of dimensional reduction results were the most obvious. After the trial was successfully completed, the neural population returned to the baseline level. The results in Figure 5 prove that the LapEig method could track the changes over time.

The other way to measure the effectiveness of the LapEig method in visualization is to compare the within-class and between-class distances. Generally a shorter within-class distance and a longer between-class distance result in better visualization performance. However, using clustering accuracy is not enough to measure the visualization quality because it can show only how much two clusters are mixed, not how far they are separated. Figure 6a shows the normalized within-class distance and between-class distance of the single-object task, and Figure 6b shows the normalized within-class distance and between-class distance of the turning table task. Figure 6 shows that the within-class distance and fluctuation of the LapEig method were very small in both tasks, which means that the neural population activity of the trials of the same gesture was very similar. The between-class distance was small in the fix and cue/light epochs because neurons did not show significant selectivity, but it increased in the plan epoch. When the motor-related neurons participated in the movement in the move epoch, the between-class distance increased rapidly and reached the peak value in the hold epoch when the gesture was fully formed. These results demonstrate that the LapEig method could not only effectively classify different grasping movements but also clearly show the changes in the distance between different actions.

Figure 6:

The normalized within-class and between-class distances of the trials. The dotted lines indicate the normalized within-class distances, and the solid lines indicate the normalized between-class distances. The plots are aligned to the go signal as indicated by the vertical lines at 0 s. The black vertical lines separate different epochs of trials. (a) The SO task. (b) The TT task.

Figure 6:

The normalized within-class and between-class distances of the trials. The dotted lines indicate the normalized within-class distances, and the solid lines indicate the normalized between-class distances. The plots are aligned to the go signal as indicated by the vertical lines at 0 s. The black vertical lines separate different epochs of trials. (a) The SO task. (b) The TT task.

### 3.5  Comparison with Other Methods

In order to demonstrate the advantage of the LapEig method in dynamics visualization applications, we ran similar experiments using the PCA method (a commonly used linear dimensionality-reduction method) and the Isomap and $t$-SNE methods (commonly used nonlinear methods).

The PCA method reduces data dimensions by trying to capture as much variation as possible from the data. Due to its simplicity, the PCA is one of the most commonly used methods in neuron dynamics visualization, but it can address only linear mapping, which often represents an oversimplification of neural dynamics data. The Isomap is a representative isometric mapping method, designed to preserve the geodesic distances defined on a manifold after mapping. In practice, it is often conducted by applying multidimensional scaling on a geodesic distance matrix consisting of pairwise path lengths of the data points on their neighborhood graph. Finally, the $t$-SNE also denotes a kind of distance-preserving method. However, instead of preserving distances directly, it tries to preserve the probability distribution of pairwise data distances. Its objective function is to minimize is the Kullback-Leibler divergence in the distance distribution before and after mapping.

In the experiment, the PCA, the Isomap, and $t$-SNE methods were used to reduce the dimension of a spike matrix $S$ in the single-object and turning table tasks. The dimensionality-reduction results are shown in Figure 7, where the left plots show the results of the single-object task and the right plots show the result of the turning table task. The dimensionality-reduction results of the PCA are shown in Figures 7a and 7b. In Figure 7a, it can be seen that in the single-object task, different epochs of the trial were roughly gathered together, but the two gestures were barely distinguishable. In the turning table task, the PCA distinguishing ability was worse than in the single-object task. Thus, the PCA failed to distinguish different motor behaviors in both tasks.

Figure 7:

Dimensionality-reduction results of the PCA, the Isomap, and the $t$-SNE methods. The left side shows the results of the SO task and the right side the results of the TT task. (a,b) Results of the PCA method. (c,d) Results of the Isomap method. (e,f) Results of the $t$-SNE method. The marks used in this figure are consistent with those used in Figure 3.

Figure 7:

Dimensionality-reduction results of the PCA, the Isomap, and the $t$-SNE methods. The left side shows the results of the SO task and the right side the results of the TT task. (a,b) Results of the PCA method. (c,d) Results of the Isomap method. (e,f) Results of the $t$-SNE method. The marks used in this figure are consistent with those used in Figure 3.

The results of the Isomap are shown in Figures 7c and 7d. In Figure 7c, the curves representing the dimensionality reduction for different gestures have noticeable separation, forming two rings. After the gesture cue was turned on, the neural activity of the two gestures started to show the difference and reached a peak after the go signal was turned on (yellow and green epochs). Then the activity gradually decreased to the baseline level, similar to the results presented in Figure 3a. However, in Figure 7c, the curves of different trials of the same gesture are scattered, and the within-class distance is large. Furthermore, as can be seen in Figure 7d, for the second task, the curves in different colors corresponding to the four gestures are obviously separated after the light-on and are gathered into four separate circular trajectories. Although the curves of the same gestures are discrete, the neural activity of different gestures still could be accurately distinguished. Figures 7e and 7f show the result of the $t$-SNE, where it can be seen that in the single-object task (see Figure 7e), the $t$-SNE distinguished two gestures after the gesture cue was turned on. Although the neural trajectory was scattered before the gesture cue was turned on, the curve of the gesture began to show a clear clustering after the gesture and the discrimination between different gestures improved, reaching a peak after the go signal. However, for the turning table task (see Figure 7f), the neural curves of different gestures were mixed together, so different gestures were indistinguishable.

Next, we quantitatively evaluated the clustering performance of the three algorithms. We calculated clustering accuracy, within-class distance, between-class distance, and their ratios for all three algorithms and compared results with that of the LapEig method. The comparison results are shown in Figure 8, where different methods are represented by different colors. The results of the single-object task and turning table task are presented on the left side and right side of Figure 8, respectively.

Figure 8:

Clustering accuracy, normalized within-class distance, and between-class distance of the four methods. The left side shows the results of the SO task, and the right side shows the results of the TT task. The results of the LapEig, $t$-SNE, PCA, and Isomap methods are shown in blue, red, magenta, and green, respectively. The dotted line denotes the within-class distance, and the solid line denotes the between-class distance. (a) Clustering accuracy. (b) Normalized within-class distance and between-class distance.

Figure 8:

Clustering accuracy, normalized within-class distance, and between-class distance of the four methods. The left side shows the results of the SO task, and the right side shows the results of the TT task. The results of the LapEig, $t$-SNE, PCA, and Isomap methods are shown in blue, red, magenta, and green, respectively. The dotted line denotes the within-class distance, and the solid line denotes the between-class distance. (a) Clustering accuracy. (b) Normalized within-class distance and between-class distance.

The clustering accuracy of the $k$NN classifier is presented in Figure 8a, where it can be seen that the accuracy of all the algorithms was close to chance level in the fix epoch, where the gesture information was not given. When the gesture information was given, the accuracy of the three nonlinear algorithms increased significantly and reached a maximum in the hold epoch after the gesture was performed. When task ended and the monkeys got the reward, the accuracy of the three nonlinear algorithms declined gradually to the chance level. The clustering quality of the PCA algorithm fluctuated around the channel level during the trial. The accuracy of the Isomap and LapEig methods in both tasks was basically the same. The accuracy of the $t$-SNE in the single-object task was similar to that of the LapEig and Isomap, but in the turning table task, its accuracy declined significantly and the peak value was about 65%, which was much lower than those of the LapEig and Isomap methods (about 95%). The clustering accuracy of the four algorithms basically agreed with the intuitive results presented in Figures 3 and 7. Therefore, the LapEig and Isomap methods achieved higher clustering accuracy.

The normalized within-class and between-class distances for the two tasks are shown in Figure 8b, where the dotted line indicates the within-class distance, and the solid line indicates the between-class distance. Figure 8b shows that the within-class distances of the LapEig method were much smaller than those of the other three methods. In fact, they were close to zero during the entire trial, leading to a thin loop pattern for each gesture in visualization, as shown in Figure 3. The within-class distances of the other three methods were much larger than that of the LapEig method. Although the within-class distance decreased after a period, it was still larger than that of the LapEig method after the gesture information was given. Analyzing the trend of the between-class distance, we found that the result of the PCA was fluctuating around zero. However, the results of the Isomap and LapEig increased slowly after the gesture cue was turned on. Then there was a sharp rise in the move epoch, and the peak was reached in the hold epoch in both tasks, so the difference between different gestures was maximized. The between-class distance of the $t$-SNE was always below 0.1 in the turning table task, and it was difficult to distinguish different gestures from the dimensional-reduction result. Besides, the LapEig method achieved the largest ratio of the between-class distance to the within-class distance, and it reached its maximum in the hold epoch and was much higher than those of the other methods. Consequently, the LapEig achieved the best clustering performance among all the tested methods.

We then evaluated the run-time efficiency of the four methods. The calculations were performed using Matlab R2012a on a Windows 7 pro 64 bit workstation with a 3.40 GHz Intel Core i5-3570 CPU and 8 GB RAM. The input data length of the four methods was from 100 s to 600 s. The computing time of the methods is shown in Figure 9. Among the three nonlinear dimensionality-reduction methods, the LapEig was several orders of magnitude faster than the other two methods. Due to the superior performance regarding clustering accuracy and visual discrimination, the LapEig method clearly outperformed all the tested methods in neural dynamics visualization, even though it was slower than the PCA method.

Figure 9:

Comparison of the running time of the four methods. The x-axis represents the input data length, and the $y$-axis represents the method running time in the exponential form. The results of the LapEig, $t$-SNE, PCA, and Isomap methods are shown in blue, red, magenta, and green, respectively.

Figure 9:

Comparison of the running time of the four methods. The x-axis represents the input data length, and the $y$-axis represents the method running time in the exponential form. The results of the LapEig, $t$-SNE, PCA, and Isomap methods are shown in blue, red, magenta, and green, respectively.

## 4  Discussion

With the development of neural recording technology, it has become possible to record signals from thousands of neurons simultaneously (Grewe, Langer, Kasper, Kampa, & Helmchen, 2010; Ahrens et al., 2012; Ahrens, Orger, Robson, Li, & Keller, 2013; Rios, Lubenov, Chi, Roukes, & Siapas, 2016). However, the extraction of effective behavioral information from the spike signals is still challenging. In this letter, we present the LapEig method, which highlights the similarity and dynamical process of neural population activity. The proposed method is based on the intrinsic properties of neural activity data and does not require any external information about motor or cognitive behaviors. Moreover, it represents a robust and credible method for neural population data analysis.

### 4.1  Advantages of the LapEig Method

The results presented in Figure 8 show that the LapEig model has better accuracy and stability in visualization applications than commonly used dimensionality-reduction methods. Neural activities reflect complex interactions between multiple variables and noise, so the information on a single sampling point can be easily disturbed. The common dimensionality-reduction methods tend to preserve the original signal characteristics as much as possible during the dimension-reduction process. Thus, the fluctuation caused by noise and other interfering factors is maintained in low-dimension information. This is why the results of the PCA, Isomap, and $t$-SNE methods show large within-class distance, especially in the fix epoch where there is no selectivity (see Figure 8). However, the LapEig method does not retain all the signal characteristics during the dimension-reduction process. The basic idea is to map the close points in a high-dimensional space to be as close as possible in a low-dimensional space. In the construction of a local neighborhood graph, only the edge weights between closely spaced points are greater than zero; the remaining weights are set to zero. In the low-dimensional embedding calculation, the distances between all connected nodes are similar in value. Therefore, if a fluctuation caused by the noise is not greater than a change caused by a behavior state, the distance between points remains close, and the difference between different motion states is more likely to be widened after dimension reduction. The results in Figure 3 show the obvious clustering performance of the LapEig method.

During the experiment, we also analyzed the classification effect of the LapEig method under different window lengths and different sampling frequencies. Figure S1 shows that when the time window for calculating the spike count was short, the fluctuation caused by the noise was relatively obvious and classification accuracy was low; when the time window was long, the fluctuation caused by the noise was relatively small and classification accuracy increased obviously. Based on these results, it can be concluded that the LapEig method not only could observe the changes in motion behaviors of the neural population effectively but could also suppress the influence of noise efficiently.

The results in Figure 9 show that in the experiment, the LapEig method was faster than the common nonlinear methods, which was due to the optimization algorithm of the LapEig model in low-dimensional embedding. In the dimension-reduction process, the LapEig method retained only the local feature information. The connection graph was generated by connecting the $k$ nearest neighbors of each point, and the connection graph matrix $W$ retained only the information on the connections. Therefore, compared with the Isomap and other methods that retained the global information, the connection graph matrix $W$ of the LapEig was a sparse matrix simple to calculate.

In addition, in the process of low-dimensional embedding of high-dimensional data, the optimization problem can be defined by $mintr(YTLY)s.t.YTDY=I$. The optimization problem can be further transformed into a generalized eigenvalue decomposition problem defined by $Ly=λDy$. Because the optimization problem is to find the eigenvectors of a sparse matrix, the operation speed of the LapEig is much faster. In general, the LapEig method is faster than the common manifold learning dimension-reduction methods.

### 4.2  Application of the LapEig Method

Neural activity usually reflects complex interactions between multiple variables and noise and often exhibits different responses at different behavior epochs (Donoghue, Suner, & Sanes, 1990; Sanes, Wang, & Donoghue, 1992; Moore, Nelson, & Sur, 1999; Hepp-Reymond, Kirkpatrick-Tanner, Gabernet, Qi, & Weber, 1999; Li, Padoa-Schioppa, & Bizzi, 2001; Tolias, Keliris, Smirnakis, & Logothetis, 2005; Stokes et al., 2013). Therefore, the analysis of the large-scale neural population data by using a single-neuron analysis method such as peristimulus time histogram is not only inefficient but also the internal connections between neural populations can be easily ignored. The commonly used population analysis methods such as motion classification (Hao et al., 2013) and motion trajectory fitting (Zhuang, Truccolo, Vargas-Irwin, & Donoghue, 2010; Pearce & Moran, 2012) can be regarded as supervised dimensionality-reduction methods. However, these methods can hardly display a single point in time. Therefore, the unsupervised dimensionality-reduction methods are increasingly applied to neural population analysis (Churchland et al., 2012, 2010, 2007; Mante et al., 2013; Lehky, Sereno, & Sereno, 2013; Rezai, Kleinhans, Matallanas, Selby, & Tripp, 2014; Mollazadeh, Aggarwal, Thakor, & Schieber, 2014; Dimitriadis et al., 2018; Vargas-Irwin et al., 2015). The proposed LapEig method also represents a kind of unsupervised dimensionality-reduction method.

Our results demonstrate that accurate motion decoding can be achieved by applying relatively simple methods, such as a $k$-nearest-neighbor-based classifier, to the LapEig method. In addition, the LapEig method not only can perform routine population analysis but also can observe more details, such as neural population change over a single trial. The LapEig method can provide a convenient visualization and quantitative evaluation of the experimental results. For instance, in Figure 3a, it can be observed that a few trials clearly deviated from the cluster center, so these trials interfered significantly and we could delete them. Moreover, in Figure 3b, the dimension-reduction results of the grasp, hook, and pinch were close, but the distance of the grip was far from the others. Consequently, we could reasonably predict that the firing pattern of a neural population when the four fingers were bent and in a circular shape with the thumb was significantly different from the pattern when the four fingers were straight and parallel with the thumb. In general, the LapEig method represents a powerful tool for neural population data analysis.

### 4.3  Limitations and Future Work

Although the LapEig method provides a visual review of the main data trends, some subtle trends are smoothed during the dimension-reduction process. Therefore, when neural activity needs to be studied further, a feasible method is to increase the dimension of retention and then analyze each dimension using the methods of a single unit analysis such as the peristimulus time histogram.

In the LapEig method implementation, the number of spikes within a time window is used as the original signal for dimensionality reduction, so the specific time information of the spikes is inevitably ignored. In the future, we will try to improve the proposed algorithm by using more accurate time signal processing methods, such as point processes and metric-space analysis, instead of the coarser spike count method to achieve more accurate neural signal acquisition, providing better dynamics visualization.

To demonstrate the application of the LapEig method in classification in this work, a simple classification method was used, and highly accurate classification results were achieved. Therefore, according to the results obtained in different experimental environments, the LapEig method can possibly be used as a brain-computer interface method to perform classification and decoding.

## 5  Conclusion

A proper understanding of neural population activity patterns is an essential step in human brain information processing. The LapEig method provides an effective tool for neural data analysis, allowing intuitive, dynamic visualization of spike signals. The proposed method is tested on two tasks using spike data, and it is shown that the LapEig method can effectively reveal the direct correspondence between neural dynamics and motion behavior. Compared with the common dimensionality-reduction methods, the LapEig method achieved higher clustering accuracy, smaller within-class clustering, and faster running speed in our experiments. Therefore, the LapEig method can be considered practical.

## Acknowledgments

We express our gratitude to Junming Zhu for his help in animal surgery, Shenglong Xiong for his help in monkey training, and Yimin Shen for his help in animal care. This work was supported by the National Key Research and Development Program of China (2017YFE0195500, 2017YFC1308500, 2014DFG32580), the National Natural Science Foundation of China (31627802, 31371001), and the Fundamental Research Funds for the Central Universities. The code we used in this letter has been uploaded to https://github.com/EveyZhang/Spike-Dimensionality-Reduction-by-Laplacian-Eigenmaps.

## References

,
D. A.
,
,
N. A.
,
Kosmidis
,
E. K.
, &
Theophilidis
,
G.
(
2010
).
Nass: An empirical approach to spike sorting with overlap resolution based on a hybrid noise-assisted methodology
.
J. Neurosci. Methods
,
190
(
1
),
129
142
.
,
D. A.
,
,
N. A.
,
Kosmidis
,
E. K.
, &
Theophilidis
,
G.
(
2012
).
In quest of the missing neuron: Spike sorting based on dominant-sets clustering
.
Computer Methods and Programs in Biomedicine
,
107
,
28
35
.
Ahrens
,
M. B.
,
Li
,
J. M.
,
Orger
,
M. B.
,
Robson
,
D. N.
,
Schier
,
A. F.
,
Engert
,
F.
, …
Portugues
,
R.
(
2012
).
Brain-wide neuronal dynamics during motor adaptation in zebrafish
.
Nature
,
485
(
7399
),
471
480
.
Ahrens
,
M. B.
,
Orger
,
M. B.
,
Robson
,
D. N.
,
Li
,
J. M.
, &
Keller
,
P. J.
(
2013
).
Whole-brain functional imaging at cellular resolution using light-sheet microscopy
.
Nature Methods
,
10
(
5
),
413
.
Alivisatos
,
A. P.
,
Chun
,
M.
,
Church
,
G. M.
,
Deisseroth
,
K.
,
Donoghue
,
J. P.
,
Greenspan
,
R. J.
, …
Yuste
,
R.
(
2013
).
The brain activity map
.
Science
,
339
(
6125
),
1284
1285
.
Belkin
,
M.
, &
Niyogi
,
P.
(
2003
).
Laplacian eigenmaps for dimensionality reduction and data representation
.
Neural Computation
,
15
(
6
),
1373
1396
.
Breakspear
,
M.
(
2017
).
Dynamic models of large-scale brain activity
.
Nature Neuroscience
,
20
(
3
),
340
.
Brunton
,
B. W.
,
Johnson
,
L. A.
,
Ojemann
,
J. G.
, &
Kutz
,
J. N.
(
2016
).
Extracting spatial-temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition
.
Journal of Neuroscience Methods
,
258
,
1
15
.
Chah
,
E.
,
Hok
,
V.
,
Della-Chiesa
,
A.
,
Miller
,
J. J. H.
,
O'Mara
,
S. M.
, &
Reilly
,
R. B.
(
2011
).
Automated spike sorting algorithm based on Laplacian eigenmaps and $k$-means clustering
.
Journal of Neural Engineering
,
8
(
1
),
016006
.
Chen
,
R.
,
Canales
,
A.
, &
Anikeeva
,
P.
(
2017
).
Neural recording and modulation technologies
.
Nature Reviews Materials
,
2
(
2
),
16093
.
Churchland
,
M. M.
,
Cunningham
,
J. P.
,
Kaufman
,
M. T.
,
Foster
,
J. D.
,
Nuyujukian
,
P.
,
Ryu
,
S. I.
, &
Shenoy
,
K. V.
(
2012
).
Neural population dynamics during reaching
.
Nature
,
487
(
7405
),
51
.
Churchland
,
M. M.
,
Cunningham
,
J. P.
,
Kaufman
,
M. T.
,
Ryu
,
S. I.
, &
Shenoy
,
K. V.
(
2010
).
Cortical preparatory activity: Representation of movement or first cog in a dynamical machine?
Neuron
,
68
(
3
),
387
400
.
Churchland
,
M. M.
,
Yu
,
B. M.
,
Sahani
,
M.
, &
Shenoy
,
K. V.
(
2007
).
Techniques for extracting single-trial activity patterns from large-scale neural recordings
.
Current Opinion in Neurobiology
,
17
(
5
),
609
618
.
Cui
,
G.
,
Jun
,
S. B.
,
Jin
,
X.
,
Pham
,
M. D.
,
Vogel
,
S. S.
,
Lovinger
,
D. M.
, &
Costa
,
R. M.
(
2013
).
Concurrent activation of striatal direct and indirect pathways during action initiation
.
Nature
,
494
(
7436
),
238
242
.
Cunningham
,
J. P.
, &
Yu
,
B. M.
(
2014
).
Dimensionality reduction for large-scale neural recordings
.
Nature Neuroscience
,
17
(
11
),
1500
1509
.
,
G.
,
Neto
,
J. P.
, &
Kampff
,
A. R.
(
2018
).
t-SNE visualization of large-scale neural recordings
.
Neural Computation
,
30
(
7
),
1750
1774
.
Donoghue
,
J. P.
,
Suner
,
S.
, &
Sanes
,
J. N.
(
1990
).
Dynamic organization of primary motor cortex output to target muscles in adult-rats. 2. Rapid reorganization following motor-nerve lesions
.
Experimental Brain Research
,
79
(
3
),
492
503
.
Ghanbari
,
Y.
,
Papamichalis
,
P.
, &
Spence
,
L.
(
2010
). Graph-spectrum-based neural spike features for stereotrodes and tetrodes. In
Proceedings of the IEEE International Conference on Acoustics Speech and Signal Processing
.
Piscataway, NJ
:
IEEE
.
Grewe
,
B. F.
,
Langer
,
D.
,
Kasper
,
H.
,
Kampa
,
B. M.
, &
Helmchen
,
F.
(
2010
).
High-speed in vivo calcium imaging reveals neuronal network activity with near-millisecond precision
.
Nature Methods
,
7
(
6
),
479
.
Hao
,
Y.
,
Zhang
,
Q.
,
Zhang
,
S.
,
Zhao
,
T.
,
Wang
,
Y.
,
Chen
,
W.
, &
Zheng
,
X.
(
2013
).
Decoding grasp movement from monkey premotor cortex for real-time prosthetic hand control
.
Chinese Science Bulletin
,
58
(
20
),
2512
2520
.
Hepp-Reymond
,
M. C.
,
Kirkpatrick-Tanner
,
M.
,
Gabernet
,
L.
,
Qi
,
H. X.
, &
Weber
,
B.
(
1999
).
Context-dependent force coding in motor and premotor cortical areas
.
Experimental Brain Research
,
128
(
1–2
),
123
133
.
Kayaert
,
G.
,
Biederman
,
I.
, &
Vogels
,
R.
(
2005
).
Representation of regular and irregular shapes in macaque inferotemporal cortex
.
Cerebral Cortex
,
15
(
9
),
1308
1321
.
Kerr
,
J. N. D.
, &
Denk
,
W.
(
2008
).
Imaging in vivo: Watching the brain in action
.
Nature Reviews Neuroscience
,
9
(
3
),
195
205
.
Kiani
,
R.
,
Esteky
,
H.
,
Mirpour
,
K.
, &
Tanaka
,
K.
(
2007
).
Object category structure in response patterns of neuronal population in monkey inferior temporal cortex
.
Journal of Neurophysiology
,
97
(
6
),
4296
4309
.
Kipke
,
D. R.
,
Shain
,
W.
,
Buzsaki
,
G.
,
Fetz
,
E.
,
Henderson
,
J. M.
,
Hetke
,
J. F.
, &
Schalk
,
G.
(
2010
).
Advanced neurotechnologies for chronic neural interfaces: New horizons and clinical opportunities
.
Journal of Neuroscience
,
28
(
46
),
11830
11838
.
Kozai
,
T. D.
,
Eles
,
J. R.
,
Vazquez
,
A. L.
, &
Cui
,
X. T.
(
2016
).
Two-photon imaging of chronically implanted neural electrodes: Sealing methods and new insights
.
Journal of Neuroscience Methods
,
258
,
46
55
.
Lehky
,
S. R.
, &
Sereno
,
A. B.
(
2007
).
Comparison of shape encoding in primate dorsal and ventral visual pathways
.
Journal of Neurophysiology
,
97
(
1
),
307
319
.
Lehky
,
S. R.
,
Sereno
,
M. E.
, &
Sereno
,
A. B.
(
2013
).
Population coding and the labeling problem: Extrinsic versus intrinsic representations
.
Neural Computation
,
25
(
9
),
2235
2264
.
Li
,
C.
,
,
C.
, &
Bizzi
,
E.
(
2001
).
Neuronal correlates of motor performance and motor learning in the primary motor cortex of monkeys adapting to an external force field
.
Neuron
,
30
(
2
),
593
607
.
Mank
,
M.
,
Santos
,
A. F.
,
Direnberger
,
S.
,
Mrsic-Flogel
,
T. D.
,
Hofer
,
S. B.
,
Stein
,
V.
, …
Griesbeck
,
O.
(
2008
).
A genetically encoded calcium indicator for chronic in vivo two-photon imaging
.
Nature Methods
,
5
(
9
),
805
811
.
Mante
,
V.
,
Sussillo
,
D.
,
Shenoy
,
K. V.
, &
Newsome
,
W. T.
(
2013
).
Context-dependent computation by recurrent dynamics in prefrontal cortex
.
Nature
,
503
(
7474
),
78
.
,
M.
,
Aggarwal
,
V.
,
Thakor
,
N. V.
, &
Schieber
,
M. H.
(
2014
).
Principal components of hand kinematics and neurophysiological signals in motor cortex during reach to grasp movements
.
Journal of Neurophysiology
,
112
(
8
),
1857
1870
.
Moore
,
C. I.
,
Nelson
,
S. B.
, &
Sur
,
M.
(
1999
).
Dynamics of neuronal processing in rat somatosensory cortex
.
Trends in Neurosciences
,
22
(
11
),
513
520
.
Op De
Beeck
,
H.
,
Wagemans
,
J.
, &
Vogels
,
R.
(
2001
).
Inferotemporal neurons represent low-dimensional configurations of parameterized shapes
.
Nature Neuroscience
,
4
(
12
),
1244
1252
.
Pearce
,
T. M.
, &
Moran
,
D. W.
(
2012
).
Strategy-dependent encoding of planned arm movements in the dorsal premotor cortex
.
Science
,
337
(
6097
),
984
988
.
Raos
,
V.
,
Umiltá
,
M. A.
,
Gallese
,
V.
, &
Fogassi
,
L.
(
2004
).
Functional properties of grasping-related neurons in the dorsal premotor area F2 of the macaque monkey
.
J. Neurophysiol.
,
92
(
4
),
1990
2002
.
Rezai
,
O.
,
Kleinhans
,
A.
,
Matallanas
,
E.
,
Selby
,
B.
, &
Tripp
,
B. P.
(
2014
).
Modeling the shape hierarchy for visually guided grasping
.
Frontiers in Computational Neuroscience
,
8
(
132
).
Richardson
,
M. J. E.
(
2008
).
Spike-train spectra and network response functions for non-linear integrate-and-fire neurons
.
Biological Cybernetics
,
99
(
4–5
),
381
392
.
Rigotti
,
M.
,
Barak
,
O.
,
Warden
,
M. R.
,
Wang
,
X.
,
Daw
,
N. D.
,
Miller
,
E. K.
, &
Fusi
,
S.
(
2013
).
The importance of mixed selectivity in complex cognitive tasks
.
Nature
,
497
(
7451
),
585
590
.
Rios
,
G.
,
Lubenov
,
E. V.
,
Chi
,
D.
,
Roukes
,
M. L.
, &
Siapas
,
A. G.
(
2016
).
Nanofabricated neural probes for dense 3-D recordings of brain activity
.
Nano Letters
,
16
(
11
),
6857
6862
.
Rolls
,
E. T.
, &
Tovee
,
M. J.
(
1995
).
Sparseness of the neuronal representation of stimuli in the primate temporal visual-cortex
.
Journal of Neurophysiology
,
73
(
2
),
713
726
.
Sanes
,
J. N.
, &
Donoghue
,
J. P.
(
2000
).
Plasticity and primary motor cortex
.
Annual Review of Neuroscience
,
23
,
393
415
.
Sanes
,
J. N.
,
Wang
,
J.
, &
Donoghue
,
J. P.
(
1992
).
Immediate and delayed changes of rat motor cortical output representation with new forelimb configurations
.
Cerebral Cortex
,
2
(
2
),
141
152
.
Song
,
Y.
,
Wang
,
Y.
, &
Viventi
,
J.
(
2017
).
Unsupervised learning of spike patterns for seizure detection and wavefront estimation of high resolution micro electrocorticographic (mu ECoG) data
.
IEEE Transactions on Nanobioscience
,
16
(
6
),
418
427
.
Spira
,
M. E.
, &
Hai
,
A.
(
2013
).
Multi-electrode array technologies for neuroscience and cardiology
.
Nature Nanotechnology
,
8
(
2
),
83
94
.
Stokes
,
M. G.
,
Kusunoki
,
M.
,
Sigala
,
N.
,
Nili
,
H.
,
Gaffan
,
D.
, &
Duncan
,
J.
(
2013
).
Dynamic coding for cognitive control in prefrontal cortex
.
Neuron
,
78
(
2
),
364
375
.
Suzuki
,
M.
(
2015
).
Manifold learning and convergence of laplacian eigenmaps
.
PhD diss., McGill University
.
Tolias
,
A. S.
,
Keliris
,
G. A.
,
Smirnakis
,
S. M.
, &
Logothetis
,
N. K.
(
2005
).
Neurons in macaque area V4 acquire directional tuning after adaptation to motion stimuli
.
Nature Neuroscience
,
8
(
5
),
591
593
.
Vargas-Irwin
,
C. E.
,
Brandman
,
D. M.
,
Zimmermann
,
J. B.
,
Donoghue
,
J. P.
,
Black
,
M. J.
(
2015
).
Spike train similarity space (SSIMS): A framework for single neuron and ensemble data analysis
.
Neural Computation
,
27
(
1
),
1
31
.
Viventi
,
J.
,
Kim
,
D.
,
Vigeland
,
L.
,
Frechette
,
E. S.
,
Blanco
,
J. A.
,
Kim
,
Y.
, …
Litt
,
B.
(
2011
).
Flexible, foldable, actively multiplexed, high-density electrode array for mapping brain activity in vivo
.
Nature Neuroscience
,
14
(
12
),
1138
1599
.
Young
,
M. P.
, &
Yamane
,
S.
(
1992
).
Sparse population coding of faces in the inferotemporal cortex
.
Science
,
256
(
5061
),
1327
1331
.
Zhuang
,
J.
,
Truccolo
,
W.
,
Vargas-Irwin
,
C.
, &
Donoghue
,
J. P.
(
2010
).
Decoding 3-d reach and grasp kinematics from high-frequency local field potentials in primate primary motor cortex
.
IEEE Transactions on Biomedical Engineering
,
57
(
7
),
1774
1784
.

## Author notes

*

Shaomin Zhang is the corresponding author.