To demonstrate the effectiveness of our model (distill-variational autoencoder [d-VAE]) in extracting behaviorally relevant signals, we conducted experiments on the synthetic dataset where the ground truth of relevant and irrelevant signals are already known (see Methods) and three benchmark datasets with different paradigms (Figure 2a, e, and i; see Methods for details), and compared d-VAE with four other distillation models, including pi-VAE (Zhou, 2020), PSID (Sani et al., 2021), TNDM (Hurwitz, 2021), and LFADS (Pandarinath et al., 2018). Specifically, we first applied these distillation models to raw signals to obtain the distilled behaviorally relevant signals, considering the residuals as behaviorally irrelevant signals. We then evaluated the decoding ${\displaystyle {\mathrm{R}}^2}$ between the predicted velocity and actual velocity of the two partition signals using a linear Kalman filter (KF) and a nonlinear artificial neural network (ANN) and measured the neural similarity between behaviorally relevant and raw signals.

Figure 2 with 1 supplement see all

Download asset Open asset

Overall, d-VAE successfully extracts effective behaviorally relevant signals that meet the three criteria outlined above on both synthetic (Figure 2—figure supplement 1) and real data (Figure 2). On the synthetic data (Figure 2—figure supplement 1), results show that d-VAE can strike an effective balance between the reconstruction and decoding performance of generated signals to extract effective relevant signals that are similar to the ground truth relevant signals, meanwhile removing effective irrelevant signals that resemble the ground truth irrelevant signals (Figure 2—figure supplement 1a–g), and outperforms other distillation models (Figure 2—figure supplement 1h–k). On the real data, specifically in the obstacle avoidance task (Figure 2a), the monkey is required to move the ball from the start point (blue) to the target point (yellow) without hitting the obstacle. For the decoding performance of behaviorally relevant signals (Figure 2b), the signals distilled by d-VAE outperform the raw signals (purple bars with slash lines) and the signals distilled by all other distillation models (PSID, pink; pi-VAE, green; TNDM, blue; and LFADS, light green) with the KF as well as the ANN. For the decoding performance of behaviorally irrelevant signals (Figure 2c), behaviorally irrelevant signals obtained by d-VAE achieves the lowest decoding performance compared with behaviorally irrelevant signals obtained by other approaches. Therefore, the combination of dVAE’s highest decoding performance for behaviorally relevant signals and lowest decoding performance for behaviorally irrelevant signals demonstrate its superior ability to extract behaviorally relevant signals from noisy signals. For the neural similarity between behaviorally relevant and raw signals (Figure 2d), the distilled signals obtained by d-VAE achieve the highest performance among competitors (p<0.01, Wilcoxon rank-sum test). Similar results were obtained for the center-out task (Figure 2e–h) and the self-paced reaching task (Figure 2i–l), indicating the consistency of d-VAE’s distillation ability across a range of motor tasks. To provide a more intuitive illustration of the similarity between raw and distilled signals, we displayed the firing rate of neuronal activity in five trials under the same condition (Figure 2m), and results clearly show that the firing pattern of distilled signals is similar to the corresponding raw signals.

In summary, d-VAE distills effective behaviorally relevant signals that preserve behavioral information maximally and are similar to raw signals. Meanwhile, the behaviorally irrelevant signals discarded by d-VAE contain a little behavioral information. Therefore, these signals are reliable for exploring the encoding and decoding mechanisms of relevant signals.

Following signal separation, we first explored how behaviorally irrelevant signals affect the analysis of neural activity at the single-neuron level. Specifically, we examined the effect of irrelevant signals on two critical properties of neuronal activity: the preferred direction (PD) (Georgopoulos et al., 1986) and trial-to-trial variability. Our objective was to know how irrelevant signals affect the PD of neurons and whether irrelevant signals contribute significantly to neuronal variability.

To explore how irrelevant signals affect the PD of neurons, we first calculated the PD of both raw and distilled signals separately and then quantified the PD deviation by the angle difference between these two signals. Results show that the PD deviation increases as the neuronal ${\displaystyle {\mathrm{R}}^2}$ decreases (red curve, Figure 3a and e and Figure 3—figure supplement 1a). It is worth noting that when using ${\displaystyle {\mathrm{R}}^2}$ to describe neurons, it indicates the extent to which neuronal activity is explained by the linear encoding model (Collinger et al., 2013; Wodlinger et al., 2015). Neurons with larger ${\displaystyle {\mathrm{R}}^2}$ (strongly linear-tuned neurons) exhibit stable PDs with signal distillation (see example PDs in the inset), while neurons with smaller ${\displaystyle {\mathrm{R}}^2}$ (weakly linear-tuned neurons) show a larger PD deviation. These results indicate that irrelevant signals have a small effect on strongly tuned neurons but a large effect on weakly tuned neurons. One possible reason for the larger PD deviation in weakly tuned neurons is that they have a lower degree of linear encoding but a higher degree of nonlinear encoding, and highly nonlinear structures are more susceptible to interference from irrelevant signals (Nogueira et al., 2023). Moreover, after filtering out the behaviorally irrelevant signals, the cosine tuning fit (${\displaystyle {\mathrm{R}}^2}$) of neurons increases (p<10^-20, Wilcoxon signed-rank test; Figure 3b and f and Figure 3—figure supplement 1b), indicating that irrelevant signals reduce the neurons’ tuning expression. Notably, even after removing the interference of irrelevant signals, the ${\displaystyle {\mathrm{R}}^2}$ of neurons remains relatively low and varies among neurons. These results demonstrate that the linear encoding model only explains a small fraction of neural responses, and neuronal activity encodes behavioral information in complex nonlinear ways.

Figure 3 with 2 supplements see all

Download asset Open asset

To investigate whether irrelevant signals significantly contribute to neuronal variability, we compared the neuronal variability (measured with the Fano factor [FF]; Churchland et al., 2010) of relevant and raw signals. Results show that the condition-averaged FF of each neuron of distilled signals is lower than that of raw signals (p<10^-20, Wilcoxon signed-rank test; Figure 3c and g), and the mean (broken line) and median FFs of all neurons under different conditions are also significantly lower than those of raw signals (p<0.01, Wilcoxon signed-rank test; Figure 3d and h), indicating that irrelevant signals significantly contribute to neuronal variability. We then visualized the single-trial neuronal activity of example neurons under different conditions (Figure 3i and Figure 3—figure supplement 2). Results demonstrate that the patterns of relevant signals are more consistent and stable across different trials than raw signals, and the firing activity of irrelevant signals varies randomly. These results indicate that irrelevant signals significantly contribute to neuronal variability, and eliminating the interference of irrelevant signals enables us to observe the changes in neural pattern more accurately.

The neural population structure is an essential characteristic of neural activity. Here, we examined how behaviorally irrelevant signals affect the analysis of neural activity at the population level, including four aspects: (1) the population properties of relevant and irrelevant signals, (2) the subspace overlap relationship between the two signal components, (3) how the two partitions contribute to raw signals, and (4) the difference in population properties between raw and distilled signals.

To explore the population properties of relevant and irrelevant signals, we separately applied principal component analysis (PCA) on each partition to obtain the corresponding cumulative variance curve in a descending variance order. Our results show that the primary subspace (capturing the top 90% variance) of relevant signals (thick red line, Figure 4a and e and Figure 4—figure supplement 1a) is only explained by a few dimensions (7, 13, and 9 for each dataset), indicating that the primary part of behaviorally relevant signals exists in a low-dimensional subspace. In contrast, the primary subspace of irrelevant signals (thick blue line, Figure 4b and f and Figure 4—figure supplement 1b) requires more dimensions (46, 81, and 59). The variance distribution of behaviorally irrelevant signals across dimensions (thick blue line, Figure 4b and f and Figure 4—figure supplement 1b) is more even than behaviorally relevant signals (thick red line, Figure 4a and e and Figure 4—figure supplement 1a) but not as uniform as Gaussian noise ${\displaystyle \mathcal{N}(\mathbf{0},\mathit{I})}$ (thin gray line, Figure 4b and f and Figure 4—figure supplement 1a), indicating that irrelevant signals are not pure noise but rather bear some significant structure, which may represent information from other irrelevant tasks.

Figure 4 with 4 supplements see all

Download asset Open asset

To investigate the subspace overlap between relevant and irrelevant signals, we calculated how many variances of irrelevant signals can be captured by relevant PCs by projecting irrelevant signals onto relevant PCs and vice versa (Elsayed et al., 2016; Rouse and Schieber, 2018; Jiang et al., 2020) (see Methods). We found that the variance of irrelevant signals increases relatively uniformly over relevant PCs (thin blue line, Figure 4a and e and Figure 4—figure supplement 1a), like random noise’s variance accumulation explained by relevant PCs (thin gray line, Figure 4a and e and Figure 4—figure supplement 1a); similar results are observed for relevant signals explained by irrelevant PCs (thin red line, Figure 4b and f and Figure 4—figure supplement 1b). These results indicate that relevant PCs cannot match the informative dimensions of irrelevant signals and vice versa, suggesting the dimensions of behaviorally relevant and irrelevant signals are unrelated. It is worth mentioning that the signals obtained by pi-VAE are in contrast to our findings. Its results show that a few relevant PCs can explain a considerable variance of irrelevant signals (thin red line, Figure 4—figure supplement 2b, f, j), which indicates that the relevant and irrelevant PCs are closely related. The possible reason is that the pi-VAE leaves many relevant signals within the irrelevant signals. Notably, Figure 4a and e and Figure 4—figure supplement 1a show that the behaviorally relevant primary subspace captures only a minor portion of the variance from irrelevant signals when they are projected onto it (9%, 12%, and 9%), indicating that the primary subspace of behaviorally relevant signals is nearly orthogonal to irrelevant space.

To investigate the composition of raw signals by the two partitions, we performed PCA on raw neural signals to obtain raw PCs, and then projected the relevant and irrelevant signals onto these PCs to assess the proportion of variance of raw signals explained by each partition. First, we analyzed the overall proportion of relevant and irrelevant signals that constitute the raw signals (the inset pie plot, Figure 4c and g and Figure 4—figure supplement 1c). The variance of the raw signals is composed of three parts: the variance of relevant signals, the variance of irrelevant signals, and the correlation between relevant and irrelevant signals (see Methods). The results demonstrate that the irrelevant signals account for a large proportion of raw signals, suggesting the motor cortex encodes considerable information that is not related to the measured behaviors. In addition, there is only a weak correlation between relevant and irrelevant signals, implying that behaviorally relevant and irrelevant signals are nearly uncorrelated.

We then examined the proportions of relevant and irrelevant signals in each PC of raw signals. We found that relevant signals (red) occupy the dominant proportions in the larger variance raw PCs (before the PC marked with a red triangle), while irrelevant signals (blue) occupy the dominant proportions in the smaller variance raw PCs (after the PC marked with a red triangle) (Figure 4c and g and Figure 4—figure supplement 1c). Similar results are observed in the accumulation of each raw PC (Figure 4d and h and Figure 4—figure supplement 1d). Specifically, the results show that the variance accumulation of raw signals (purple line) in larger variance PCs is mainly contributed by relevant signals (red line), while irrelevant signals (blue line) primarily contribute to the lower variance PCs. These results demonstrate that irrelevant signals have a small effect on larger variance raw PCs but a large effect on smaller variance raw PCs. This finding eliminates the concern that irrelevant signals would significantly affect the top few PCs of raw signals and thus produce inaccurate conclusions. To further validate this finding, we used the top six PCs as jPCA (Churchland et al., 2012) did to examine the rotational dynamics of distilled and raw signals (Figure 4—figure supplement 3). Results show that the rotational dynamics of distilled signals are similar to those of raw signals.

Finally, to directly compare the population properties of raw and relevant signals, we plotted the cumulative variance curves of raw and relevant signals (Figure 4—figure supplement 4). Results (upper left corner curves, Figure 4—figure supplement 4) show that the cumulative variance curve of relevant signals (red line) accumulates faster than that of raw signals (purple line) in the preceding larger variance PCs, indicating that the variance of the relevant signal is more concentrated in the larger variance PCs than that of raw signals. Furthermore, we found that the dimensionality of primary subspace of raw signals (26, 64, and 45 for datasets A, B, and C) is significantly higher than that of behaviorally relevant signals (7, 13, and 9), indicating that using raw signals to estimate the neural dimensionality associated with behaviors leads to an overestimation.

The results presented above regarding PDs (Figure 3 and Figure 3—figure supplement 1) demonstrate that irrelevant signals significantly impact smaller ${\displaystyle {\mathrm{R}}^2}$ neurons and weakly impact larger ${\displaystyle {\mathrm{R}}^2}$ neurons. Under the interference of irrelevant signals, it is difficult to explore the amount of behavioral information in neuronal activity. Given that we have accurately separated the behaviorally relevant and irrelevant signals, we explored whether irrelevant signals would mask some encoded information of neuronal activity, especially for smaller ${\displaystyle {\mathrm{R}}^2}$ neurons.

To answer the question, we divided the neurons into two groups of smaller ${\displaystyle {\mathrm{R}}^2}$ (${\displaystyle {\mathrm{R}}^2}$<=0.03) and larger ${\displaystyle {\mathrm{R}}^2}$ (${\displaystyle {\mathrm{R}}^2}$>0.03), and then used decoding models to assess how much information is encoded in raw and distilled signals. As shown in Figure 5a, for the smaller ${\displaystyle {\mathrm{R}}^2}$ neuron group, both KF and ANN decode behavioral information poorly on raw signals, but achieve high decoding performance using relevant signals. Specifically, the KF decoder (left plot, Figure 5a) improves the decoding ${\displaystyle {\mathrm{R}}^2}$ significantly from 0.044 to 0.616 (improves by about 1300%, Wilcoxon rank-sum test) after signal distillation; the ANN decoder (right plot, Figure 5a) improves from 0.374 to 0.753 (improves by about 100%, Wilcoxon rank-sum test). For the larger ${\displaystyle {\mathrm{R}}^2}$ neuron group, the decoding performance of relevant signals with ANN does not improve much compared with the decoding performance of raw signals, but the decoding performance of relevant signals with KF is significantly better than that of raw signals (p<0.01, Wilcoxon rank-sum test). Similar results are obtained with datasets B (Figure 5d) and C (Figure 5—figure supplement 1a). These results indicate that irrelevant signals mask behavioral information encoded by neuronal populations, especially for smaller ${\displaystyle {\mathrm{R}}^2}$ neurons with a higher degree of nonlinearity, and that smaller ${\displaystyle {\mathrm{R}}^2}$ neurons actually encode rich behavioral information.

Figure 5 with 2 supplements see all

Download asset Open asset

The fact that the smaller ${\displaystyle {\mathrm{R}}^2}$ neurons encode rich information seems unexpected, and interestingly, we cannot obtain this rich information solely by distilling smaller ${\displaystyle {\mathrm{R}}^2}$ neurons. This observation gives rise to two alternative scenarios. The first is that larger ${\displaystyle {\mathrm{R}}^2}$ neurons introduce additional signals to smaller ${\displaystyle {\mathrm{R}}^2}$ neurons, which they do not inherently possess, resulting in an excessive amount of behavioral information within the smaller ${\displaystyle {\mathrm{R}}^2}$ neurons. The second is that the smaller ${\displaystyle {\mathrm{R}}^2}$ neurons inherently possess a substantial amount of information, and larger ${\displaystyle {\mathrm{R}}^2}$ neurons utilize their neural activity, which is correlated with that of small ${\displaystyle {\mathrm{R}}^2}$ neurons, to aid in restoring the small ${\displaystyle {\mathrm{R}}^2}$ neurons’ original appearance; this process is analogous to image denoising, where damaged noisy pixels necessitate the assistance of their correlated, clean neighboring pixels to recover their original appearance. We initially tested the first scenario and found it to be unsupported for two key reasons. First, our model enforces that distilled neuronal activity closely resembles the corresponding original neuronal activity, effectively preventing the generation of arbitrarily shaped neuronal activity, such as that of other neurons. As shown in Figure 3i and Figure 3—figure supplement 2, our distilled relevant neuronal activity exhibits a high degree of similarity to the corresponding raw neuronal activity. To assess whether the distilled neurons exhibit the highest similarity to the corresponding raw neurons, we compared the neural similarity (${\displaystyle {\mathrm{R}}^2}$) of each distilled neuron to all raw neurons. The results indicate that 78/90 (87%, dataset A), 153/159 (96%, dataset B), and 91/91 (100%, dataset C) distilled neurons are most similar to the corresponding neurons. The remaining distilled neurons rank among the top four in similarity to the corresponding neurons, further confirming the close resemblance of distilled neuronal activity to the corresponding raw neuronal activity. Second, as we emphasized in the section on validating behaviorally relevant signals with the second criterion, if this large amount of information is compensated by other neurons, the residuals should also contain a large amount of information. However, as illustrated in Figure 2c, g, and k, the residuals contain only little information. Therefore, based on these two reasons, the first scenario is rejected. Then, we tested the second scenario. To verify this scenario, we conducted experiments using synthetic data with known ground truth (see Methods). In this dataset, small ${\displaystyle {\mathrm{R}}^2}$ neurons inherently contained a substantial amount of information but were obscured by noise, making them undecodable. We aimed to assess whether d-VAE could recover the lost information and restore the damaged neuronal activity. The results demonstrate that, with the assistance of large ${\displaystyle {\mathrm{R}}^2}$ neurons, d-VAE effectively recovers a significant amount of information that is obscured by noise (Figure 5—figure supplement 2a). Additionally, the distilled signals exhibit a remarkable improvement in neural similarity to the ground truth signals compared to the raw signals (p<0.01, Wilcoxon rank-sum test; Figure 5—figure supplement 2b). Therefore, these results support the second scenario and collectively confirm that smaller ${\displaystyle {\mathrm{R}}^2}$ neurons indeed contain rich behavioral information, and this finding is not a by-product of d-VAE.

Given that both smaller and larger ${\displaystyle {\mathrm{R}}^2}$ neurons encode rich behavioral information, it is worth noting that the sum of the decoding performance of smaller ${\displaystyle {\mathrm{R}}^2}$ neurons and larger ${\displaystyle {\mathrm{R}}^2}$ neurons is significantly greater than that of all neurons for relevant signals (red bar, Figure 5a and d and Figure 5—figure supplement 1a), demonstrating that movement parameters are encoded very redundantly in neuronal population. In contrast, we cannot find this degree of neural redundancy in raw signals (purple bar, Figure 5a and d and Figure 5—figure supplement 1a) because the encoded information of smaller ${\displaystyle {\mathrm{R}}^2}$ neurons are masked by irrelevant signals. Therefore, these smaller ${\displaystyle {\mathrm{R}}^2}$ neurons, which are usually ignored, are actually useful and play a critical role in supporting neural redundancy. Generally, cortical redundancy can arise from neuronal correlations, which are critical for revealing certain aspects of neural circuit organization (Yatsenko et al., 2015). Accordingly, we visualized the ordered correlation matrix of neurons (see Methods) for both raw and relevant signals (Figure 5b and e and Figure 5—figure supplement 1b) and found that the neuronal correlation of relevant signals is stronger than that of raw signals. These results demonstrate that irrelevant signals weaken the neuronal correlation, which may hinder the accurate investigation of neural circuit organization.

Considering the rich redundancy and strong correlation of neuronal activity, we wondered whether the neuronal population could utilize redundant information from other neurons to exhibit robustness under the perturbation of neuronal destruction. To investigate this question, we evaluated the decoding performance of dropping out neurons from larger ${\displaystyle {\mathrm{R}}^2}$ to smaller ${\displaystyle {\mathrm{R}}^2}$ on raw and relevant signals. The results (Figure 5c and f and Figure 5—figure supplement 1c) show that the decoding performance of the KF and ANN on raw signals (purple line) decreases steadily before the number of neurons marked (vertical gray line), and the remaining smaller ${\displaystyle {\mathrm{R}}^2}$ neurons decode behavioral information poorly. In contrast, even if many neurons are lost, the decoding performance of KF and ANN on relevant signals (red line) maintains high accuracy. This finding indicates that behaviorally relevant signals are robust to the disturbance of neuron drop-out, and smaller ${\displaystyle {\mathrm{R}}^2}$ neurons play a critical role in compensating for the failure of larger ${\displaystyle {\mathrm{R}}^2}$ neurons. In contrast, this robustness cannot be observed in raw signals because irrelevant signals mask neurons’ information and weaken their correlation. Notably, the ANN outperforms the KF when only smaller ${\displaystyle {\mathrm{R}}^2}$ neurons are left (Figure 5c and f and Figure 5—figure supplement 1c), suggesting that smaller ${\displaystyle {\mathrm{R}}^2}$ neurons can fully exploit their nonlinear ability to cope with large-scale neuronal destruction.

The results presented above regarding subspace overlap (Figure 4 and Figure 3—figure supplement 1) show that irrelevant signals have a small impact on larger variance PCs but dominate smaller variance PCs. Therefore, we aimed to investigate whether irrelevant signals would mask some encoded information of neural population, especially signals composed of smaller variance PCs.

To answer the question, we compared the decoding performance of raw and distilled signals with different raw PC groups. Specifically, we first divided the raw PCs into two groups, i.e., smaller variance PCs and larger variance PCs, defined by ratio of relevant to irrelevant signals in the raw PCs (the red triangle, see Figure 4c and g and Figure 3—figure supplement 1c). Then, we projected raw and distilled signals onto these two PC groups and got the corresponding signals. Results show that, for the smaller variance PC group, both KF and ANN achieve much better performance on distilled signals than raw signals (p<0.01, Wilcoxon rank-sum test, for ANN), whereas for the larger variance PC group, the decoding performance of relevant signals does not improve a lot compared with the decoding performance of raw signals (see Figure 6a and d and Figure 5—figure supplement 1d). These results demonstrate that irrelevant signals mask the behavioral information encoded by different PC groups, especially for signals composed of smaller variance PCs (smaller variance PC signals), and smaller variance PC signals actually encode rich behavioral information.

Figure 6 with 1 supplement see all

Download asset Open asset

The above results are based on raw PCs. However, raw PCs are biased by irrelevant signals and thus cannot faithfully reflect the characteristics of relevant signals. As we have successfully separated the behaviorally relevant signals, we aimed to explore how behavioral information of distilled signals is distributed across relevant PCs. To do so, we used decoding models to evaluate the amount of behavioral information contained in cumulative PCs of relevant signals (using raw signals as a comparison). The cumulative variance explained by PCs in descending and ascending order of variance and the dimensionality corresponding to the top 90% variance signals (called primary signals) and the last 10% variance signals (called secondary signals) are shown in Figure 4—figure supplement 4.

Here, we first investigated secondary signals’ decoding ability solely by accumulating PCs from smaller to larger variance. The results show that, for relevant signals, KF can hardly decode behavioral information solely using secondary signals (red line; left plot, Figure 6b and e and Figure 5—figure supplement 1e), but ANN can decode rich information (red line; right plot, Figure 6b and e and Figure 5—figure supplement 1e). These results indicate that smaller variance PC signals encode rich information in complex nonlinear ways. In contrast, when using raw signals composed of the same number of dimensions as the secondary signals (purple line, Figure 6b and e and Figure 5—figure supplement 1e), the amount of information identified by ANN is significantly smaller than that of relevant secondary signals (p<0.01, Wilcoxon rank-sum test). These results demonstrate that signals composed of these neural dimensions actually encode rich behavioral information, and irrelevant signals make them seem insignificant, indicating that behavioral information is distributed in a higher-dimensional subspace than expected from raw signals.

We then investigated the effect of superimposing secondary signals on primary signals by accumulating PCs from larger to lower variance. The results (Figure 6c and f and Figure 5—figure supplement 1f) show that secondary signals improve the decoding performance of ANN a little but improve the decoding performance of KF a lot. The discrepancy between the two decoders reflects their different abilities to utilize the information within the signal. KF cannot use the nonlinear information in primary signals as effectively as ANN can and thus require secondary signals to improve decoding performance. Notably, KF shows steady growth in decoding performance on relevant signals across 10–30 dimensions, and requires approximately 30–40 dimensions to achieve performance saturation. These results demonstrate that these smaller variance PC signals actually encode behavioral information, and suggest that behavioral information exists in a higher-dimensional subspace than anticipated from raw signals. Interestingly, we can find that although secondary signals nonlinearly encode behavioral information and are decoded poorly by linear decoders, they considerably improve KF performance by superimposing on primary signals (left plot, Figure 6c and f and Figure 5—figure supplement 1f); and the sum of the sole decoding performance of primary and secondary signals is lower than the decoding performance of full signals. These results indicate that the combination of smaller and larger variance PCs produces a synergy effect (Narayanan et al., 2005), enabling secondary signals that cannot be linearly decoded to improve the linear decoding performance.

Finally, considering the substantial enhancement in KF decoding performance when superimposing the secondary signals on the primary ones, we explored which aspect of movement parameters was most improved. In BMIs, directional control has achieved great success (Georgopoulos et al., 1986; Hochberg et al., 2012), but precise speed control, especially at lower speeds such as hold or stop, has always been challenging (Wodlinger et al., 2015; Inoue et al., 2018). Thus, we hypothesized that these signals might improve the lower-speed velocity. To test this, we divided samples into lower-speed and higher-speed regions and assessed which region improved the most by superimposing the secondary signals (see details in Methods). After superimposing the secondary signals, the absolute improvement ratio of the lower-speed region is significantly higher than that of the higher-speed region (p<0.05, Wilcoxon rank-sum test; Figure 6—figure supplement 1a, b, and c). Furthermore, we visualized the relative improvement ratio of five example trials for the two regions, and the results (Figure 6—figure supplement 1d) demonstrate that secondary signals significantly improve the estimation of lower speed. These results demonstrate that the secondary signals enhance the lower-speed control, suggesting that smaller variance PC signals may be involved in regulating precise motor control.

Understanding the readout mechanism of the motor cortex is crucial for both neuroscience and neural engineering, which remains unclear. By filtering out the interference of behaviorally irrelevant signals, we found a stunning result: the linear decoder KF achieves comparable performance to that of the nonlinear decoder ANN (p=0.10, 0.15, and 0.55 for datasets A, B, and C, Wilcoxon rank-sum test; Figure 2b, f, and j). Considering the decoding performance and model complexity (the simplicity principle, also called Occam’s razor), movement behaviors are more likely to be generated by the linear readout, suggesting linear readout may be performed in the motor cortex.

Given the significant improvement in linear decoding performance, one might doubt that it is our distillation model that makes signals that are inherently nonlinearly decodable become linearly decodable. In practice, this situation does not hold for two reasons. First, our criterion that irrelevant signals should contain minimal information can effectively exclude this situation. Specifically, if this situation occurs, the model would significantly modify the structure of the generated signals, causing a deviation from the structure of the ground truth signals. Consequently, these uncharacterized or modified ground truth signals would remain within the residuals, resulting in residuals that contain a substantial amount of information. To illustrate this, consider an example where ${\displaystyle z=x+y=n^2+y}$, with ${\displaystyle z}$, ${\displaystyle x}$, ${\displaystyle y}$, and ${\displaystyle n}$ representing raw signals, relevant signals, irrelevant signals, and behavioral variables, respectively. If the distilled relevant signals are ${\displaystyle \hat{x}=n}$, the corresponding irrelevant signals are ${\displaystyle n^2-n+z}$. Clearly, the distilled signal can be linearly decoded, but this results in the residuals containing a large amount of information. However, as demonstrated in Figure 2c, g, and k, the irrelevant signals obtained by d-VAE only contain little information, thus excluding this situation. Second, our synthetic experiments offer additional evidence supporting the conclusion that d-VAE does not make inherently nonlinearly decodable signals become linearly decodable ones. As depicted in Figure 5—figure supplement 2c, there exists a significant performance gap between KF and ANN when decoding the ground truth signals of smaller ${\displaystyle {\mathrm{R}}^2}$ neurons (p<0.01, Wilcoxon rank-sum test). KF exhibits notably low performance, leaving substantial room for compensation by d-VAE. However, following processing by d-VAE, KF’s performance of distilled signals fails to surpass its already low ground truth performance and remains significantly inferior to ANN’s performance (p<0.01, Wilcoxon rank-sum test). These results collectively confirm that our approach does not convert signals that are inherently nonlinearly decodable into linearly decodable ones.

In summary, these findings demonstrate that behaviorally irrelevant signals significantly complicate the readout of behavioral information and provide compelling evidence supporting the notion that the motor cortex may use a linear readout mechanism to generate movement behaviors.

Behaviorally relevant patterns can be extracted either at single-neuron or latent neural population levels. In our study, we focused on the single-neuron level, aiming to preserve the underlying properties of individual neurons. By maintaining the properties of each neuron, researchers can investigate how the neuronal population performs when one of the neurons is destroyed. This kind of analysis is particularly useful in closed-loop stimulation experiments that use electrophysiological (Sun et al., 2022) or optogenetic (Zhang et al., 2023) interventions. Furthermore, behaviorally relevant signals also allow for population-level analysis and provide clean benchmark signals to test and compare the variance capture ability of different hypothesis-driven models.

At the single-neuron level, it is common practice to use trial-averaged neuronal responses of the same task parameters to analyze neural mechanisms (Kobak et al., 2016; Rouse and Schieber, 2018). However, trial averaging sacrifices single-trial information, thereby providing an incomplete characterization of neural activity. Furthermore, trial-averaged responses still contain a significant amount of behaviorally irrelevant signals caused by uninstructed movements (Musall et al., 2019), which can lead to a contaminated version of behaviorally relevant signals. In contrast, our model is capable of extracting clean behaviorally relevant neural activity for every single trial. At the latent population level, existing latent variable models (Sani et al., 2021; Pandarinath et al., 2018; Yu, 2008; Zhou, 2020; Hurwitz, 2021) focus on modeling some specific properties of latent population representations, such as linear or nonlinear dynamics (Sani et al., 2021; Pandarinath et al., 2018; Churchland et al., 2012; Hurwitz, 2021), temporal smoothness (Yu, 2008), and interpretability (Zhou, 2020). Since these models make restrictive assumptions involving characterizing specific neural properties, they fail to capture other parts of behaviorally relevant signals that do not meet their assumptions, providing no guarantee that the generated signals preserve behavioral information maximally. In contrast, our objective is to extract accurate behaviorally relevant signals that closely approximate the ground truth relevant signals as much as possible. To ensure this, we deliberately impose constraints on the model, ensuring that it generates signals that retain neuronal properties while preserving behavioral information to the highest degree possible. Notably, the pivotal operation of striking a balance between the reconstruction and decoding performance of generated signals to extract relevant signals is a distinctive feature absent in other models. At the population level, dimensionality reduction methods aided by task parameters (Kobak et al., 2016; Schneider et al., 2023) are another important way to discover the latent neural embeddings relevant to task parameters, which may provide new insight into neural representations. In contrast with this class of methods, our model focuses on the signal level, not the latent embedding level.

Although we made every effort, our model is still not able to perfectly extract behaviorally relevant neural signals, resulting in a small amount of behavioral information leakage in the residuals. Nevertheless, the signals distilled by our model are reliable, and the minor imperfections do not affect the conclusions drawn from our analysis. In the future, better models can be developed to extract behaviorally relevant signals more accurately, such as incorporating multiple time step information (Pandarinath et al., 2018; Sani et al., 2021; Hurwitz, 2021) and contrastive learning (Schneider et al., 2023) or metric learning (Li et al., 2021) techniques into models.

Studying neural mechanisms through noisy signals is akin to looking at flowers in a fog, which makes it difficult to discern the truth. Thus, removing the interference of irrelevant signals is necessary and beneficial for analyzing neural activity, whether at the single-neuron level or population level.

At the single-neuron level, trial-to-trial neuronal variability poses a significant challenge to identifying the actual neuronal pattern changes. The variability can arise from various sources, including meaningless noise (Faisal et al., 2008), meaningful but behaviorally irrelevant neural processes (Musall et al., 2019), and intrinsic components of neural encoding (Walker et al., 2020). However, it is still unclear to what extent each source contributes to the variability (Faisal et al., 2008). By separating behaviorally relevant and irrelevant parts, we could roughly determine the extent to which these two parts contribute to the variability and explore which type of variability these two parts may contain. Our results demonstrate that behaviorally irrelevant signals are a significant contributor to variability, which may include both meaningless noise and meaningful but behaviorally irrelevant signals as behaviorally irrelevant signals are not pure noise and may carry some structures (thick blue line, Figure 4b and f and Figure 4—figure supplement 1b). Notably, behaviorally relevant signals also exhibit some variability, which may arise from intrinsic components of neural encoding and provide the neural basis for motor learning (Dhawale et al., 2017). Moreover, eliminating the variability caused by irrelevant signals enables us to better observe and compare actual neuronal pattern changes and may facilitate the study of learning mechanisms (Sadtler et al., 2014; Hennig et al., 2021).

At the population level, the dimensionality of neural manifolds quantifies the degrees of freedom required to describe population activity without significant information loss (Lee and Verleysen, 2007; Altan et al., 2021). However, determining the dimensionality of neural manifolds associated with specific behaviors from raw signals is challenging since it is difficult to discern how many variances correspond to irrelevant signals, which often depend heavily on signal quality. A previous study (Altan et al., 2021) demonstrated, through simulation experiments involving different levels of noise, that such noise makes methods overestimate the neural dimensionality. Our results, consistent with theirs, indicate that using raw signals which include many irrelevant signals will cause an overestimation of the neural dimensionality (Figure 4—figure supplement 4). These findings highlight the need to filter out irrelevant signals when estimating the neural dimensionality. Furthermore, this perspective of signal separation has broader implications for other studies. For instance, researchers can isolate neural signals corresponding to different behaviors and explore their shared and exclusive patterns to uncover underlying common and unique mechanisms of different behaviors (Gallego et al., 2018).

At the single-neuron level, previous studies (Carmena et al., 2005; Narayanan et al., 2005) have shown that neuronal ensembles redundantly encode movement behaviors in the motor cortex. However, our results reveal a significantly higher level of redundancy than previously reported. Specifically, prior studies found that the decoding performance steadily declines as neurons drop out, which is consistent with our results drawn from raw signals. In contrast, our results show that decoders maintain high performance on distilled signals even when many neurons drop out. Our findings reinforce the idea that movement behavior is redundantly encoded in the motor cortex and demonstrate that the brain is robust enough to tolerate large-scale neuronal destruction while maintaining brain function (Alstott et al., 2009).

At the population level, previous studies have proposed that motor control is achieved through low-dimensional neural manifolds, with analyses typically using between 6 and 15 PCs (Churchland et al., 2012; Kaufman et al., 2014; Elsayed et al., 2016; Sadtler et al., 2014; Golub et al., 2018; Gallego et al., 2017; Gallego et al., 2020). However, our results challenge this idea by showing that signals composed of smaller variance PCs nonlinearly encode a significant amount of behavioral information, and the number of useful PCs ranges from 30 to 40, far exceeding the usual number analyzed. These results suggest that behavioral information is distributed in a higher-dimensional neural space than previously thought. Interestingly, although smaller variance PC signals nonlinearly encode behavioral information, their behavioral information can be linearly decoded by superimposing them onto larger variance PC signals. This result is consistent with the finding that nonlinear mixed selectivity can yield high-dimensional neural responses and thus allow linear readout of behavioral information by downstream neurons (Rigotti et al., 2013; Fusi et al., 2016). Moreover, we found that smaller variance PC signals can improve precise motor control, such as lower-speed control. Analogously, recent studies have found that smaller variance PCs of hand postures are task-dependent and relate to the precise and complex postures (Yan et al., 2020). These findings suggest that neural signals composed of lower variance PCs may be involved in the regulation of precise motor control.

In the motor cortex, in what form downstream neurons read out behavioral information is still an open question. Previous studies have shown that nonlinear readout is superior to linear readout on raw signals (Naufel et al., 2019; Glaser et al., 2020; Willsey et al., 2022). However, by filtering out the interference of behaviorally irrelevant signals, our study found that accurate decoding performance can be achieved through linear readout, suggesting that the motor cortex may perform linear readout to generate movement behaviors. Similar observations involving raw signals have been reported across various cortices, including the inferotemporal cortex (Majaj et al., 2015), perirhinal cortex (Pagan et al., 2013), and somatosensory cortex (Nogueira et al., 2023). These observations support the hypothesis that linear readout might serve as a general principle in the brain. However, further experiments are needed to verify this hypothesis across a wider range of cortical regions. In motor cortex, different neurons encode behavioral information with varying degrees of nonlinearity, exhibiting complex and heterogeneous response patterns. Despite this complexity of neural encoding, these responses allow for a linear readout of behavioral information. This phenomenon suggests that the complexity of encoding mechanisms may underlie the simplicity of decoding mechanisms.

About studying decoding mechanisms, recent studies (Pitkow et al., 2015; Ganmor et al., 2015; Yang et al., 2021) have focused on investigating how the brain decodes task information in the presence of noise. Unlike previous works, our research specifically explores the decoding mechanisms of behaviorally relevant signals rather than raw signals. We assume that the brain filters out irrelevant signals before decoding the relevant ones. This leads to the question of whether the brain actually adopts this strategy to access relevant signals. Given the existence of behaviorally relevant signals, it is reasonable to assume that the brain has intrinsic mechanisms to differentiate between relevant and irrelevant signals. There is growing evidence suggesting that the brain utilizes various mechanisms, such as attention and specialized filtering, to suppress irrelevant signals and enhance relevant signals (Sreenivasan and Fiete, 2011; Schneider et al., 2018; Nakajima et al., 2019). Therefore, it is plausible that the brain filters before decoding, thereby effectively accessing behaviorally relevant signals. Furthermore, our study reveals that irrelevant signals are the most critical factor affecting accurate and robust decoding, and achieving accurate and robust linear decoding requires weak neural responses. These findings have two important implications for developing accurate and robust BMIs: designing preprocessing filtering algorithms or developing decoding algorithms that include filtering out behaviorally irrelevant signals, and paying attention to the role of weak neural responses in motor control. More generally, our study provides a powerful framework for separating behaviorally relevant and irrelevant signals, which can be applied to other cortical data to uncover more hidden neural mechanisms.

Three datasets with different paradigms are employed, including obstacle avoidance task dataset (Wang et al., 2017), center-out reaching task dataset (Dyer et al., 2017), and self-paced reaching task dataset (O’Doherty, 2017).

The first dataset (dataset A) is the obstacle avoidance dataset. An adult male Rhesus monkey was trained to use the joystick to move the computer cursor to bypass the obstacle and reach the target. Neural data were recorded from the monkey’s upper limb area of the dorsal premotor (PMd) using a 96-electrode Utah array (Blackrock Microsystems Inc, USA). Multi-unit activity (MUA) signals are used in the present study. The corresponding behavioral data (velocity) were simultaneously collected. There are 2 days of data (20140106 and 20140107), and each day contains 171 trials on average. All animal handling procedures were authorized by the Animal Care Committee at Zhejiang University, China, and conducted following the Guide for Care and Use of Laboratory Animals (China Ministry of Health).

The second dataset (dataset B) is publicly available and provided by Kording Lab (Dyer et al., 2017). The monkey was trained to complete two-dimensional eight-direction center-out reaching tasks. We used 2 days of data from subject C (20161007 and 20161011). Each day contains 190 trials on average. Neural data are spike-sorted PMd signals. The behavioral data were simultaneously collected in instantaneous velocity.

The third dataset (dataset C) is publicly available and provided by Sabes Lab (Zenodo dataset) (O’Doherty, 2017). An adult male Rhesus monkey was trained to finish self-paced reaching tasks within an 8-by-8 square grid. There are no inter-trial intervals during the experiment. Neural data were recorded from the monkey’s primary motor cortex (M1) area with a 96-channel silicon microelectrode array. The neural data are the MUA signals. Hand velocity was obtained from the position through a discrete derivative. The recording period for the data (20170124 01) is about 10 min.

For all datasets, the neural signals were binned by a 100 ms sliding window without overlap. As a preprocess, we smoothed the neural signals using a moving average filter with three bins. We excluded some electrode activities with low mean firing rates (<0.5 Hz mean firing rates across all bins) and did not perform any other pre-selection to select neurons. For the computation of the FF, we chose 12 and 14 points as the thresholds of trial length for datasets A and B, respectively; trials with a length less than the threshold were discarded (discard about 7% and 2% trials for datasets A and B), trials longer than the threshold were truncated to threshold length from the starting point. Since dataset C has no trial information, FF is not calculated for this dataset. For the analysis of datasets A and B, we selected 1 day of these two datasets for analysis (20140107 for dataset A and 20161011 for dataset B).

The synthetic dataset is used to demonstrate that d-VAE can extract effective behaviorally relevant signals that are similar to the ground truth signals. The specific process of generating synthetic data is as follows. First, we randomly selected nine larger ${\displaystyle {\mathrm{R}}^2}$ neurons from neurons that ${\displaystyle {\mathrm{R}}^2}$ is greater than 0.1, and three smaller ${\displaystyle {\mathrm{R}}^2}$ neurons from neurons that ${\displaystyle {\mathrm{R}}^2}$ is lower than 0.01 of dataset B (20161011). Second, we used deep neural networks to learn the encoding model between movement kinematics (movement velocity of dataset B) and neural signals using onefold train data. The details of the networks are demonstrated as follows. The networks use two hidden layer multilayer perceptron (MLP) with 500 and 500 hidden units. The activation function is ReLU. A SoftPlus activation function follows the last layer of the networks. The reconstruction loss is the Poisson likelihood function. After learning the encoding model, we used the learned encoding model to generate the ground truth of behaviorally relevant signals from all kinematics data of dataset B. Then, we added white Gaussian noise to the behaviorally relevant signals such that the noisy signals have a signal-to-noise ratio of 7 dB. After adding noise, the ${\displaystyle {\mathrm{R}}^2}$ of the three smaller ${\displaystyle {\mathrm{R}}^2}$ neurons is lower than 0.03. We regarded the noisy signals as raw signals and the added Gaussian noise as behaviorally irrelevant signals. Finally, we separated the synthetic data into five folds for cross-validation model evaluation.

Notation: ${\displaystyle \mathit{x}\in {\mathbb{R}}^n}$ denotes raw neural signals. ${\displaystyle {\mathit{x}}_r\in {\mathbb{R}}^n}$ represents behaviorally relevant signals. ${\displaystyle {\mathit{x}}_i=\mathit{x}-{\mathit{x}}_r\in {\mathbb{R}}^n}$ represents behaviorally irrelevant signals. ${\displaystyle \mathit{z}\in {\mathbb{R}}^d}$ denotes the latent neural representations. ${\displaystyle {\mathit{z}}_{prior}\in {\mathbb{R}}^d}$ denotes the prior latent neural representations. ${\displaystyle \mathit{y}\in {\mathbb{R}}^k}$ represents kinematics. ${\displaystyle f:{\mathbb{R}}^n\to {\mathbb{R}}^d}$ represents the inference model (encoder) of VAE. ${\displaystyle g:{\mathbb{R}}^d\to {\mathbb{R}}^n}$ represents the generative model (decoder) of VAE. ${\displaystyle m:{\mathbb{R}}^k\to {\mathbb{R}}^d}$ represents the mapping from kinematics to prior latent representations. ${\displaystyle h:{\mathbb{R}}^d\to {\mathbb{R}}^k}$ represents an affine mapping from latent representations to kinematics.

d-VAE is a generative model based on VAEs (Kingma, 2013), specially designed to extract behaviorally relevant signals from raw signals. The generative model of d-VAE is

where ${\displaystyle p_m(\mathit{z}|\mathit{y})}$ denotes the conditional prior distribution of latent variables given the kinematics parameterized by feedforward neural networks ${\displaystyle m}$, ${\displaystyle p_g(\mathit{x}|\mathit{z})}$ denotes the conditional prior distribution of raw signals given the latent variables parameterized by feedforward neural networks ${\displaystyle g}$, ${\displaystyle p_{\theta }(\mathit{x},\mathit{z}|\mathit{y})}$ represents the joint distribution of raw signals and latent variables given the kinematics parameterized by parameters ${\displaystyle \theta =(m,g)}$, and ${\displaystyle p_{\theta }(\mathit{x}|\mathit{y})}$ is the marginal distribution of raw signals parameterized by parameters ${\displaystyle \theta }$.

To learn the model, we need to maximize the evidence lower bound (ELBO) of ${\displaystyle p_{\theta }(\mathit{x}|\mathit{y})}$:

where the first term on the right-hand side of Equation 3 is called reconstruction term, the second term is called regularization term, ${\displaystyle q_{\varphi }(\mathit{z}|\mathit{x},\mathit{y})}$ denotes inference model parameterized by parameters ${\displaystyle \varphi }$, and ${\displaystyle D_{KL}(\cdot |\cdot )}$ denotes the Kullback-Leibler (KL) divergence. In d-VAE, we set ${\displaystyle q_{\varphi }(\mathit{z}|\mathit{x},\mathit{y})=q_f(\mathit{z}|\mathit{x})}$, where ${\displaystyle f}$ is the inference model parameterized by feedforward neural networks. Because during the test stage, we cannot obtain the ground truth of kinematics, and we need to use only raw signals to extract relevant signals. Note that d-VAE aims to extract behaviorally relevant signals from raw signals, not generate signals that are too similar to raw signals. Therefore, we modified the objective loss function based on ${\displaystyle {\mathcal{L}}_{ELBO}(\mathit{x})}$ (see Equation 10).

To distill behaviorally relevant neural signals, d-VAE utilizes the trade-off between the decoding and reconstruction abilities of generated behaviorally relevant signals ${\displaystyle {\mathit{x}}_r}$. The basic assumption is generated behaviorally relevant signals that contain behaviorally irrelevant signals harms their decoding ability. Our approach for distilling behaviorally relevant signals consists of three steps, including identifying latent representations ${\displaystyle \mathit{z}}$, generating behaviorally relevant signals ${\displaystyle {\mathit{x}}_r}$, and decoding behaviorally relevant signals ${\displaystyle {\mathit{x}}_r}$.

For each model, we use the fivefold cross-validation manner to assess performance. Specifically, we divide the data into five equal folds. In each experiment, we take one fold for testing and use the rest for training (taking three folds as the training set and one as the validation set). The reported performance is averaged over test sets of five experiments. The validation sets are used to choose hyperparameters based on the averaged performance of fivefold validation data. To avoid overfitting, we apply the early stopping strategy. Specifically, we assess the criteria (loss for training distillation methods, ${\displaystyle {\mathrm{R}}^2}$ for training ANN and KF) on the validation set every epoch in the training process. The model is saved when the model achieves better validation performance than the earlier epochs. If the model cannot increase by 1% of the best performance previously obtained within 10 epochs, we stop the training process.

As previously mentioned, the hyperparameter ${\displaystyle \alpha }$ of d-VAE plays a crucial role in balancing the trade-off between reconstruction and decoding loss. Once the appropriate value of ${\displaystyle \alpha }$ is determined, we can use this value to obtain accurate behaviorally relevant signals for subsequent analysis. To determine the optimal value of ${\displaystyle \alpha }$, we first enumerated different values of ${\displaystyle \alpha }$ to guide d-VAE in distilling the behaviorally relevant signals. Next, we used ANN to evaluate the decoding ${\displaystyle {\mathrm{R}}^2}$ of behaviorally relevant (${\displaystyle {\mathrm{D}}_{\mathrm{r}\mathrm{e}}}$) and irrelevant (${\displaystyle {\mathrm{D}}_{\mathrm{i}\mathrm{r}}}$) signals generated by each ${\displaystyle \alpha }$ value. Finally, we selected the ${\displaystyle \alpha }$ value with the criteria formula ${\displaystyle 0.75\times {\mathrm{D}}_{\mathrm{r}\mathrm{e}}+0.25\times (1-{\mathrm{D}}_{\mathrm{i}\mathrm{r}})}$. The ${\displaystyle \alpha }$ value that obtained the highest criteria score on the validation set was selected as the optimal value.

Note that we did not use neural similarity between behaviorally relevant and raw signals as a criterion for selecting behaviorally relevant signals. This is because determining the threshold for neural similarity is challenging. However, not using similarity as a criterion does not affect the selection of suitable signals because the decoding performance of behaviorally irrelevant signals can indirectly reflect the degree of similarity between the generated behaviorally relevant signals and the raw signals. Specifically, if the generated behaviorally relevant signals are dissimilar to the raw signals, the behaviorally irrelevant signals will contain many useful signals. In other words, when the neural similarity between behaviorally relevant and raw signals is low, the decoding performance of behaviorally irrelevant signals is high. Therefore, the decoding performance of irrelevant signals is a reasonable alternative to the neural similarity.

Regarding the ratio between ${\displaystyle {\mathrm{D}}_{\mathrm{r}\mathrm{e}}}$ and ${\displaystyle {\mathrm{D}}_{\mathrm{i}\mathrm{r}}}$, any ratio greater than or equal to 3:1 is suitable, and we recommend opting for a higher ratio. This recommendation is based on the observation that when the model is biased toward reconstruction (associated with lower ${\displaystyle \alpha }$ values), the decoding performance of relevant signals improves as ${\displaystyle \alpha }$ increases, yet it has not reached saturation. At the same time, the decoding performance of irrelevant signals remains low, but their fluctuations are larger than the improvements in the decoding performance of relevant signals. Consequently, setting the ratio too low poses a risk of selecting an ${\displaystyle \alpha }$ value where both irrelevant and relevant signals demonstrate low decoding performance. This situation fails to satisfy our requirement that relevant signals should exhibit high decoding performance.

For other generative models, we iterate through a range of hyperparameters, generating the corresponding behaviorally relevant neural signals, and subsequently evaluate these signals using ANN. The hyperparameter associated with the signals that exhibit the highest ANN decoding performance is then selected. In other words, the signals corresponding to this particular hyperparameter are chosen as the selected behaviorally relevant neural signals.

All the VAE-based models use the Poisson observation function. The details of different methods are demonstrated as follows:

• d-VAE. The encoder ${\displaystyle f}$ of d-VAE uses two hidden layer MLP with 300 and 100 hidden units. The activation function of the hidden layers is ReLU. The dimensionality of the latent variable is set to 50. The decoder ${\displaystyle g}$ of d-VAE is symmetric with the encoder. The last layer of the decoder is followed by a SoftPlus activation function. The prior networks ${\displaystyle m}$ use one hidden layer MLP with 300 units. The ${\displaystyle \beta }$ is set to 0.001. The ${\displaystyle \alpha }$ is set to 0.3, 0.4, 0.7, and 0.9 for datasets A, B, and C. We perform a grid search for ${\displaystyle \alpha }$ in ${\displaystyle \{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9\}}$, and ${\displaystyle \beta }$ in ${\displaystyle \{0.001,0.01,0.1,1\}}$, and latent variable dimension in ${\displaystyle \{10,20,50\}}$. For the synthetic dataset A and B experiments, the ${\displaystyle \beta }$ and the latent variable dimension are directly set to 0.001 and 50. For synthetic dataset the ${\displaystyle \alpha }$ is set to 0.9, which is grid searched in ${\displaystyle \{0.01\text{--}0.09,0.1\text{--}0.9,1\text{--}9\}}$, with intervals of 0.01, 0.1, and 1, respectively.

• pi-VAE. The original paper utilizes label information (as shown in Equation 6) to approximate the posterior of the latent variable and performs Monte Carlo sampling for decoding during the test stage. However, in our signal generation setting, it is inappropriate to use label information (kinematics) for extracting behaviorally relevant signals during the test stage. As a result, we modified the model to exclude the use of label information in approximating the posterior. The encoder, decoder, and prior networks of our pi-VAE are kept the same as those in d-VAE.

• LFADS. The hidden units of the encoder for the generator’s initial conditions, the controller, the generator are set to 200, 200, 200, and 100 for datasets A, B, and C and the synthetic dataset. The dimensionality of latent factor is set to 50 for all datasets. We perform a grid search for hidden units in ${\displaystyle \{100,200\}}$, and latent factor dimensions in ${\displaystyle \{10,20,50\}}$. The dimensionality of inferred inputs is set to 1. The Poisson likelihood function is used. The training strategy follows the practice of the original paper. For datasets A and B, a trial length of 18 is set. Trials with lengths below the threshold are zero-padded, while trials exceeding the threshold are truncated to the threshold length from their starting point. In dataset A, there are several trials with lengths considerably longer than that of most trials. We found that padding all trials with zeros to reach the maximum length (32) led to poor performance. Consequently, we chose a trial length of 18, effectively encompassing the durations of most trials and leading to the removal of approximately 9% of samples. For dataset B (center-out), the trial lengths are relatively consistent with small variation, and the maximum length across all trials is 18. For dataset C, we set the trial length as 10 because we observed the video of this paradigm and found that the time for completing a single trial was approximately 1 s. The segments are not overlapped.

• TNDM. The hidden units of the encoder for the generator’s initial conditions, the controller, and the generator are set to 64, 64, 100, and 64 for datasets A, B, and C and the synthetic dataset. The dimensionality of relevant latent factors is set to 50, 25, 50, and 50 for datasets A, B, and C and the synthetic dataset. We set the dimensionality of irrelevant latent factors as the same as that of relevant latent factors. The behavior weight is set to 5, 5, 0.2, 0.5 for datasets A, B, and C and the synthetic dataset. We perform a grid search for relevant latent factor dimensionality in ${\displaystyle \{10,25,50\}}$, the hidden units in ${\displaystyle \{64,100,200\}}$, and the behavior weight in ${\displaystyle \{0.2,0.5,1,5,10\}}$. The Poisson likelihood function is used. The other hyperparameter setting and the training strategy follow the practice of the original paper. TNDM uses the same trial length and data as LFADS.

• PSID. PSID uses several (horizon size) past neural data to predict behavior at the current time without using neural observations at the current time. For a fair comparison, we let PSID see current neural observations by shifting the neural data one sample into the past relative to the behavior data. We perform a grid search for the horizon hyperparameter in ${\displaystyle \{2,3,4,5,6,7\}}$. Due to the relevant latent dimension should be lower than the horizon times the dimensionality of behavior variables (two-dimensional velocity in this paper), we just set the relevant latent dimension as the maximum. The horizon number of datasets A, B, C, and synthetic datasets is 7, 6, 6, and 5, respectively. And thus the latent variable dimension of datasets A, B, and C and the synthetic dataset is 14, 12, 12, and 10, respectively.

• ANN. ANN has two hidden layers with 300 and 100 hidden units. The activation function of the hidden layers is ReLU.

• KF. The matrix parameters of observation and state transition process are optimized with the least square algorithm. KF is a linear-Gaussian state-space model designed to provide an optimal estimate of the current state (kinematics in this paper). It does so by considering both the current measurement observations (neural signals in this paper) and the previous state estimate. The KF operates in a recursive and iterative manner, continually updating its state estimate as new observations become available.

We applied PCA to behaviorally relevant and irrelevant signals to get relevant PCs and irrelevant PCs. Then, we used the percentage variance captured (also called alignment index) to quantify how many variances of irrelevant signals can be captured by relevant PCs by projecting irrelevant signals onto the subspace composed of some relevant PCs and vice versa. The percentage of variance captured is

where ${\displaystyle \mathit{D}\in {\mathbb{R}}^{N\times d}}$ is the top ${\displaystyle d}$ PCs of relevant signals (irrelevant signals). ${\displaystyle \mathit{C}\in {\mathbb{R}}^{N\times N}}$ is the covariance matrix of irrelevant signals (relevant signals), and ${\displaystyle \mathrm{Tr}\left(\mathit{C}\right)=\sum _{i=1}^N{\lambda }_i}$, where ${\displaystyle {\lambda }_i}$ is the ith largest eigenvalue for ${\displaystyle \mathit{C}}$. The percentage variance is a quantity between 0% and 100%.

Suppose ${\displaystyle \mathit{x}\in {\mathbb{R}}^{1\times T}}$, ${\displaystyle \mathit{y}\in {\mathbb{R}}^{1\times T}}$, and ${\displaystyle \mathit{z}\in {\mathbb{R}}^{1\times T}}$ are the random variables for a single neuron of behaviorally relevant signals, behaviorally irrelevant, and raw signals, where ${\displaystyle \mathit{z}=\mathit{x}+\mathit{y}}$, ${\displaystyle T}$ denotes the number of samples. The composition of raw signals’ variance is as follows:

where ${\displaystyle \mathbb{E}}$, ${\displaystyle \mathrm{V}\mathrm{a}\mathrm{r}}$, and ${\displaystyle \mathrm{C}\mathrm{o}\mathrm{v}}$ denote expectation, variance, and covariance, respectively. Thus, the variance of raw signals is composed of the variance of relevant signals ${\displaystyle \mathrm{V}\mathrm{a}\mathrm{r}(\mathit{x})}$, the variance of irrelevant signals ${\displaystyle \mathrm{V}\mathrm{a}\mathrm{r}(\mathit{y})}$, and the correlation between relevant and irrelevant signals ${\displaystyle 2\mathrm{C}\mathrm{o}\mathrm{v}(\mathit{x},\mathit{y})}$. If there are ${\displaystyle N}$ neurons, calculate the composition of each neuron separately and then add it up to get the total composition of raw signals.

The correlation matrix is reordered with a simple group method. The order of reordered neurons is determined from raw neural signals, which is then used for behaviorally relevant signals.

The steps of the group method are as follows:

• Step 1: We get the correlation matrix in original neuron order and set a list A that contains all neuron numbers and an empty list B.

• Step 2: We select the neuron with the row number of the largest value of the correlation matrix except for the diagonal line and choose nine neurons in list A that are most similar to the selected neuron. We selected the value nine because it offers a good visualization of neuron clusters.

• Step 3: Remove these selected neurons from list A, and add these selected neurons in descending order of correlation value in list B.

• Step 4: Repeat steps 2 and 3 until list A is empty.

The code is available at https://github.com/eric0li/d-VAE, copy archived at eric0li, 2024.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We thank Yiwen Wang for sharing the monkey obstacle avoidance reaching data. We thank Yuxiao Yang and Huaqin Sun for valuable discussions. This work was supported in part by the National Natural Science Foundation of China under Grants 62336007, in part by the Key R&D Program of Zhejiang under Grant 2022C03011, in part by the Starry Night Science Fund of Zhejiang University Shanghai Institute for Advanced Study under Grant SN-ZJU-SIAS-002, and in part by the Fundamental Research Funds for the Central Universities.

All animal handling procedures were authorized by the Animal Care Committee at Zhejiang University, China (No. ZJU20220142), and conducted following the Guide for Care and Use of Laboratory Animals (China Ministry of Health).

1. Sent for peer review: May 8, 2023
2. Preprint posted: May 9, 2023
3. Reviewed Preprint version 1: July 14, 2023
4. Reviewed Preprint version 2: February 2, 2024
5. Reviewed Preprint version 3: July 18, 2024
6. Version of Record published: August 9, 2024

You can cite all versions using the DOI https://doi.org/10.7554/eLife.87881. This DOI represents all versions, and will always resolve to the latest one.

© 2023, Li et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.