A Unified, Scalable Framework for Neural Population Decoding
Authors: Mehdi Azabou
Vinam Arora
Venkataramana Ganesh
Ximeng Mao
Santosh B Nachimuthu
Michael Jacob Mendelson
Blake Aaron Richards
Matthew G Perich
Guillaume Lajoie
Eva L Dyer
Venue: NeurIPS 2023 poster (2023)
Keywords: neural population
brain decoder
transformer
tokenization
sequence-to-sequence
electrophysiology
brain-computer interfaces
TL;DR: This paper introduces a scalable and unified learning framework for efficient, large-scale neural decoding across diverse multi-lab multi-animal multi-session neural recordings.

Abstract: Our ability to use deep learning approaches to decipher neural activity would likely benefit from greater scale, in terms of both the model size and the datasets. However, the integration of many neural recordings into one unified model is challenging, as each recording contains the activity of different neurons from different individual animals. In this paper, we introduce a training framework and architecture designed to model the population dynamics of neural activity across diverse, large-scale neural recordings. Our method first tokenizes individual spikes within the dataset to build an efficient representation of neural events that captures the fine temporal structure of neural activity. We then employ cross-attention and a PerceiverIO backbone to further construct a latent tokenization of neural population activities. Utilizing this architecture and training framework, we construct a large-scale multi-session model trained on large datasets from seven nonhuman primates, spanning over 158 different sessions of recording from over 27,373 neural units and over 100 hours of recordings. In a number of different tasks, we demonstrate that our pretrained model can be rapidly adapted to new, unseen sessions with unspecified neuron correspondence, enabling few-shot performance with minimal labels. This work presents a powerful new approach for building deep learning tools to analyze neural data and stakes out a clear path to training at scale for neural decoding models.

A Unified, Scalable Framework for Neural Population Decoding
Mehdi Azabou':s Vinam Arora', Venkataramana Ganesh's Ximeng Mao
Santosh Nachimuthu', Michael J. Mendelson', Blake Richards
Matthew G. Perich Guillaume Lajoic Eva L. Dyer
Georgia Tech, 2 Mila, 3 Université de Montréal, McGill University
Abstract
Our ability to use deep learning approaches to decipher neural activity would likely benefit from greater scale, in terms of both model size and datasets. However, the integration many neural recordings into one unified model is challenging, as each recording contains the activity of different neurons from different individual animals. In this paper, we introduce a training framework and architecture designed model the population dynamics of neural activity across diverse, large-scale neural recordings. Our method first tokenizes individual spikes within the dataset to build efficient representation of neural events that captures the fine temporal structure of neural activity. We then employ cross-attention and a PerceiverIO backbone further construct latent tokenization of neural population activities. Utilizing this architecture and training framework, we construct a large-scale multi- session model trained on large datasets from seven nonhuman primates, spanning over 158 different sessions of recording from over 27,373 neural units and over 100 hours of recordings. In a number of different tasks, we demonstrate that our pretrained model can be rapidly adapted to new, unseen sessions with unspecified neuron correspondence, enabling few-shot performance with minimal labels. This work presents a powerful new approach for building deep learning tools to analyze neural data and stakes out a clear path to training at scale.
Introduction
Recent advances machine learning, particularly in the context large-scale pretrained models like GPT [1, 2, 3, 4], have showcased the immense potential of scaling up, both the terms of the size of datasets and models [5, 6]. Similarly, in neuroscience, there is growing need for foundational model that can bridge the gaps between diverse datasets, experiments, and individuals, allowing for more holistic understanding of brain function and information processing [7]. The development of such model would allow researchers uncover underlying patterns and interactions within neural populations [8]. and potentially allow for more robust decoding of brain states.
However, creating a large-scale neural decoding model that can effectively combine spiking datasets from various sources is complex challenge [9]. One of the central challenges is the lack of shared "vocabulary" in neural recordings. Unlike the case for text text-wherein every document written in a given language shares basic lexicon for tokenization- between neurons in different individuals. As such, every recording from a different individual involves a unique set of neurons that cannot be easily aligned with another set. An additional core challenge lies in the inherent variability of neuron sets observed across different days [10]. Even when monitoring the same individual, variability in electrode/tissue interface can lead to distinct
no one-to-one correspondence
CCntacct: [mazabou,evadyer/ @ @gatech.edu. Project page and code: httpss://pyye brain, 1o
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

neuron sets which, despite advanced sorting methods, can lead to inconsistencies in input channels across sessions [11, 12, 9]. Overall, this lack of correspondence across recordings complicates the integration of information from different experiments and individuals, ultimately hampering efforts to construct a unified perspective on population-leval interactions and dynamics in the brain.
In response to these challenges, we propose a new framework for large-scale training on neural spiking data called POYO (Pre-training On man' neurOns). This framework is designed to enable scalable and efficient training across multiple sessions of neural recordings, even when spanning different sets of neurons with no known correspondence. Our approach centers around a novel tokenization scheme that transforms individual neural action potentials, or "spikes". into discrete tokens, thereby preserving the neural code's finest temporal structure while simultaneously enhancing computational efficiency. The resulting tokenization not only allows for a more effective representation of neural activity but also paves the way for training on larger volumes of data. We combine this input tokenization method with an architecture that builds on the PerceiverlO [13] to compress the input spikes into a latent space and learns interactions across spikes in time and across neurons.
We evaluate the performance of our proposed approach on data from over 158 sessions from open electrophysiology datasets from seven non-human primates (NHPs), spanning over 27,373 units and 100 hours of recordings. We demonstrate that through pretraining on large amounts of data, we can transfer with very few samples (few-shot learning) and thus improve overall brain decoding performance. Our work not only presents an innovative framework for training large models on neuroscience datasets, but also offers insights into the scaling laws that govern decoding from neural populations. By enabling the development of large pretrained models for neural decoding, our approach advances the field of brain-machine interfaces and other decoding applications.
The main contributions of this work include:
A framework for large-scale training on neural recordings: We present a novel framework for training transformer models end-to-end on multi-session and across-individual elec- trophysiology datasets derived from neural populations, enabling efficient and effective decoding from diverse range of neural recordings.
Innovative spike-based tokenization strategies: We introduce a fundamentally different way tokenize neural population activity. Our approach tokenizes individual spikes (events) across neural populations, preserving fine temporal structure and enhancing computational efficiency by adopting a sparse representation of the data.
Pre-trained models for neural decoding: We build two large pretrained models (POYO-1, POYO-mp) that can be fine-tuned on new sessions and across recordings from different animals and new behavioral tasks. We will make the weights and code available, and provide both pretrained models as a resource to the community.
2 Approach
The transformer architecture [14], originally introduced in the context of natural language processing (NLP), has shown remarkable flexibility and effectiveness in various domains, especially in the presence of large and diverse datasets [15, 3]. In this work, we explore how to leverage this versatility in the neural data domain.
2.1 Tokenizing neural population activity
Neurons communicate asynchronously using electrical impulses called spikes. The timing and frequency of spikes encode signals that convey information about the external world and coordinate the internal dialogue within the brain. In many neuroscience experiments, neural activity is recorded from the brain through multiple electrodes, and then processed [16] to extract the spiking events for a set of neural "units" [17].3 The resulting data are multi- aurti-variate event-baser sequences that are typically very sparse relative to the number of recorded time points.
'Poyo is the exclamation used by Kirby, who has insatiable appetite. Similarly, POYO consumes spikes from many neural datasets and combines them into one unitied model.
WWe use the term units instead of neurons because recordings may include both single-units (isolated neurons) and multi-units (contain spikes from multiple neurons).

---

unit ids
Behavior Queries
Raw Spikes
unit
dataset
sample
samplc
simple
Batch
POYO
Learned Latents (7v)
Predicted Behavior
Self Attention
Cross Attention
Zr
Cross Attention
Linear
Vs.
x(2
Figure I: Overview our approach. Input spike tokens are compressed into smaller latent tokens which are processed through multiple self-attention blocks, finally time-van varying outputs are predicted by querying the latent space. All tokens in this model are assigned timestamp, which are used towards rotary position encoding.
Neurons and their spikes are of course not independent of one another. Rather, they are part of complex, interwoven tapestry of neural activity, with each spike contributing to a collective dialogue across billions of neurons. The true challenge and opportunity lie in interpreting these spikes not in isolation, but in the context of this broader "conversation"
When trying to scale up and train on many sessions of data, we are faced with the challenge of having recordings of neural conversations with completely different of "speakers" for which don t know the identity or their functional tuning (or what they respond to). This is because each time we record from the brain, we are tuning into a conversation between new set of speakers. However, as with language, there is some reason to think that neurons are ultimately communicating on similar topics, e.g. sensations from the external world, internal states of the body, muscle commands, In other words, the lexicon may be different, but everyone is talking about similar subjects. Despite the variability of our data, our aim to decipher these "conversations" neural datasets without fixed or known correspondence across their inputs.
way that generalizes to new
The neural tokenizer. Based upon these motivations, we propose novel approach for tokenizing neural population activity where each spike is represented as token (see Figure 1). In this case, each token can be defined in terms of which unit came from (via learnable embedding) and the time that the spike event was detected. By representing our data in this way, we never define or fix the expected number units, and the model can ingest populations of arbitrary size, and thus can be trained across many datasets. At the same time, this approach also avoids having specify specific time resolution for the neural code, as is the case with binning, and gets rid of the accumulation of lot of sparse tokens containing events.
More concretely, we assign each unit a unique identifier and a corresponding D-dimensional learnable embedding. Let Unit denote a lookup table that associates each unit to its unit embedding. Akin to word embeddings that encode semantic meaning and relationship between words, the unit embedding space encodes something about the "speaker's ss identity and role in neural computation. unit fires sequence of spikes in a context window of time [0, T]. Each spike will be represented by token characterized by (x;, (2), where
X; UnidEmbbed(spike i's unit)
the learned embedding associated with the unit that emitted the spike and t; is the event timestamp.
Collectively, and for an arbitrary set of units, we combine all the spikes into sequence of length M. The neural population activity represented by [xi,...,X xM] and their corresponding times [/,,.., l.). While the size of our context window T stays fixed, M will vary depending on the number of units the population and their firing rates (see Figure 1). Note that all spikes from specific unit have the same unit embedding, and only differ in their timing.

---

2.2 Building the latent space
Compressing the input sequence. Rather than processing the input sequences with selffattention, which is quadratic in sequence length and can be expensive for long sequences, we use a perceiver encoder module [18] that summarizes the input sequence into a shortened "latent" sequence. Let's
consider a sequence of M input tokens X = Xx,..... XM]
and sequence of learned latent tokens Zo [2..1.... N. where 70, E TR. and N < M. We use cross-attention to pull information from the input sequence into a compressed latent sequence. The queries Q = W Zo are a projection of the learned latent tokens Zo and the keys and values are projection of the input tokens: K
q
W and V = X, respectively, The cross-attention operation can be expressed as follows:
Zo
QK
2 f- Cross-Attn
Q
K
softmax
V.
dd
We use the standard transformer block with pre-norman pre-normalization layers and feed-forward nets.
Attention in the latent space. Following the compression of the input sequence through cross- attention, we apply multiple self-attention blocks on the latent token sequence, which now costs O(N2) instead of (M²) with N < M. Let Zi denote the embedding of the latent tokens at the /th layer. The self-attention operation can be expressed as follows:
Zr
QK
Z++1 + Self-Attn( Q.K.
K
V = softmax
V
(2)
Vdg
We denote the final latent sequence obtained after L layers as Z [ZL,1,.... ZL.N].
2.3 Encoding relative timing information
To incorporate timing information, we leverage rotary position encoding (RoPE) [19] across all attention layers in the architecture. Unlike traditional position encoding methods that inject absolute position information into the token's embedding, RoPE applies a rotation operation in the query, key embedding space, allowing each token in a sequence to attend to others based on its relative positional information.
Recall that each input token has an associated timestamp t4, We will also assign a timestamp to each of the latent tokens: we divide the latent tokens into groups of equal size and spread each group uniformly over the context window [0, 7]- This method allows us to capture temporal relationships between the latent tokens and the input tokens, enabling temporal understanding of the encoded sequence. By distributing the latent tokens evenly within the context window [0, T), we create a structured temporal representation, which preserves and propagates temporal information all the way through the model. A detailed formulation of the relative encoding can be found in Appendix A.
2.4 Querying the latent space
Having built latent space which can encode and model any population of units, we now want flexible way to readout behavioral variables. The output of the latent encoder is the latent sequence of size N, while the desired output can be any arbitrary sequence of length P. Let us consider the task of hand velocity decoding for example, in the context window [0, T). the length of the output sequence datasets sourced from various
will depend on the sampling frequency. Since aspire to train labs, the sampling rate can differ significantly. We thus need a flexible mechanism for predicting outputs of varying lengths and querying from the neural activity at specific points in time.
We define sequence of output tokens Y, .yo. Po, where the number of output time points, which can change across sequences. Each output token is defined by (yo,s, tqut). Initially all the output tokens within a session are set to the same learned embedding Yo,4 Yo E g Vi. The output Y is obtained by querying the latent space through cross-attention:
Yo
Zr
QK
Y +- Cross-Attm Q
K.
V
= softmax
V
(3)

---

Since we use rotatary embeddings, and both the latent and output tokens have assigned timestamps, the querying mechanism can leverage the relative position of latent and output tokens to extract the temporally relevant context that enable prediction of the behavioral variable of interest.
Session embeddings. To account for variability of the experimental setups in real world settings, propose to define a learnable session embedding which captures the hidden experimental variables. This information is injected into the output query yo. where again use a lookup table register each new session we encounter. To illustrate how we probe the model to generate an output, we can think about it as asking a question of the form:
"Predict [BEHAVIOR] at time [TIMESTAMP] under the experimental conditions of [SESSION IDJ".
This produces a of output tokens of the dimension of the number of queries which we then pass through an MLP to generate the behavioral variables of the desired dimension (e.g., 2D for hand velocities in planar movement task).
2.5 Unit identification in new sessions
Our design of the unit embedding space allows model to learn latent information about the units encounters, as well as capture the relationship between units in the population. Given a new recording with unidentified units, we transfer model by mapping these units into the unit embedding space. To do this, we introduce an approach that we call "unit identification". which leverages gradient descent learn the embeddings of new units. In this approach, freeze all existing weights of the model and simply add new rows to the UnitEmbed(()) lookup table for each of the new units, well as a new session embedding. Notably, the bulk of our model which maps the neural population activity to behavior is unchanged and is simply transferred to the new dataset. In our experiments, we find that this approach is surprisingly effective and allows us rapidly integrate new datasets into the same underlying model. Further details and analysis on the robustness of this approach are presented in Appendix C.2.
Experiments
In this section, we demonstrate the promise of our approach for large-scale training and examine the benefit of scaling in neural population decoding.
3.1 Datasets and experiment setup
One of the key advantages of our approach is its ability to scale to handle large amounts of neural data, including sessions from different numbers of neurons, across different tasks and recording setups, and from different animals. Thus we set out to build a diverse dataset large enough to test our approach. We curated a multi-lab dataset with electrophysiological recordings from motor cortical regions, where neural population activity has been extensively studied [10]. and deep learning tools and benchmarks have recently been established [27]. In total, aggregated 178 sessions worth of data, spanning 29,453 units from the primary motor (MI), premotor (PMd). and primary somatosensory (S1) regions in the cortex of nonhuman primates (see Table 1). We place this in the context of standard analyses within single lab or paper which typically involve 10's of sessions and a few hundred neurons.
All of these neural recordings were collected while the animals performed various motor tasks that vary in their inherent complexity (see Figure 2A-B). The center-out (CO) task is relatively stereotyped, with the animal making a reach to one of eight targets after receiving a go cue, and
Study
Regions MI, PMd
# Indiv
Sess
Units
# In
# Out
Tasks
Perich et al. [20, 21, 22, 23]
4
I17
11,557
143M
20M
CO. RT
Churchland ct al. [24]
M1
2
9
1,728
706M
87M
CO
Makin et al. [25]
M1,
2
47
14,899
123M
15M
RT
Flint et al. [26]
M1
1
5
957
7.9M
0.3M
CO
NLB- LLB-Maze [27]
MI
1
1
182
3.6M
6.8M
Maze
NLB RTT [27]
MI
1
130
1.5M
2.8M
RT
Table 1: Datasets used in this work. CO:
Center- RT; Random Target.
5

---

A
Munkey
UN/2/2014
Monkey M
192220015
C
Pardhotsl Morkopuc Mii Tade
%
Moskey
Manbory 2s
0492//016
Center-out
Controllen: Manipulandium
P010 Scoket
Ri 1.0 ceepae- mose)
50
92
24
B
Monkwy C
1220/2013
Monkey M
01/7E/2314
feicheral MarkeguC
Dank
09
Monkey
Nankay
BOG3--2144
Random target
Controller: Manipulandum
50
appl Socahie
100 (iiggl Spocion mode
Figure 2: Building multi-session model spanning multiple animals and tasks. (A) Center-out reach and (B) Random target task [20], along with examples of true and predicted behavior (x.y velocities). (C-D), we show the decoding performance for single-session models (gray) and the POYO-m, multi-session model (green).
then returning to the center. In contrast, the random target (RT) task significantly more complex. The animals make continuous and self-paced movements with new targets appearing in succession random locations the workspace. In addition the greater exploration of the workspace and heterogeneity in movements, this task allows individual to plan their next movement while finishing the execution of the current movement leading to greater complexity in neural dynamics. Other sources of variability that we find across labs include the: choice of pre-prowe pre-processing algorithms (spike sorting or threshold crossing analysis), type controller (manipulandum, touch screen), and sampling rate when recording behavior (100Hz to 1kHz). We do not re- e-prrcesss the data or attempt to standardize it across labs or tasks. Further details on the datasets are provided Appendix
Experiment setup. Throughout all of our experiments, we use a context window of 1s and do not segment data into trials during training. We only use the trial structure when reporting the decoding performance, in particular, center-out sessions are evaluated during the reaching movement. We train the model with N 512 latent tokens and a dimension D 128. We use the LAMB optimizer [28], and employ cosine decay of the learning rate at the end of training. For every session, we holdout 20% of the trials for testing, and 10% for validation. We use 1-GPU and 8-GPU setups for single-session and multi-session models, respectively. All details can be found in Appendix A.
3.2 Testing the model single sessions
To first investigate the performance of our architecture when trained on single session, we trained single-session models on 100 different recording sessions acquired from three nonhuman primates (Monkey C, Monkey M, Monkey Ja), each containing anywhere from to 106 minutes data [29]. In all of these recordings, the same behavioral task setup, behavioral recording apparatus (10 ms temporal resolution), and spike sorting procedure was used.
Across all 100 single-session models, we obtained an average R2 of 0.9347 on CO and 0.8402 for RT sessions. When we compared these single-session results with existing models for neural decoding. including a Wiener Filter, GRU and MLP [30], we found that our approach consistently outperforms these baselines, with even greater improvements observed on the RT task. Additionally, we found that our single-session session models are stable for wide range of hyperparameters (Appendix C.3).
3.3 POYO-mp: Building a large pretrained across-animal, multi-session model
To investigate the question of how training with more sessions of data can improve brain decoding, we trained a large model (POYO-m) 24 layers) on all 100 sessions that we studied in our single-session analysis. total, we used 9,789 units and 4,367 neuron-hours (number of neurons amount of time recorded) to train our POYO model.

---

Same animal, New day
New animal
Method
Monkey C- CO (2)
Monkey T CO (6)
Monkey T RT (6)
From scratch
Wiener Filter
0.8860 + 0.0149
0.6387 + 0.0283
0.5922 + 0.0901
GRU
0,9308 + 0.0257
0.8041 + 0.0232
0.6967 - 0.1011
MLP
0,9498 + 0.0119
0.8577 + 0.0242
0.7467 = 0.0771
P0Y0-JSingle-session] POYO-mp Unit ID
0.9682 + 0.0111
0,9194 + 0.0185
0.7800 = 0.0702
0,9675 + 0.0079
0,9012 i 0.0271
0.7759 Σ 0.0471
Pre- trained
POYO-mp mp Finetune
0.9708 + 0.01
0.9379 + 0.0193
0.8105 + 0.0561
POYO- Unit ID
0,9677 + 0.0096
0,9028 0.0248
0.7788 + 0.0548
PΟYO-1 1 Finetune
0.9683 + 0.0118
0.9364 + 0.0132
0.8145 + 0.0496
Table 2: Behavioral decoding results across neural recordings from two nonhuman primates performing two different tasks. All the baselines and the single-session model are trained from scratch, while POYO-mp and POYO-1 are pretrained. The standard deviation is reported over the sessions. The number of sessions in each dataset is contained in (.) and the top performing models in each category are indicated boldface (first) and underlined (second).
This model achieves an average test R2 of 0.9512 on the center-out tasks and 0.8738 on the random target tasks, which is marked improvement over the single-session average. When we compared the performance of our multi-session model to single-session models head-to-head (see Figure 2C-D), we found that across the board, multi-session training improves over single-sessions, and we observed even greater improvements for datasets with fewer trials (indicated by the size of the circles in the plot). On RT sessions, we observe an even bigger improvement over the single session models. These results suggest that there are significant benefits in woint-training across multiple sessions, especially when decoding more complex behaviors like the random target task.
Scaling analysis. By enabling multi- session training, we can start to ask ques- tions about the extent to which having more data or more parameters can improve de-
Center- Out
Random Target
Random Target
35.5
68
86
24 lysrs
Average R!x 100
Swiing model
More parameters
95.0
67
87
coding (Figure 3). Thus, we studied the improvements obtained for three different depths of models, with L=6 (3,8M), L=14
94.5
06
86
AT
Training
Add
(7.4M), and L=26 (13M) layers (param- eters). For both CO and RT tasks, we find that even at the same depth (L=6), our 1 multi-session model shows an im- provement over the average performance of 100 singl-ssesson models. We see further improvements with increasing the model depth, with the RT task benefiting even fur- ther from scaling. When we fix the model depth I==,, and study how training jointly with more data contributes to improvement in performance, we find an improvement in performance when comparing single-session training (average 1.2M tokens) training jointly on all RT sessions (26M tokens). When we also train jointly with other CO sessions (134M tokens in total), and increase model depth to help in accommodating the more growing and diverse training set, we continue seeing improvement in the decoding performance on RT sessions, accumulating to 4% over our single-sessin baseline.
Nult-Gasa on Single sesdin
Vult Geelnn Singe Sesdion
Singlo seision
93.0
83
83
3.8M
7.0M
13M
3.8M
7 AM
13M
1.
1.2M
26M
134M
Num parametors
Num takens
Figure 3: Scaling curves. R is evaluated and averaged over the test splits on both CO and RT tasks. The single-session performance the mean over 100 different models.
3.4 Transferring to new sessions
After pretraining, we can then test on new sessions with unknown neurons using either the (i) unit identification approach we described in Section 2.5, or (ii) full firetuning of the model weights. Recall that when we use unit identification approach, we freeze the weights of the model and only learn new unit and session embeddings.
Results on held out sessions from the same animals. We first tested the unit identification approach on held-out sessions from Monkey C that are unseen during training (Table 2, Left). In this case, we don't 't have correspondence between units in the training and the testing conditions. Surprisingly, we find that we can achieve comparable performance with unit identification on the pretrained model (0.9675) with that of the single-session models trained fully from scratch (0.9682). With further finetuning of the rest of the model's weights, we can improve the accuracy further (0.9708) to go beyond the accuracy of the single-session models. This highlights the robustness
7

---

Sample Efficiency
Compute Efficency
Sample Efficleccy
Compute Emficencc
1.0
1.00
0.8
0.75
0.8
0.50
steps
2 0.6
-.00 steps
0.7
0.6
Sne Serdin made
0.5
Single Searin mode
Led Mertifaaai
0.5
0.4
16
32
64
128 167
1000
2000
3000
8
16
52
64
86
1000
2000
5000 4000
Number of sample triats
Number of training steps
Number sample trials
Number of training stops
(a) Monkey C, heldout session
[b) Monkey T (new subject), RT session
Figure 4: Sample and compute efficiency for training, unit identification, and fine-tuning approaches. On the left, we show the sample and compute efficiency for heldout CO session. On the right. we plot the sample and compute efficiency for a new animal not seen during training.
of our model and its flexibility to accommodate fresh data with even a simple input mapping (unit identification).
We compare with the single session model and a GRU baseline (Figure 4A) with our multi-session model trained on the same number of trials. Both our single-session and finetuning approach achieve good performance scaling as we increase the number of samples, with the finetuning maintaining a gap over the single session models over the entire sampling space.
Results on new animal performing the same tasks, Next we tested our model on 12 sessions from a completely new animal (6 CO sessions, 6 RT sessions from Monkey T) that was not included in the training of the model (Table 2, Right). When we applied our unit identification approach on the multi-session model (POYO-mp mp + Unit ID), we see bit of a dip in overall accuracy from the single-session model (POYo-[Single-een -eesgion-sssion):: however, when we fine-tune the weights of the model further, we achieve an accuracy of 0.9379 on CO, which is a significant improvement over all of the other baselines. For the RT task, the single session model is 0. 0.7569 and with unit identification we achieve 0.7669 and with full fine-tuning we get up to 0.7916 (Figure 4B). These results are promising and show that we can use our pretrained model on new animals with only a few minutes of labeled data.
Performance on the Neural Latents Bench- mark (NLB). To understand how well our pre-trained model performs on data collected from new animals performing novel tasks with different equipment (example: touch screen vs. manipulandum), we applied our pretrained model to the MC-Maze (Monkey L) and MC-RIT (Monkey I) datasets from the NLB (Table [)]] The NLB serves as a benchmark for neural representation learning and decoding and thus we can include other single-session baselines to comparisons, including self supervuee models AutoLFADS [31] and NDT [32] which produce denoised firing rate estimates over which we fit a linear layer (+ Linear), a supervised variant of NDT [33] (NDT-Sup), and EIT [33]. Both datasets contain single sessions from new animals performing movement tasks that we haven't 't seen during training.
Method
NLB-Mat L-Ba
NLB-RTT
Wiener Filter
0.7485
0.5438
GRU
0.8887
0.5951
From scratch
MLP
0.8794
0.6953
AutoLFADS Linear [31]
0.9062
0,5931
NDT + Linear [32]
0.8929
0,5895
NDT-Sup [33]
0.8708
0.4621
[33]
0.8791
0.4691
POYO-ISingle-session]
0.9470
0.6850
Pretrained
POYO-mp Unit ID
0.8962
0.7107
POYO-mp Finetune
0.9466
0.7318
POYO-1 + Unit ID
0.9329
0.7294
POYO-1 + Finetune
0.9482
0.7378
Table 3: Behavioral decoding results for datasets from the Neural Latents Benchmark [27]. Best performing model is in bold and second best model is underlined.
On the NLB-Maze dataset, we obtain a R2 of 0.8952 after unit identification, which is competitive with the baselines. These results are surprising since we do not modify the model's weights, and yet our pretrained model yields competitive results on a dataset collected under very different conditions. When finetune the model, we boost the performance even further establishing a 4.4% gap over the best baseline. Similar trends can be observed for the RTT task (Monkey I), with even larger (2%) improvement after finetuning.
3.5
POYO-1: A multi-lab, multi-task model for neural decoding
Given the impressive transfer results of POYO-mp to datasets from different labs, we ask whether we can use our approach to build model that spans even more diverse recording setups that we expect
8

---

A
Monleay
04711/2016
Monkey
C.i
Mankeju V.J Tark
Mankepi Task RT
Moniage LL Taok
90
100
Marshey
Narkey
014807517
80
Random target Controller: screen
Singe
Dhygk
B
Martory
Morkey
C.ii
Chuchland stal
Monkapu
Loannd Sesien Emberddinn FOHO
Monieyy
Nankay
Maze
Contralior: Touch
X00
Figure 5: Scaling up to more tasks and diverse neural and behavioral recording conditions. (A) Random target task from Makin et al. included in training of POYO-1 and (B) Maze task from NLB heldout for transfer testing. In (C) decoding accuracy for the single-session (gray) and the POYO-1 model (blue). (D) PCA projection of the learned session embeddings.
to encounter when trying to unify data from many sources. In total, we used datasets from seven non- human primates spanning three different labs, with total of 27,373 units and 16,473 neuron-hour for training our model. We call this pretrained multi-lab, multi-task model POYO-1
Even in light of the high amounts of variability across these different datasets, POYO-1 provides consistent improvements over the single-session models (Figure 5C). When tested on number of different transfer tasks (Tables 2, 3), we again find that the unit identification and finetuning approaches provide effective strategies for adapting our pretrained model new datasets. We obtain notable performance on the NLB-Maze, where we find that we obtain R of 0.9329 with only unit identification remapping. an almost 4% improvement over the unit identification result for our POYD-mp pretrained model, suggesting that we cover more of the task space with POY0-1
When comparing the POYO-1 model with PΟΥO-mp model, it is clear that both methods have their strengths. PΟYO-mp excels on the datasets sourced from the same lab, with 0.8788 on RT sessions compared to POYO-1's 0.8664. On the other hand, both models exhibit great transfer capabilities, with POYO-1 having the edge especially when using unit identification, indicating its ability to generalize better across diverse experimental conditions. This flexibility and scalability make these methods promising tools for future research in analyzing neural data from diverse sources.
Visualizing the session embedding space. We visualized the learned task embeddings see if the model learns relationships between sessions seen during training (Figure 5D), even though it was not explicitly trained with this information. This analysis revealed clusters corresponding to different data sources and tasks, suggesting that our model has not only learned to identify patterns within each session but also recognized overarching structures across sessions. In particular, we find that the datasets collected in the Miller Lab (C,M.J) used to train the POYO-mp model are mapped into similar regions of the session latent space, and (I and L) sessions are mapped to their own distinct clusters.
4 Related Work
Transformer architectures for time-series data. Transformers have emerged powerful model for processing time-series data due to their ability to capture long oongerange dependencies [34, 35]. The key to their success lies in the self-attention mechanism [14], which allows the model weigh the importance of each time step in the input sequence when producing an output. A typical approach applying transformers to time-series data involves discretizing the continuous time domain into discrete bins [35], and treating the binned data as sequence of tokens. Each bin is then linearly projected into high-dimensional embedding, which is input to the transformer model. The Perceiver framework [18, 13] moves away from the idea of patches, by directly processing the input bytes,

---

using cross-attention layers to reduce the complexity of attention. In addition to regular timeseries, there has been interest in applying transformers to model irregular timeseries including event stream and point process data [36, 38]. However, many of these models work on univariate event streams and when they extend to multivariate cases, they assume fixed and known input channels.
Transformers applied to neural data. With the recent advances in transformers for sequence modeling in many different time-series, transformers have also recently found successful applications in building representations of neural population activity [32, 33, 39]. In these models, the spiking activity is first binned and then the estimated bin counts are tokenized. In the neural data transformer (NDT) model [32]. the firing rates of all neurons in the population are embedded jointly in one token (time step or bin). In the embedded interaction transformer (EIT) [40]. neurons are considered independently in one stage of processing and at population level in second stage. and thus the whole dataset is tokenized over both neurons and time. the spatiotemporal (STNDT) model [39], two different tokenizations are also considered, one in space and one in time, and two representations are learned jointly for both tokenizations. In all cases, binned data are used and the models are trained on single session and fixed set of neurons.
Multi-session training and alignment. The idea of decoding across multiple sessions has been explored in previous work [41, 42, 43, 44]. In many of these works, an initial baseline representation formed on one day and alignment- based approaches are used transfer a model trained on one session across recording days [23, 45, 12, 46, 47J. A subset of these methods [41, 43] can be trained on many sessions jointly. but rely on the assumption of shared dynamies or structure of single task to achieve alignment. To the best of our knowledge, our work is the first to demonstrate multi-session transfer across subjects performing different tasks, and the first to demonstrate scaling across different data and model sizes.
Discussion
In this paper, we introduce a novel framework for training transformers on large multi-session, multi-task neural activity datasets. To tackle the challenges of training such large heterogeneous sources, we introduce a novel spike-level tokemization scheme and architecture that enables the model to learn from populations with varying numbers of neurons. We show that training single unified model on multiple recordings is possible, and find that it leads to improved decoding performance. Finally, we build two large pretrained models (POYO-1, POYO-mp) that be efficiently fine-tuned on new datasets, and make them available as a resource to the community.
In contrast to models trained on single dataset, the pretrained models that we have developed provide a potential way to compare and contrast datasets, and also understand common motifs of activity and dynamics that may be shared across different sessions, tasks, and individuals. Thus, it will be critical to develop tools to probe the patterns and motifs learned by such models and characterize the neural mechanisms underlying different tasks and computations. In particular, we look to understand how spike tokens are grouped across different latent tokens and how the dynamics of the population are modeled in this latent space. Additionally, our proposed unit embedding space allows to map units into a high- sigh-dimensional space; thus understanding how unit projections are organized might help reveal the similarities between different neurons and the nature of their interactions. Similarly, can analyse the session embeddings to glean insights into inter-session and across -animal differences.
Our work shows how pretraining on diverse data, including datasets from animals performing different tasks and across different laboratories, can all help to improve our ability to decode from novel and unseen neural datasets. Already, our results demonstrate the positive effect of scale for neural data analysis. However, to scale this approach further and integrate even more diverse brain regions and tasks, it will be critical to move toward a self-supervised objective. Thus, current architecture and multi-session framework could be also extended to self-supervised tasks like generative next-event prediction or masked modeling to allow for even larger datasets to be ingested.
This framework has the potential to advance neuroscience research in several ways. By enabling the development of large pretrained models that can be fine-tuned for various downstream tasks, our work can accelerate progress brain-machine interfaces and disease modeling applications. The ability to quickly and effectively decode neural activity from diverse datasets and modalities can have significant implications for improving our understanding of the brain and developing new therapeutic interventions.
10

---

Acknowledgements
We would like to thank Patrick Mineault for providing a lot of great feedback on the work, and Jingyun Xiao for helping to run baselines for NLB. This project was supported by NIH award IROIEB029852. 01, ISF award IIS-2146072 as well as generous gifts from the Alfred Sloan Foundation (ED), the McKnight Foundation (ED), and the CIFAR Azrieli Global Scholars Program (ED). BAR acknowledges support from NSERC (Discovery Grant: RGPIN-2020-0N 105; Discovery Accelerator Supplement: RGPAS-2022-00033 2020-00031; Arthur B. McDonald Fellowship: 566355 2022) and CIFAR (Canada AI Chair; Learning in Machine and Brains Fellowship). MGP acknowledges support the Fonds de recherche du Quebec Santé (Grant: Chercheurs-boursee acknowledge support from the Canada CIFAR Al Chair program, the Canada Research Chair for Neural Computation and Interfacing, and NSERC Discovery grant HEPIN-2018-04822
en intelligence artificielle). GL
References
[1] A. Radford, K. Narasimhan, T. Salimans, I. Sutskever, el al., "Improving language understanding by generative pre-training"" 2018. 1
[2] A. Chowdhery, S. Narang, J. Devlin, M. Bosma, G. Mishra, A. Roberts, P. Barham, H. W. Chung, C. Sutton,
S. Gehrmann, al., "Palm: Scaling language modeling with pathways." arXiv preprint arXiv2220.003311, 2022.
[3] M. Dehghani, J. Djolonga, B. Mustafa, P. Padlewski, J. Heek, J. Gilmer, A. Steiner, M. Caron, R. Geirhos. I. Alabdulmohsin, et al. "Scaling vision transformers to 22 billion parameters." arXiv preprint 2xXiv:3322.05422 2023. 1.2
[4] A. Radford, Kim, Xu, G. Brockman, C. McLeavey, and Sutskever, "Robust speech recognition via large-scale weak supervision," arXiv preprint arXiv:22/2.043356, 2022.
[5] R. Bommasani, D. A. Hudson, F. Adeli, R. Altman, S. Arora, S. von Arx, M. S. Bernstein, J. Bohg, A. Bosselut, E. Brunskill, et al. "On the opportunities and risks of foundation models." arXiv preprint arXiv:2108.0724 07258, 2021. 1
[6] J. Kaplan, S. McCandlish, T. Henighan, B. Brown, B. Chess, R. Child, S. Gray, A. Radford, J. Wu, and D. Amodei, "Scaling laws for neural language models," arXiv preprint arXiv:2001.0883 2020.
[7] A. E. Urai, B. Doiron, A. M. Leifer, and A. K. Churchland, "Large-scale neural recordings call for new insights to link brain and behavior," Nature neuroscience, vol. 25, no. 1, pp. 1 11-19, 2022.
[8] C. Hurwitz, N. Kudryashova, A. Onken, and M. H. Hennig, "Building population models for large-scale neural recordings: Opportunities and pitfalls," " Current Opinion Neurobiology, vol. 70, pp. 64-73, 2021. Computational Neuroscience.
[9] M. Dabagia, K. P. Kording, and E. L. Dyer, "Aligning latent representations of neural activity." Nature Biomedical Engineering, pp. 1-7, 2022. 1,2
[10] C. Pandarinath, K. C. Ames, A. A. Russo, A. Farshchian, L. E. Miller, E. L. Dyer, and J. C. Kao, "Latent factors and dynamics in motor cortex and their application to wriin-achhnn interfaces." " Journal of Neuroscience, vol. 38, no. 44, pp. 9390-9401, 2018. 1, 5, 16
[11] W. M. Grill, S. E. Norman, and R. V. Bellamkonda, "Implanted neural interfaces: biochallenges and engineered solutions," Annual review of biomedical engineering, vol. 11, pp. 1 1-24, 2009. 2
[12] J. Jude, M. Perich, L. Miller, and M. Hennig, "Robust alignment of cross-session recordings of neural population activity by behaviour via unsupervised domain adaptation." " in Proceedings of the 39th Interna- tional Conference on Machine Learning (K. Chaudhuri, S. Jegelka, L. Song, C. Szepesvari, G. Niu, and S. Sabato, eds.), vol. 162 of Proceedings of Machine Learning Research, pp. 10462-10475, PMLR, 17-23 Jul 2022. 2, 10
[13] A. Jaegle, S. Borgeaud, J.-B. Alayrac, C. Doersch, C. Ionescu, Ding, S. Koppula. D. Zoran, A. Brock, E. Shelhamer, et al., "Perceiver io: general architecture for structured inputs & outputs," arXiv preprint mxXiv:2107.14995 2021. 2, 15
[14] A. Vaswani. N. Shazeer, N. Parmar. J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, "Attention is all you need," Advances neural information processing systems, vol. 30, 2017. 2,
11

---

[15] T. Brown, B. Mann. N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry. A. Askell, el al., "Language models are few-shot learners," Advances in neural information processing systems, vol. 33, pp. 1877-1901, 2020. 2,
[16] H. G. Rey. C. Pedreira, and R. Q. Quiroga. "Past, present and future of spike sorting techniques," Brain research bulletin, vol. 119, pp. 106-117, 2015. 2
[17] M. Pachitariu, N. A. Steinmetz, S. N. Kadir, M. Carandini, and K. D. Harris, "Fast and accurate spike sorting of high-channel count prohes with kilosort," Advances in neural information processing systems, vol. 29, 2016."
[18] A. Jaegle. F. Gimeno, A. Brock, O. Vinyals, A. Zisserman, and Carreira, "Perceiver: General perception with iterative attention," in International conference on machine learning, pp. 4651-4664, PMLR, 2021. 4. 9, 14
[19] Su, Lu. S. Pan, A. Murtadha, B. Wen, and Liu, "Roformer: Enhanced transformer with rotary position embedding," arXiv preprint arXiv:2104.09833 2021. 4, 14
[20] G. Perich, J. A. Gallego, and E. Miller, neural population mechanism for rapid learnin."" " Newron, vol. 100, no. 4, pp. 964 976, 2018. 5, 6, 16, 17
[21] M. G. Perich and L. Miller, "Altered tuning in primary motor cortex does not account for behavioral adaptation during force field learning," Experimental brain research, vol. 235, no. 9, pp. 2689-2704, 2017.
5, 17
[22] I. Glaser, M. G. Perich, Ramkumar. L. Miller, and K. P. Kording, "Population coding of conditional probability distributions in dorsal premotor cortex." Nature communicutions, vol. 9, no. 1, p. 1788, 2018. 5, 11
[23] J. A. Gallego, M. G. Perich, R. H. Chowdhury, S. A. Solla, and L. E. Miller, "Long-ten term stability of cortical population dynamics underlying consistent behavior," Nature neuroscience, vol. 23, no. 2, pp. 260-270000000 260-270, 2020. 5, 10, 17
[24] M. M. Churchland, J. P. Cunningham, M. T. Kaufman, J. D. Foster, P. Nuyujukian, S. Ryu, and K. V. Shenoy, "Neural population dynamics during reaching," Nature, vol. 487, no. 7405, pp. 51-56, 2012. 5, 18
[25] J. G. Makin, J. E. O'Doherty, M. M. Cardoso, and P. N. Sabes, "Superior arm-movement decoding from cortex with a new, unsupervised-eeaminn algorithm." Journal of neural engineering, vol. 15, no. 2. p. 026010, 2018. 18
[26] R. D. Flint, E. W. Lindberg, L. R. Jordan, L. E. Miller, and M. W. Slutzky, "Accurate decoding of reaching movements from field potentials in the absence of spikes," Journal of neural engineering, vol. 9, no. p. 046006, 2012. 5, 18
[27] F. Pei, J. Ye, D. Zoltowski, A. Wu, R. H. Chowdhury, H. Sohn, J. E. O'Doherty, K. V. Shenoy, M. T. Kaufman, Churchland, al., "Neural latents benchmark 21: evaluating latent variable models neural population activity," arXiv preprint (tv:2109.04466, 2021. 5, 8, 16, 18
[28] Y. You, Li, S. Reddi, J. Hseu, S. Kumar, S. Bhojanapalli, X. Song, Demmel, K. Keutzer, and C.-J. Hsiih, "Large batch optimization for deep learning: Training bert in 76 minutes," arXiv preprint arXiv: reee::9049990 2019. 6, 15
[29] G. Perich, L. Miller, M. Azabou, and Dyer, "Long-term recordings of motor and premotor cor- tical spiking activity during reaching in monkeys." https:// daniiarchiveeeeee 2023. 6, 17
0000688
[30] J. I. Glaser, A. S. Benjamin, R. H. Chowdhury, M. G. Perich, L. E. Miller, and K. P. Kording, "Machine learning for neural decoding," Eneuro, vol. 7, no. 4, 2020. 6, 18
[31] M. R. Keshtkaran, A. R. Sedier, R. H. Chowdhury, R. Tandon, D. Basrai, S. L. Nguyen, H. Sohn, M. Jazayeri, L. E. Miller, and C. Pandarinath, "A large-scale neural network training framework for generalized estimation of single-trial population dynamics." Nature Methods, pp. 1-6, 2022. 8, 19
[32] J. Ye and C. Pandarinath, "Representation learning for neural population activity with neural data trans- formers," arXiv preprint arXiv:2108.01210, 2021. 8. 10
[33] R. Liu, M. Azabou, M. Dabagia, J. Xiao, and E. Dyer, "Seeing the forest and the tree: Building represen- tations of both individual and collective dynamics with transformers." Advances in Neural Information Processing Systems, vol. 35, pp. 2377-2391, 2022. 8, 10
12

---

[34] G. Zerveas, S. Jayaraman, D. Patel, A. Bhamidipaty, and C. Eickhoff, "A Gansformer-based framework for multivariate time series representation learning." in Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pp. 2114-2124, 2021.
[35] Y. Nie. N. H. Nguyen, P. Sinthong, and J. Kalagnanam. "A time series is worth 64 words: Long- forecasting with transformers." " International Conference on Learning Representations. 2023."
[36] X. Shou, T. Gao, D. Subramanian, D. Bhattachariya, and K. Bennett, "Influence-aware attention for multivariate temporal point processes," " in 2nd Conference on Causal Learning and Reasoning, 2023,
[37] X. Shou, I. Gao, S. Subramaniam, D. Bhattachariya, and K. Bennett, "Concurrent multi-label prediction in event streams," " in AAAI Conference on Artificial Intelligence, 2023, 10
[38] S. Zuo, H. Jiang, Z. Li, T. Zhao, and H. Zha, "Transformer hawkes process," in International conference on machine learning, pp. 11692-1 11702, PMLR, 2020.
[39] T. Le and E. Shlizerman, "Stadt: Modeling neural population activity with spatiotemporal transformer," arXiv preprint arXiv:2206.04727, 2022. 10
[40] R. Liu, M. Azabou, M. Dabagia, C.-H. Lin, M. Gheshlaghi Azar, K. Hengen, M. Valko, and E. Dyer, "Drop, swap, and generate: supe-super approach generating neural activity." Advances in neural information processing systems, vol. 34, pp. 10587-10599, 2021. 10
[41] C. Pandarinath, D. J. O'ha,,, Collins, R. Jozefowicz, D. Stavisky, J. C. Kao, M. Trautmann, M. T. Kaufman, S. I. Ryu, L. R. Hochberg, et al., "Inferring single-trial neural population dynamics using sequential auto encoders,, Nature methods, vol. 15, no. 10, pp. 805 815, 2018. 10
[42] A. Farshchian, J. A. Gallego, J. P. Cohen, Y. Bengio, L. E. Miller, and S. A. Solla, "Adversarial domain adaptation for stable brain-machine interfaces," in International Conference on Learning Representations, 2019. 10
[43] S. Schneider. J. H. Lee, and M. W. Mathis, "Learnable latent embeddings for joint behavioural and neural analysis," Nature, pp. 1-9, 2023. 10
[44] E. L. Dyer, M. Gheshlaghi Azar, M. G. Perich, H. L. Fernandes, S. Naufel, L. E. Miller, and K. P. Körding, "A cryptography based approach for movement decoding," Nature biomedical engineering, vol. 1, no. 12, pp. 967-977779 976, 2017. 10
[45] A. D. Degenhart, W. E. Bishop, E. R. Oby, C. Tyler-Kabara S. M. Chase, A. P. Batista, and B. M. Yu, "Stabilization of brmin-computer interface via the alignment low-dimensional spaces neural activity," 2020, 10
Nature biomedical engineering, vol. 4, no. 7, pp.
[46] X. Ma, Rizzoglio, K. L. Bodkin, E. Perreault, L. E. Miller, A. Kennedy, "Using adversarial networks to extend brain computer interface decoding accuracy over time," eLife. vol. 12, p. e84296, aug 2023. 10
[47] Safaie, C. Chang. J. Park, E. Miller, T. Dudman, M. G. Perich, and A. Gallego, "Preserved neural population dynamics across animals performing similar behaviour," bioRxiv, 2022. 10
[48] Azabou, M. Azar, Liu, C.-H. Lin, E. Johnson, K. Bhaskaran-Nait M. Dabagia, B. Avila- Pires, L. Kitchell, K. B. Hengen, W. Gray-Roncal, M. Valko, and E. L. Dyer, "Mine your own view: Self-supervised learning through across-sample prediction," arXiv preprint 20X 2021. 16
[49] H. Banville, O. Chehab, A. Hyvarinen, D. Engemann, and A. Gramfort, "Uncovering the structure of clinical EEG signals with self self-super learning." Journal of Neural Engineering, 2020. 16
13

---

Appendix
A
Additional Model Details
A.1 Rotary position encoding
Rotary Position Embeddings (RoPH) used to incorporate relative positional information into transformer models in a memory and compute efficient manner[19) Unlike absolute position embeddings which are added appended the input embeddings, these are directly incorporated in the attention mechanism.
First, we define 2 x 2 rotation matrix R2x2(t, T) for given timestamp and time- tine-period T:
[cos(2 (2Tl/T) sin ( cos(2xl//7)
sin ( (T)
R2x2 T)
R2x2(t, (t, T)) provides time/positional encoding with time period T. Intuitively, the sinusoidal values in this rotation matrix can help the model learn and distinguish variations in timestamps at a timescale T.
Given the head dimension of the attention block D, we create D x rotary embedding matrix R(t) for timestamp t by constructing block-diagonal matrix of rotation matrices R2x2, where the time-periods are varied logarithmically between Twin and T man
(t,
0
(t, Ti)
R(1)-
0
R2x2 (1, To
Given D, T Tinss and Tmazi we generate T; by the following equation:
T Tmaz
2i/D
T: = T Tii
Tiin
In our models, T Tmin and T maz were set Ims and respectively.
This rotation matrix is incorporated into the attention mechanism by replacing the attention score calculation with the following equation:
a.j = softmax (R((t)):) (R(t,)k,)
Since R(t,)"R(tt) = R(t, -t), tz), the attention score a, depends only on the relative timing between the query and key tokens. This is what makes RoPE relative in nature.
What is explained above has been introduced in [19], and use successfully in many models [18]. In experiments, we found an interesting extension of RoPE. In addition to incorporating RoPE in values inside the attention mechanism. We do this to encode the relative timing between values and queries in the output of the attention block. A direct, but computationally inefficient, way of doing so would be:
calculating the attention scores, we also incorporate
the reiigheed-summatin
23 R3 + Σ a R(4 4,)x
This is inefficient since it requires us to explicitly compute R(t, -ty) for all N values of (t, t2). use the fact that R(t: - (,) can be factorized as R(-t,) to reach the following, much more efficient, implementation:
2 - Zj R R(-1;)) Σ a R(,)(
That is, we first apply rotation matrices R(1,) to all values v; before passing them to inputs the attention block. Then apply the rotation matrices R(-t,) to its outputs.
14

---

This value- rotationn mechanism captures Sllativettiming information more explicitly in the output of the attention block, allowing the feedforward network to take advantage of timing information. In the absence of this mechanism, the only way timing information is added to the latent tokens is through the attention scores. We observed in early experiments that our model performed better with value-rotation enabled for the input cross-attention and all self-attention layers.
As discovered in the Perceiver framework [13], we observed an improvement in performance upon rotating only half of the head dimensions. That is, half of the R2x2 matrices in R. were replaced by 2 x identity matrices.
A.2 Time-based context window
We consider an arbitrary window of time [t, t + 7], meaning that we do not use the trial structure when training our model. In our experiments, we set T equal to 1 second. At inference time, we use a sliding window of step 500ms and average the predicted output in overlapping segments.
In language models, and more generally transformer-based models, the context window is defined with respect to the number of tokens. The neural data modality is a bit more complex, as the number of spikes scales with the number of units: Across multiple datasets, we record variable number of units, ranging from tens of units to hundreds of units. Defining the context window with respect to number of tokens will result in a variable time duration, long for small populations and short for larger ones. Our context window defined with respect to time (fixed to 1s), each sequence will have a variable number of tokens (M). When creating batches, we pad to the largest sequence in the batch, and appropriately use an attention mask in the first cross-attention layer. GPU hardware is optimized to process fixed-size tensors, to make computation more efficient, we introduce a distributed training load balancing mechanism, which deals with variable-length input sequences (Appendix A.4.4).
A.3 Further details tokenization
Delimiters. When extracting the spike data in a window of time [t, 1 + T), we include delimiters for each unit, indicating the model that we tuning into the activity of that unit between times t and t + T. This choice is important because: 1) we are using relative position encoding, so the beginning or the end of the window is not clear, and 2) our tokenization only captures activity and not inactivity. If we consider a unit that is inhibited and does not fire in that window of time, the model would be unaware of the presence of that unit, unless we add [START] and [END] tokens. We learn two embedding vectors for these delimiters Xstart and Xend- For each unit in the population, we add these two tokens with embedding equal to the sum of the delimiter embedding and the unit embedding, and assign each token timestamps and t + T respectively.
Latent tokens. Recall that we divide the latent tokens into groups of and spread each group uniformly over the context window. Specifically, we use a total of 256 latent tokens and divide them into groups of 8, this means that there will be 32 tokens that share the same timestamp. We experiment with weight-sharing, where we set the learned latent tokens in the same group to have the same weights (Co. ... C31 in Fig 1). This means that tokens in the same group (size 8) will have the same embedding Zo,; e( mod 32). but different timestamps (t)....,tt 1 in Fig 1). We believe this a reasonable inductive-bias to architecturally enforce since, intuitively, we want all sets of latent tokens that share a timestamp to query the spike tokens in the exact same manner, just with a different timestamps. These tokens will evolve independently once the information is pulled from the input space (after the first cross-attention).
A.4 Training details
The model is trained using the LAMB optimizer [28] with weight decay. The learning rate is held constant, then decayed towards the end of training (last 25% of epochs), using a cosine decay schedule. Single-session models are trained on single GPU with batch size of 128 while large models are trained with 8 GPUs with a total batch size of 1400. Note that we didn't see any benefits increasing the batch size when training single-session models.
15

---

A.4.1 Compute
The large models are trained on machine with an AMD EPYC 7452 32-Core Processor and 8 Nvidia A40 GPUs (48Gb memory). PD was trained for 2 days, and POYO-1 was trained for days (both for a total of 400 epochs). Single-session models are trained with single Nvidia GeForce RTX 3090 GPU, and take less than an hour to train. Finetuning models is also done with a single GPU. Unit identification converges very quickly and takes less than minute on a single GPU or a few minutes on the CPU.
A.4.2 Data augmentation
Previous work [48, 49] explores the use of augmentations when training neural population models. In our experiments, we use the unit dropout augmentation which randomly samples subset of the population. We set a minimum population size of 30 for this augmentation to ensure that it is not too destructive (note that our model is still trained with any number of units, this limit is only imposed for the augmentation). We did not explore the use of other augmentations, but believe it could be promising direction to further improve the capabilities of the model.
A.4.3 Loss
We train our model using mean-squared error loss over the hand velocity sequence, We normalize the velocity data using z-scoring the lab level, this means that all velocity data in datasets from the same lab is rescaled using the same scalar.
Dealing with variable sampling rate. When training on mixtures of datasets from multiple labs, any given batch will contain behavior sequences that have different sampling rates. In a window of 1s, some sequence will have more timepoints over which we will evaluate the MSE compared to other sequences. To avoid the over-representation of samples that simply have a high sampling rate, we evaluate the error on a randomly sampled 100 timepoints (these points are randomly sampled at each pass).
Weight in training mix. While behavior during random target tasks can be more complex and noisy, behavior during center-out reaching tasks is more structured, the monkey usually gets preparation phase, where the movement can be planned [20]. The latter task has been studied extensively for that reason. Knowing that neural activity is very salient during center-out reaching, weight of the prediction loss during the reaching segments by a factor of 5. The idea of over- ever-weighing samples of "good quality" is not novel and is commonly used when training language models [15].
increase the
A.4.4 Load Balancing
We distribute our training over multiple GPUs. Note that, in our experiments, the input sequence length (M) can vary anywhere between 1K and 20k tokens. Because of the variability in the length of the input sequences, we distribute sequences that have close length same GPU node.
We now provide further details on our load balancing strategy:
Each GPU node has local batch size. A GPU that is assigned large sequences will have smaller local batch size than one that assigned shorter sequences.
We create number of data buckets equal the number of GPUs. Each bucket i is assigned sequences that are of length M E [MMinn Minal, and local batch size B,.
We select these parameters in such way that we effectively utilize all available compute capacity, and minimize the amount of padding used.
A.5 Evaluation Details
During training, we train on any Is segment of the recording, but at evaluation we report the decoding
performance using the same evaluation strategy used in previous work [27, 10]: for center-out reaching datasets, we report the decoding score during the reaching movement (only when the trial is completed successfully), for random target datasets that are trialized (hold period followed by 3/4 random reaches), report the decoding score during the reaching as well, and finally for all other continuous random target datasets, we report the decoding score for all segments. Note that none of the segments used for testing seen during training.
16

---

CS
Monkey C (Left M1 & PMd)
IPS
AS
Monkey Cp (Right M1)
Perich, et al.
Monkey (Left M1 PMd)
Monkey Ja (Left M1)
PMd
S1
PCD
Monkey T (Left PMd)
M1
Monkey Je (Left M1)
Churchland, et al.
Midline
Monkey (Left M1)
M1
Makin et al.
Monkey (Right M1 S1)
PMd
(estimated locetions)
Monkey L (Left M1 & S1)
Flint, et al.
Monkey Mi (Right M1)
Figure 6: Diversity of neural recordings used to train POYO. We visualized the array locations for all of the datasets used in our current study. POYO-mp is trained on recordings from PMD and MI (CL, CR. M, Ja) and POYD-1 is trained recordings from (CL, CR, M, Ja, I, L, Mi) spanning PMD, M1, and across both hemispheres and in seven animals.
A.6 Fine-tuning Details
When finetuning the model, perform gradual unfreezing of the weights. For the first few steps, optimize the unit embeddings and the session emhedding only (unit identication) and then unfreeze the rest of the weights. We use the LAMB optimizer with batch size 128, learning rate 10-1 and weight decay 10-4 We use the exact same hyper parameters and finetuning process for all experiments we report. This strongly suggests that our method general, and easy adapt new data, without the need for expert neural network training knowledge, or expensive hyperparameter search.
B Datasets
To train and evaluate our framework, we curated large number of open and publicly available datasets that represent range of neural and behavioral recording platforms.
B.1 Datasets used for training
Datasets used train POYO-mp. To build our first large-scale model (POYO-mp), we acquired 100 sessions of unique recordings from three non-human primates performing two different movement tasks that previously studied in four publications [21, 22, 20, 23]. We make these datasets publicly available on the DANDI open repository [29]. In all of these datasets, nonhuman primate seated in primate chair and executes movements using custom 2-D planar manipulandum to control a cursor on a computer screen. They performed one of two tasks within a session.
Center-out Task (CO): They performed a standard center-out reaching task involving eight outer targets arranged around a circle with radius of 8 cm. The monkeys were required to move to the center target, wait for a variable hold period, and then move to one of the outer targets upon receiving auditory go cue. The task aimed to study neural activity during movement planning, preparation, and execution.
Random Target Task (RT): The setup is similar, the targets are not arranged around a circle but randomly placed. The monkeys required to move between 4 targets upon receiving an auditory go cue.
After extensive training, the monkeys were implanted with chronic multi-electrode arrays in the primary motor cortex (MI) and the dorsal premotor cortex (PMd). Neural activity recorded using
17

---

a Blackrock Cerebus system, and the data was manually processed offline to identify single neurons. To ensure only independent channels were included, steps were taken to identify and disable potential candidates with high crosstalk and to exclude cells with high percentage of coincident spikes.
Datasets used train the multi-lab, multi-session model (POYO-1). When training POYO-1, we acquire more data as shown in Table 1, adding all of the datasets outlined from Churchland et al., [24], Makin et al., [25], and Flint et al., [26]. We still hold out the 12 sessions from Monkey sessions from Monkey C, and the two Neural Latents Benchmark (NLB) [27] datasets from training.
As depicted in Figure instead of using a manipulandum like in the previous experiments. these new datasets use a touch screen and have different sampling rates for their behavioral outputs. addition, in the Churchland, et al. datasets, they threshold crossing based processing rather than applying spike sorting to identify single units. We refer the reader to the listed publications for more details about how these datasets were collected.
We would like to reiterate that there are major differences in the way the datasets were collected. Yet, in spite of this variability, it is important that our method works without having to standardize the datasets across sites further and apply specialized techniques to one dataset but not other. Thus we did not attempt to homogenize or standardize the data. The only processing that we perform is on the behavior to clip outliers in the behavior based upon times of high acceleration. We use a simple threshold heuristic to avoid training on periods with extremely high acceleration.
This means that:
1. We do not filter units based presence ratio, reject multi-units.
2. We do not re-apply the same spike sorting algorithm on all datasets, and use the data as is. In [24], the spiking events obtained through threshold crossings (i.e. all units are multi-units), and in [25], have both threshold crossing units and spike sorted isolated single units, and we use both. Each algorithm used in each lab was tuned differently, but we consider this a good example of diversity that general model needs to deal with.
3. We do not resample process the velocity timeseries.
Minimizing the amount of processing needed to integrate new datasets into our model is key for more easily scaling to more datasets, and also providing an accessible model that the community at large can use out of the box, without having to adapt to specific standard.
B.2 Datasets for evaluation and fine-tuning
In addition to holding out 20% of the data from each session, we hold out all sessions from Monkey T (6 CO, 6 RT). We also use two datasets from the Neural Latents Benchmark [27], MC-Maze and RTT. The Maze dataset consists of recordings from primary motor and dorsal premotor cortices of a monkey performing reaches with instructed delays visually presented targets while avoiding the boundaries of a virtual maze. This dataset offers a variety of behavioral configurations, with each configuration using different combination of target position, number of virtual barriers, and barrier positions, resulting in a variety of straight and curved reach trajectories. With thousands of trials, the Maze dataset allows for rich investigation into the structure of population activity, while its instructed delay paradigm enables clear separation of neural processes related to preparation and execution. On the other hand, the RTT dataset introduces different modeling challenges, as it contains continuous, point-to- to-nnnt reaches that start and end in various locations, have highly variable lengths, and few repetitions. We report the performance of our models on the same segments of data used to evaluate models in the NLB benchmark.
C
Additional Results
C.1 Comparison with single-session baselines
In our comparisons Table 2, we compared PΟYO against existing baseline models that are commonly used for neural decoding [30]. In particular, we applied a Wiener Filter (WF), Feed Forward Network (MIP), and a Gated Recurrent Unit (GRU). These models predict the behavior at each timestamp t given a history window of neural activity [t T, t). For the NLB datasets, we also compare with
18

---

Auto-LFADS [31], unsupervised approach that can be used to obtain single trial smoothed rate estimates for a population of neurons before applying linear layer on top of the inferred rates to obtain an estimate of the behavior. While this model isn' 't designed for decoding per se, it is frequently used for denoising for BMI decoding applications.
For the NLB, we use a 5 ms binning of the neural activity for training models and for evaluation of the behavioral output as this the finest resolution used in the benchmark. For the rest of the MP datasets, we use a 10 ms binning rate for the neural data and behavior as this corresponds to the behavior sampling rate for these datasets.
The hyperparameters we used for training our single-session POYO model can be found in Table 4.
Hyperparameter
Value
Embedding Dimension
128
Head Dimension
64
Number of Latents
128
Depth
6
Number Heads
8
FFN Dropout
0.3
Linear Dropout
0.3
Attention Dropout
0.3
Weight Decay
le-4
Learning Rate
4e-3
Batch Size
128
Table 4: Hyperparameters used for training all POYO single-session models
C.2 Robustness of Unit-identification
Unit-identifiiation, introduced in Section 2.5, can be seen as an approach for identifying units from new dataset. To evaluate the robustness of this process, we ran unit-identificttio was already seen during pre-training. We test whether finetuning on this "new" set of units, will map them close to their true embedding (found during pre-training). We use our largest model, POYO-1, start with randomly initialized unit embeddings and run unil-identification.
session that
We perform this experiment on two different animals. In Figure 7, we show the evolution of unit- embeddings during training. To compare the newly calibrated set of units with its pre-trained version, we normalize the umit-embeddings and report the cosine similarity averaged over the set of units. Note that normalization done because POYO's first cross-attention layer applies layer norm to these embeddings. We also report the "unit identification accuracy" which we define as the ratio of tuned unit embeddings that have their true unit embedding as their nearest neighbor. Results for these experiments are present in Table 5.
Dataset
Num of Units
Unit-Id Accuracy
Unit-Id Cosine-Similane
Monkey C, CO, 2013/10/03
73
1.000
0.845
Monkey M, CO, 2014/02/03
116
0.966
0.802
Table 5: Unit re-identification results. Cosine-similarii similarity measure is averaged over all units in dataset. Accuracy measured for neighoottn classifier.
Overall, this analysis suggests that the unit-identification approach is robust and reliable for identifying units. We believe that understanding the functional similarities of units that are mapped close each other in the unit embedding space will be an interesting avenue for future work.
C.3 Hyperparameter sensitivity
When initially training the single-session models, we selected subset datasets (out the 100 sessions) and performed a random hyperparameter search and selected the set of hyperparameters with the best performance on the validation set. These experiments quickly revealed that the model was very stable to wide range of hyperparameters. In Figure 8, report some examples of this for
19

---

Monkey M, CO 201 4/02/03
Epoch
Epoch 100
Fasch 853
60-
-02
-04
Monkey 2013/10/03
Epoch
Epach
Spech 100
-22
Figure 7: Evolution of the unit embeddings with unit-identification sessions from two different animals. From left to right, the finetuned embeddings (orange) are visualized 10, 100, and 400 epochs, relative their true embedding in the pre-trained model (blue).
two sessions. We perform 300 runs for each session, then report the average performance for each hyperparameter.
Given these findings, we selected the best set hyperparameters for the subset and used the same parameters when training all of the 100 single-session models that study in Figure 2.
Validation R^2 averaged over multiple runs
Monkey C 2016/10/13
Monkey T 2013/08/19
1.00
0.95
0.85
0.80
10
200
400
100
200
10
15
Model depth
Batch size
Embedding dimension
Number of heads
Figure 8: Single-session model performance for wide range of values for different hyperparameters. The R2 score reported on the validation sets for both datasets.
20


Reviews:
[{"review_id": "mEenZuJx7D", "reviewer": "Reviewer_c1Lq", "summary": "The authors describe a novel method for the task of neural decoding: using the time series of activity recorded from a population of neurons to predict the activity of scientifically relevant target variables. The describe their approach, called POYO, based upon the tokenization of spike data, the application of transformer based models to this data to generate pre-trained models, and the use of these pre-trained models for neural decoding. They demonstrate that this approach allows for accurate generalization across different sessions of neural recordings from the same animal, and across different animals. ", "questions": "- How many steps/epochs are fine tuning and unit identification run for? I was unable to find this information in the appendix. \n- Line 55: the phrase \"neural scaling laws\" is ambiguous. Can it be changed to \"scaling laws\"? \n- Line 349-350: \"AutoLFADS extends the LFADS framework with Population Based Training to perform hyperparameter tuning, but has not been demonstrated for multi-session alignment.\"Is there any reason to believe that AutoLFADS is not capable of performing multi-session alignment? To my knowledge, Auto-LFADS is a optimization routine for LFADS. I would assume that it works for multi-session recordings as well. \n", "limitations": "As noted by the authors, generalizing this framework across other brain regions (as well as many other conditions) is a key limitation. In general, I have reservations about the use of the phrase \"foundation model\" without qualification in the context of neural decoding. Is the suggestion of this work that neuroscientists should aim to build a general purpose model that performs decoding regardless of model organism, brain area, or recording modality? While the results of this study are impressive, I do not believe that such claims are justified, even as future work. ", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The authors present a novel approach for the decoding of neural data, that appears promising in its ability to generalize over different animals and experimental sessions, which is a major problem in neuroscientific data analysis. The performance of their method is impressive, especially in the setting of pretrained models that are finetuned to novel datasets. The ability to combine various datasets in the way described here has the potential to be impactful. ", "weaknesses": "My main concerns are with regard to scientific rigor, and the scope of the claims made in this work.\n- Generalization across brain regions. I believe there needs to be a more thorough discussion of the role of brain regions in tasks where POYO is asked to generalize across individuals with recordings from different brain areas. To what degree will region-specific neural activity impact the ability of POYO to generalize when performing few shot learning? Concretely, what would behavioral decoding results in Table 2 look like if they are broken out per individual brain areas? \n- Lines 37-39: \"Overall, this lack of correspondence across recordings and channels complicates the integration of information from different experiments and individuals, ultimately hampering efforts to construct a unified perspective on population-level interactions and dynamics in the brain.\"It is not clear to me how POYO will address these challenges. While it is clear that POYO allows for highly performant decoding, the lack of recurrent structure within POYO makes it difficult to understand how this model would assist these challenges. Perhaps analysis of the attention mechanism could provide some answers to this question, but that point is not discussed in this work. \n- Comparison with existing baselines. In Table 1, it is unclear to me why the performance of AutoLFADS + Linear only provided for the NLB datasets. It would be useful to see the performance of this model on all datasets, as it is a more powerful and commonly used decoding model than the others provided. Comparison with other transformer based neural decoding approaches like NDT (Ye and Pandarinath 2021) would also provide a better perspective on the value of the methodological advances proposed in this particular work. \n\nI am open to reconsidering my score if my concerns here are addressed. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments and questions, and for finding that our work \"has the potential to be impactful\"! In what follows, you will find our point-by-point response to your main concerns, and results from new experiments that we ran to address your feedback. \nPlease let us know if there’s anything we can clarify further.\n\n> 1. “I believe there needs to be a more thorough discussion of the role of brain regions in tasks where POYO is asked to generalize across individuals with recordings from different brain areas.”\n\n**Reply:** Thank you for your suggestion. In the final paper, we plan to add a detailed discussion of the functional roles of the areas studied, and provide a new visualization of the recording sites across all the datasets that we study (see Figure 2 in attached PDF). Our current efforts have focused on neural recordings in three regions—M1, PMd, S1—that are intimately involved in planning and executing movements. While their functional role in movement and sensation are unique, many studies including ours have shown that they can all provide information that is used to guide decoding of movement. Yet, while different regions serve specific functions during voluntary movement (e.g., contracting muscles and sensing the state of the limb), they may perform these functions using common principles in the neural code.\n\n> 2. “To what degree will region-specific neural activity impact the ability of POYO to generalize when performing few shot learning?”\n\n**Reply:**  The datasets we consider in this paper broadly sample reach-related activity across PMd, M1 and S1. These three regions are highly distinct both functionally and based on cytoarchitecture and anatomy. For example, M1 and S1 have been shown to have highly distinct dynamics underlying their population activity. Furthermore, even within a region such as M1, our recording sites are distributed across varying functional sub-regions (Figure 2) and represent a great deal of diversity in terms of the underlying neural functions and firing patterns measured.\n\nIn terms of across-region transfer of the model, we believe it is an important question that requires very thorough investigation, and a careful selection of pre-training and fine-tuning datasets which is outside of the scope of this work. We will make sure to include many of these points in our discussion.\n\n> 3. “Lines 37-39: \"[...] The lack of recurrent structure within POYO makes it difficult to understand how this model would assist these challenges.”\n\n**Reply:** While we have not yet tackled this question, as the reviewer notes, analysis of the attention heads provides one possible way of using POYO to understand population-level interactions. We will add some discussion of this possible use to the manuscript.  \nWith regards to the comment that there is a “lack of recurrent structure within POYO”, the transformer architecture that we leverage provide powerful and general building blocks for modeling complex sequential data and time-series. Perhaps we misunderstood your point. If you can clarify, we will do our best to address your concern!\n\n> 4. “In Table 1, it is unclear to me why the performance of AutoLFADS + Linear only provided for the NLB datasets.” \n\n**Reply:** AutoLFADS is one of the main baselines in the NLB benchmark and the authors provide code and hyperparameters for the NLB datasets. Thus, we were able to reproduce these numbers. However, we were not able to obtain good results without further hyperparameter tuning on the other datasets in Table 2.\n\nSince the original submission, we have been able to run our datasets on NeuroCAAS, which is a cloud service containing the AutoLFADS tuning procedure, to obtain results for 4 of the datasets below.\n\n| | Monkey C, CO 10/13 | Monkey C, CO 10/21 | Monkey T, RT 08/20 | Monkey T, RT 09/06 |\n| - | - | - | - | - | \n| AutoLFADS + Linear | 0.9292 | 0.9519 | 0.4116 | 0.5061 |\n| POYO | 0.9603 | 0.9759 | 0.7986 | 0.8306 |\n\nAs can be seen, POYO outperforms AutoLFADS on all of the datasets tested, with large gaps on RT. We will continue these runs as they are compute intensive, and we have 14 additional datasets that we report results on in Table 2. \n\n> 5. Comparison with other transformer based neural decoding approaches like NDT would also provide a better perspective on the value of the methodological advances proposed in this particular work.\n\n**Reply:** Thank you for this suggestion. We were able to reproduce NDT on both NLB datasets, and run other transformer models (see Table 1 in the General Response). We extensively tuned the hyperparameters of all baselines to ensure a fair comparison. Our results suggest that POYO (both single-session or pretrained) outperforms other transformer approaches. \n\n** We would also like to note that we attempted to train NDT on the rest of the datasets, but we faced a similar challenge with hyperparameter tuning. We will continue to tune these baselines and hope to include them on all the evaluation data in Table 2.\n\n> 6. How many steps/epochs are fine tuning and unit identification run for? \n\n**Reply:** Fine tuning and unit identification are run for 50 epochs, with a batch size of 128. We have updated the appendix to include these details.\n\n> 7. “Is there any reason to believe that AutoLFADS is not capable of performing multi-session alignment? To my knowledge, Auto-LFADS is a optimization routine for LFADS.”\n\n**Reply:** Auto-LFADS has not been demonstrated in the multi-session condition and code for multi-session HPO has not been released. While, in theory, AutoLFADS can be applied to the multi-session stitching condition, it is unclear how to tackle the initialization of the alignment matrix for each session, or what HPO strategy to use with the presence of multiple datasets. \nWe also note that multi-session stitching has only been done in sessions that are recorded from the same subject in the same brain region (and same chronic implant) during a stereotyped center-out task."}, {"review_id": "ZzZC22YmCP", "reviewer": "Reviewer_Vmcn", "summary": "This paper introduces a novel method for developing models that can predict the activity of neural populations by learning from data recorded across different sessions, tasks, and animals. The method is uses a custom tokenization procedure for spikes and a deep neural network architecture based on PerceiverIO. The paper presents the core method, then describes a technique to reuse a large pre-trained model to learn efficiently to predict a new neuronal population/condition, and finally presents a series of validation experiments.\n\n", "questions": "- line 277: \"we see a bit of a dip in overall accuracy\": how much is a bit? is this data shown anywhere?\n\n- in figure 4, the blue/green colors are hard to see in grayscale. Can you choose colors that render better in grayscale?\n\n- in figure 4, the compute efficiency plots don't seem to be cited anywhere in the text?\n\n- line 289: \"is competitive with other baselines\". What other baselines?\n\n- in the paragraph \"Performance on new animals performing new tasks\", it is not clear how much data was using for training and how much for testing. Can you clarify?\n", "limitations": "Limitation of the present work are appropriately discussed. I see no potential issue with societal impact.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- the proposed technique is technically solid, elegant, and performs well in the experiments.\n\n- this technique addresses a real need in neuroscience labs, making the paper potentially very high impact.\n\n- I found the scaling analysis (lines 237-257) particularly instructive, and especially relevant because for a new method like this one the reader may wonder what type of performance they may expect on their own data, as a function of the architectural choices they would have to take when deploying the method. This is of course a very hard thing to assess correctly, but this type of analysis is a great starting point.\n\n- the paper is very well written.\n", "weaknesses": "- I have noticed a few minor issues (mostly related to clarity and plots) in section 3.4. See below under \"questions\".\n\n- The discussion of related literature seems to imply that all spike train analysis has always been done by binning spikes. It would be good to include at least a cursory mention of kernel methods for spikes and other binless measures: see for instance Paiva, A.R.C., Park, I., Príncipe, J.C., 2010. A comparison of binless spike train measures. Neural Comput & Applic 19, 405–419.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review of the paper! We are very excited to hear that you found the proposed technique to be “technically solid and elegant” and also agree that it “addresses a real need in neuroscience labs”.\n\nIn what follows, we will provide point-by-point replies to your questions.\n\n> 1. “The discussion of related literature seems to imply that all spike train analysis has always been done by binning spikes. It would be good to include at least a cursory mention of kernel methods for spikes and other binless measures”\n\n**Reply:**  Thanks for the suggestion! We will update the related work to include more context and mention other binless methods, such as the paper you cited and related work on kernel methods. When contrasting our spike-based tokenization scheme with other binning-based approaches, we were referring specifically to the deep learning for neural data literature here, where binning has been the approach to-date, because binned representations regularize the structure of spiking data and enable compatibility with existing machine learning frameworks, like MLPs, RNNs, and Transformers. \n\n> 2. “Line 277: \"we see a bit of a dip in overall accuracy\": how much is a bit? Is this data shown anywhere?”\n\n**Reply:**  This is a comment on the results from Table 2. More specifically, the two columns corresponding to Monkey T performing CO and RT tasks respectively. We use POYO-MP, which is pre-trained on multiple sessions from Monkey C and M, and transfer it to the new sessions from Monkey T, using unit identification (meaning that the transformer weights are frozen). We show that despite tuning a small number of parameters (unit embeddings and session embedding), we achieve an accuracy that is only 1% away from the accuracy we obtain when we train all weights from scratch, the latter setting corresponding to POYO-[Single-session]. This is the “dip” we are referring to. We will update the description to be more specific.\n\n> 3.  In figure 4, the blue/green colors are hard to see in grayscale. Can you choose colors that render better in grayscale?\n\n**Reply:** Yes, we will update the colors to work well in grayscale in a revision.\n\n> 4. In figure 4, the compute efficiency plots don't seem to be cited anywhere in the text?\n\n**Reply:** Thank you for pointing this out. We have revised the text to include a reference to Figure 4 and will include further discussion on the significance of this result in a final version of the paper.\n\n> 5.  Line 289: \"is competitive with other baselines\". What other baselines?\n\n**Reply:** We are referring to the baselines in Table 2, including the Wiener Filter, MLP, and AutoLFADS. We will make sure to revise this section to make it more clear.\n\n> 6. In the paragraph \"Performance on new animals performing new tasks\", it is not clear how much data was used for training and how much for testing. Can you clarify?\n\n**Reply:** We use the splits defined in the NLB benchmark. For NLB Maze, there are 1721 training samples and 574 testing samples. For NLB RTT, there are 810 training samples and 270 testing samples.\n"}, {"review_id": "3l0CRtnVdX", "reviewer": "Reviewer_BY6C", "summary": "Deep learning and transformer models have shown great promise in identifying structure from large datasets. With recent advancements in neural recording methods, it is now possible to generate rich and heterogeneous recordings from large populations of neurons across multiple brain regions and experimental conditions. The proposed work bridges this gap by introducing a new approach to modeling neural dynamics across these recordings. The authors leverage transformer models by tokenizing neural activity, preserving the temporal structure by treating each spike as a discrete token. This strategy allows training across multiple sessions and individuals. The success of the approach was extensively tested in several datasets across multiple labs, brain regions, and tasks, demonstrating increased decoding performance compared to alternative methods.", "questions": "I would be interested in applying this method to even more diverse recordings throughout the brain and tasks to explore universal computational principles.\nThe models allows for neural variability, but still assumes that all neurons are functionally similar. Could this be a problem when looking at other brain regions? Specially if the proportion of excitatory cells is smaller.\n", "limitations": "The authors mention some limitations and future work in their paper. However, it would be important to acknowledge the computational costs, training times, and data demands associated with the proposed method. ", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "The paper is presented clearly and is technically sound. The method was tested and shown to work well across multiple datasets. The proposed strategy leverages the success of transformer models and applies it to neuroscience by introducing a novel way to tokenize neural activity. Extracting shared variability across multiple datasets is crucial for neuroscience and biomedical applications, and this approach could provide a framework for analyzing these emerging datasets. The authors demonstrated the success of the approach when trained on multiple neural recordings and tasks, with the ability to pool data across sessions and animals. Decoding results showed that the method outperforms alternative approaches, including lesion versions of the introduced model, when tested within and across sessions and animals. They also investigated the impact of different architectures and numbers of parameters on performance. Importantly, they showed that minimal fine-tuning allows the model to improve decoding performance on other tasks. The proposed strategy could lay the foundation for studying unifying principles of neural computation.", "weaknesses": "However, it should be noted that this strategy heavily relies on having access to multiple large neural recordings, which may not be feasible for most basic neuroscience research where datasets are often limited. While the authors tested the method across multiple brain regions and tasks, all of them were motor tasks in motor-related areas. While this choice likely facilitated the identification of shared structure and enabled transfer learning, it would be even more intriguing if this approach could be applied across categorically different tasks, such as sensory encoding and decision-making, spanning from sensory areas to higher cortical regions. This could potentially help identify universal coding principles of neural computation as well as task-specific coding principles related to motor function.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review of the paper! We are very excited to hear that you found the paper to be “presented clearly and is technically sound” and also agree that “extracting shared variability across multiple datasets is crucial for neuroscience”.\n\nIn what follows, we will provide point-by-point replies to your questions.\n\n> 1. “it should be noted that this strategy heavily relies on having access to multiple large neural recordings, which may not be feasible for most basic neuroscience research where datasets are often limited.”\n\n**Reply:** We want to clarify two points regarding POYO’s ability to help in settings where data is limited. First, we have a single-session variant of our model (POYO-[single session]) which is not pre-trained on any data but is randomly initialized and trained end-to-end on a single session worth of data. This model, even without access to additional recordings for pretraining, outperforms other single-session methods (Table 2). In addition to working well on single sessions, the method is very sample efficient: In Figure 4, we show that our model can be trained from scratch from as little as 8 trials worth of data (which is less than 2 minutes), and is more data-efficient than a GRU model. \n\nSecond, we agree that neural datasets are usually limited; this is a core motivation for our work and the creation of a large-scale pretrained model. This work shows that it is possible to combine multiple heterogeneous datasets that are collected in different labs, for different behavioral tasks, and under different experimental conditions, to provide a pretrained model that can transfer well on new sessions (Figure 3). Hence, we advocate for leveraging such an approach to make best out of available datasets. For example, a smaller lab with limited access to data could take a POYO model trained on much larger, open datasets (such as those available on DANDI) and then fine-tune it on their smaller dataset. We believe that expanding our method to decode various behavioral signals within smaller datasets will be an important objective for future work in order to demonstrate this potential use of POYO.\n\n> 2. “I would be interested in applying this method to even more diverse recordings throughout the brain and tasks to explore universal computational principles. The model allows for neural variability, but still assumes that all neurons are functionally similar.”\n\n**Reply:** The question of across region generalization and applying the model to larger-scale across-region questions is an exciting open area for future work! Even now in POYO-1, there is a lot of diversity in the location and areas for units in the different recordings, which span M1, PMd, and S1 (see Figure 1 in the attached pdf). These regions are long known to have distinct firing properties, cell types, and underlying dynamics. Yet, we are able to capture and exploit the activity of these diverse cells in a common model. We think POYO will be an extremely valuable tool in future work to better understand commonalities and differences in the neural code in different regions, both by testing generalization across regions and by studying the properties of models jointly trained on many regions.\n\n> 3. Could this be a problem when looking at other brain regions? Especially if the proportion of excitatory cells is smaller.\n\n**Reply:** In preliminary experiments, we have started to test the model on datasets from the visual cortex with good success. We plan to expand the model and training datasets in a next publication. In principle, we don’t believe that the number of excitatory cells will be a limiting factor, as long as the model is pre-trained on a diverse set of brain regions with different numbers of excitatory cells.\n\n\n> 4. The authors mention some limitations and future work in their paper. However, it would be important to acknowledge the computational costs, training times, and data demands associated with the proposed method.\n\n**Reply:**  Thank you for bringing this up. We discuss some of the computational requirements of our method in Appendix A.4.1, but we plan to include a paragraph in the discussion on this topic. Currently, we show that our largest model, trained on billions of tokens, uses 8 GPUs over a few days. As we increase the number of datasets and the complexity of the model, we expect that computational requirements will increase, but we are optimistic about recent advances in model scaling that have been enabled by breakthroughs in NLP. In terms of data demands, we have shown that large recordings are not a requirement to achieve performance that is competitive with existing methods, but our results suggest that we get significant boosts when a lot of data is leveraged. In settings where a lot of data has not been collected yet, we think that an important question for future work is whether we can leverage learned patterns of the neural code across various brain regions.\n\nThat being said, we see the ability to train on a large-scale model, and use it in new settings via. unit identification or finetuning to be the key objective here, as the large model can be trained once but leveraged many more times. To give a concrete example, we plan to publicly release POYO-1 which can be tuned by anyone that has recordings from motor areas of non-human primates. Since the model is pre-trained, it can be tuned with a few samples (sample-efficient, Fig 4), it converges in a small number of steps (training time, Fig 4), it can be done with a single GPU or even just a CPU in less than an hour (Appendix A.4.1), and it does not require extensive hyperparameter tuning (all our transfer results in Table 2 were obtained using the same hyperparameters despite major differences across datasets), making our model accessible. We will highlight these advantages in the manuscript.\n"}, {"review_id": "8MRRxJ1Lhf", "reviewer": "Reviewer_Mfux", "summary": "Neural population decoding refers to inferring the behavioral output of organisms from a recording of a population of neuronal recordings. This paper introduces a transformer architecture to improve the accuracy of decoding. The efficiency of the system is increased by representing input spike trains with event-based tokens and by switching to a latent representation (to decrease the size of attention matrices).\n\nIn the presence of behavioral data, the model can identify new neurons (units) that are potentially obtained from different experiments. This is done by freezing the parameters of the network during gradient descent except for the encoding of the new units. This is a novel and interesting approach although it depends on the presence of behavioral recordings.\n\nThe method is tested on multiple experiments where neural activity recordings are obtained from motor cortical regions of monkeys while the monkeys are engaged in tasks involving arm movements.", "questions": "- What is the relationship between the context window T (1s) and Tmax (4s) used in positional encoding? How does Tmax change with respect to T?\n\n- Could you provide more details on extending the model to self-supervised tasks (mentioned at the end of the paper)?", "limitations": "Yes.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- A novel framework for neural population decoding\n- Efficient representation of the recordings via event-based tokenization\n- Identification of units from new experiments in a way that is consistent with previous experiments. i.e., embedding a new unit such that the embedding location is consistent with the locations of other units in previous experiments whose activities/roles are similar to the new unit.\n- Use of multiple experiments done in different labs to demonstrate the method.\n- The inference accuracy of the method seems to improve significantly upon strong baselines.", "weaknesses": "- I oppose alluding to 'foundation models', which have become popular in natural language processing. Importantly, such models are trained in an unsupervised manner and demonstrate emergent abilities to solve multiple downstream tasks. The proposed architecture depends strongly on the presence of joint behavioral recordings. Similarly, it wouldn't make sense to apply this trained model to recordings from, say, the visual cortex.\n\n- [line 165] Are commands given via natural language? If not, please provide mathematical description instead (or in addition).\n\n- One simple experiment to probe the fidelity of \"unit identification\" could be: (i) take a recording, learn embeddings of units, etc. (ii) shuffle the order of units and repeat the experiment by adding rows to Embed as described in text. (iii) Compare and report the similarity between the two embeddings.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your feedback and questions! We are very excited to hear that you found the paper to be “novel” and also appreciate our use of “multiple experiments done in different labs”.\n\nIn what follows, we will provide point-by-point replies to your questions.\n\n> 1. “I oppose alluding to 'foundation models', which have become popular in natural language processing. Importantly, such models are trained in an unsupervised manner and demonstrate emergent abilities to solve multiple downstream tasks. The proposed architecture depends strongly on the presence of joint behavioral recordings. Similarly, it wouldn't make sense to apply this trained model to recordings from, say, the visual cortex.”\n\n**Reply:**  We agree with the reviewer that our current model wouldn’t be considered a foundation model, and we will clarify that in the text. However, POYO makes important advances towards a foundation model for neuroscience on multiple fronts.\n\n1. Scalability:  Building foundation models for neural data requires training on large amounts of data, and a key innovation of our work is developing a framework where integration and training on many sessions is made possible. \n\n2. Pre-training and transfer learning:  A key advantage of large language and other foundation models is that they enable fine tuning on smaller datasets. In our experiments, we show how our pretrained models POYO-MP and POYO-1 can be used on a number of diverse recordings from non-human primates, even if the data is from a new animal, their behavioral task is different (Figure 2), and when the preprocessing of the data is different (POYO-MP to NLB-Maze). In our few-shot learning experiments, we show that finetuning can quickly integrate new sessions from different labs with fewer than 32 labeled trials comprising less than 2 minutes of recordings.\n\n3. Democratizing AI:   In addition to providing pretrained models and transfer strategies, POYO also can be very easily adapted and trained without extensive hyperparameter tuning (which usually requires expertise), and without need for large compute resources. Unit identification can be done on CPU in a few minutes, or finetuning can be done on a single GPU. We note that the same hyper parameters are used across all our results in table 2, showing its robustness.\n\nMoreover, as noted in the paper, there is no reason that the POYO tokenization scheme and architecture could not be used for self-supervised learning. We will mention this in the text when we raise the issue of foundation models.\n\n> 2. “Are commands given via natural language? If not, please provide a mathematical description instead (or in addition).”\n\n**Reply:** The command is not given via natural language, we will add this for clarity. We will update the text to include a clear mathematical formulation, to avoid any confusion:\nA query token is defined by its embedding $y_{0, i} = Embed(session\\ id) + Embed(behavior\\ id) $ and its timestamp $t_i$.\n\n> 3. One simple experiment to probe the fidelity of \"unit identification\" could be: (i) take a recording, learn embeddings of units, etc. (ii) shuffle the order of units and repeat the experiment by adding rows to Embed as described in text. (iii) Compare and report the similarity between the two embeddings.\n\n**Reply:**  Thank you for suggesting this experiment! We provide an in depth discussion of the results in the General Response (see Table 2). Further visualizations are provided in the attached PDF (See Figure 1).\n\nAs you can see, when we tested your idea using POYO-1 weights, we indeed found that the tuned embeddings converge to their original unit embeddings (in a nearest-neighbor sense). We are excited to refine this analysis and include these results in the final paper. Thanks again for your suggestion!\n\n\n> 4.  What is the relationship between the context window T (1s) and Tmax (4s) used in positional encoding? How does Tmax change with respect to T?\n\n**Reply:** This is a great question that we explored in our initial experiments, where we did a sweep over various values of both Tmin and Tmax for POYO-SS (single session) and found that various values of Tmax provided good performance. We hypothesize that as long as Tmax >= T, our model should perform well. We will update appendix C.3 to include these insights.\n\n> 5. “Could you provide more details on extending the model to self-supervised tasks (mentioned at the end of the paper)?”\n\n**Reply:** The PerceiverIO architecture that we build upon is flexible and allows for different predictive targets, which could be adapted to other objectives like the prediction of the firing activity of held out and masked neurons as in the NDT.  Contrastive learning can also be used to generate various views (different sub-populations for example) and then learn an embedding that maps these views to similar parts of the latent space (i.e., using BYOL or SimCLR). We plan to add more discussion on these potential directions for future work in the revised manuscript."}]

Author rebuttal:
We would like to thank all of the reviewers for their great feedback and suggestions! The reviewers agreed on the impact of the work and acknowledged the innovations behind the work for multi-session neuroscience.

Some highlights and praise from the reviewers:
- **Method and Approach:** “the proposed technique is technically solid, elegant, and performs well in the experiments.” (Vmcn) “a novel way to tokenize neural activity” (BY6C) “Efficient representation of the recordings via event-based tokenization” (Mfux)
- **Impact:**  “The ability to combine various datasets in the way described here has the potential to be impactful.” (c1LQ) “this technique addresses a real need in neuroscience labs, making the paper potentially very high impact.”  (Vmcn) “a major problem in neuroscientific data analysis” (c1LQ) “Extracting shared variability across multiple datasets is crucial for neuroscience” (BY6C) “proposed strategy could lay the foundation for studying unifying principles of neural computation”
- **Scaling analysis and rigor of experiments:** “I found the scaling analysis (lines 237-257) particularly instructive” (Vmcn) “Use of multiple experiments done in different labs” (Mfux) “accuracy of the method seems to improve significantly upon strong baselines” (Mfux)
- **Writing and Presentation:** “the paper is very well written.” (Vmcn) “The paper is presented clearly and is technically sound.” (BY6C)

Based upon reviewer comments, we ran a number of new experiments, and are currently working on the following revisions to our original submission: 
- **Additional Baselines and Evaluations (c1LQ):** In response to rev (c1LQ)’s comments about including further single-session baselines, we conducted a number of new experiments to compare with the requested baselines (AutoLFADS, NDT [1]), in addition to a supervised-variant of NDT and another supervised neural data transformer baseline, EIT [2]. Across the board, we find that POYO-[Single-session] outperforms the other transformer-based approaches that use binning-based tokenization of the neural activity. Finetuning POYO-MP provides even further improvements in performance beyond the single-session models in the RTT task, demonstrating the power of having more data to pretrain on.

| | NLB-Maze | NLB-RTT |
| - | - | - |
| NDT + Linear | 0.8929 | 0.5895 |
| NDT-Supervised | 0.8708 | 0.4621 |
| EIT | 0.8791 | 0.4691 |
| AutoLFADS + Linear | 0.9062 | 0.5931 |
| POYO-[Single-session] | 0.9470 | 0.6850 |
| POYO-MP | 0.9466 | 0.7318 |
_Table 1: Behavioral decoding performance on NLB Datasets._

- **Neuron shuffling and identification experiment (Mfux):**  In response to reviewer Mfux’s great suggestion to conduct a neuron shuffling experiment, we ran unit-identification on a session that was already seen during pre-training. We test whether finetuning on this “new” set of units, will map the “new” units close to their true embedding (found during pre-training). We use our largest model, POYO-1, start with randomly initialized unit embeddings and run unit-identification. To compare the newly calibrated set of units with its pre-trained version, we normalize the unit-embeddings and report the cosine similarity averaged over the set of units. Note that normalization is done because POYO’s first cross-attention layer applies layer norm to these embeddings. We also report the “unit identification accuracy” which we define as the ratio of tuned unit embeddings that have their true unit embedding as their nearest neighbor.

| Dataset | # of units | Unit-Identification Accuracy | Unit-Identification Cosine-Similarity |
| - | - | - | - |
| Monkey C, CO, 2013/10/03 | 73 | 1.000 | 0.845 |
| Monkey M, CO, 2014/02/03 | 116 | 0.966 | 0.802 |
_Table 2: Unit re-identification results._

We include, in Figure 1 of the accompanying pdf, a visualization of the unit embeddings over multiple training steps during the unit-identification stage. We show how the new unit embeddings, which are initially random, converge towards their true embeddings, determined during pre-training.

Overall, this analysis suggests that the unit-identification approach is robust and reliable for identifying units. We believe that understanding the functional similarities of units that are mapped close to each other in the unit embedding space will be an interesting avenue for future work. 

**Impact of the work:** The proposed framework achieves a number of new firsts. POYO is the first to integrate and decode population activity across multiple animals using a unified, common model, and is also the first to jointly integrate datasets from different labs and behavioral tasks. The diversity of datasets that we have demonstrated our approach on is staggering: we have data from 9 nonhuman primates and over 150 sessions, train on diverse datasets from different spike sorting algorithms and threshold crossings, and across diverse behavioral tasks from 4 different research labs that are measured and executed through completely different manipulandum and at different sampling rates. 

---

References:

[1] Ye et al., "Representation learning for neural population activity with Neural Data Transformers." Neurons, Behavior, Data analysis, and Theory 2021.

[2] Liu et al., “Seeing the forest and the tree: Building representations of both individual and collective dynamics with transformers” NeurIPS 2022


Meta-review: This paper received unanimous support from all the reviewers. There is general agreement that the problem targeted is important and that the performance of the neural decoding algorithm is impressive. There were initial concerns regarding some of the claims and the experimental evaluation that appeared to have been well addressed by the rebuttal. The AC thus recommends the paper to be accepted.
Decision: Accept (poster)