Identifying interactions across brain areas while accounting for individual-neuron dynamics with a Transformer-based variational autoencoder
Authors: Qi Xin
Robert Kass
Venue: NeurIPS 2025 poster ()
Keywords: functional connectivity
neural spike train modeling
trial-to-trial variation
point process
generalized linear models
deep latent variable models
Primary area: neuroscience_and_cognitive_science
TL;DR: We introduce GLM-Transformer, a hybrid deep learning and GLM framework that identifies interpretable cross-area neural interactions while accounting for trial-varying individual-neuron dynamics via a Transformer-based VAE..

Abstract: Advances in large-scale recording technologies now enable simultaneous measurements from multiple brain areas, offering new opportunities to study signal transmission across interacting components of neural circuits. However, neural responses exhibit substantial trial-to-trial variability, often driven by unobserved factors such as subtle changes in animal behavior or internal states. To prevent evolving background dynamics from contaminating identification of functional coupling, we developed a hybrid neural spike train model, GLM-Transformer, that incorporates flexible, deep latent variable models into a point process generalized linear model (GLM) having an interpretable component for cross-population interactions. A Transformer-based variational autoencoder captures nonstationary individual-neuron dynamics that vary across trials, while standard nonparametric regression GLM coupling terms provide estimates of directed interactions between neural populations. We incorporate a low-rank structure on population-to-population coupling effects to improve scalability. Across synthetic datasets and mechanistic simulations, GLM-Transformer recovers known coupling structure and remains robust to shared background fluctuations. When applied to the Allen Institute Visual Coding dataset, it identifies feedforward pathways consistent with established visual hierarchies. This work offers a step toward improved identification of neural population interactions, and contributes to ongoing efforts aimed at achieving interpretable results while harvesting the benefits of deep learning.

Identifying interactions across brain areas while accounting for individual-neuron dynamics with a wransformer-based variational autoencoder
Qi Xin Carnegie Mellon University Pittsburgh, PA 15213 xinqi0511@gmail. com
Robert E. Kass Carnegie Mellon University Pittsburgh, PA 15213 kassestat. cmu, edu
Abstract
Advances in large-scale recording technologies now enable simultaneous mea- surements from multiple brain areas, offering new opportunities to study signal transmission across interacting components of neural circuits. However, neural responses exhibit substantial trial-to-trial variability, often driven by unobserved factors such as subtle changes in animal behavior internal states. To prevent evolving background dynamics from contaminating identification of functional coupling, developed hybrid neural spike train model, GI.--Tranformme that incorporates flexible, deep latent variable models into point process generalized linear model (GLM) having an interpretable component for cross-population in- teractions. A Transformer-based variational autoencoder captures nonstationary individual-neuron dynamics that vary across trials, while standard nonparamet- ric regression GLM coupling terms provide estimates of directed interactions between neural populations. We incorporate low-rank structure on population-to too- population coupling effects to improve scalability. Across synthetic datasets and mechanistic simulations, GI.M G..M--rrasffrrme recovers known coupling structure and remains robust to shared background fluctuations. When applied to the Allen Institute Visual Coding dataset, it identifies feedforward pathways consistent with established visual hierarchy. This work offers step toward improved identifica- tion of neural population interactions, and contributes to ongoing efforts aimed at achieving interpretable results while harvesting the benefits of deep learning.
Introduction
A central goal of systems neuroscience is to understand how multiple populations of neurons, across different brain areas, interact as a coordinated circuit to produce perception and behavior. Recent advances in electrophysiological recording technologies, which have enabled simultaneous recordings from thousands neurons across anatomically distinct parts of the brain, have offered exciting new opportunities for progress [] 2 4. Because neural action potentials, spikes, the dominant sources information transmission, much of modern neurophysiology has focused on sequences of spikes, that is, neural spike trains. Neural networks in living animals, however, also exhibit continually evolving dynamics. Thus, a major contemporary challenge is to develop methods for identifying interactions among spike trains from large populations of neurons while accounting for changing background dynamics.
One of the main collections of methodologies for analyzing neural spike trains has come under the heading of generalized linear models (GLMs), These models are based on Poisson nonparametric regression and often incorporate latent variables. They have been remarkably adept at identifying functional interactions between neurons due to their flexibility, interpretability, and rigorous statistical
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

In these models, a neuron's 's spiking activity is treated as point process whose intensity function depends on behaviorally relevant variables such as external stimuli, the recent spiking history of the neuron itself, and also the recent spiking history of other neurons. The latter is used to represent interaction effects as coupling terms in the intensity function.
foundations 6
10.
GLMs typically use a trial-invariant function of time, such as a function of time since stimulus onset, model nonstationary effects that are not explicitly attributed cross-neuron coupling or the neuron's 's own spike history [5
6
g 7
However, even under identical stimulus conditions,
[6]. This structured variability can be driven by behaviors such as locomotion, which induce strong and widespread background fluctuations in firing rates throughout the brain, including in primary sensory areas [17] [8] 19 20 21 22 23 24 [25 26, (27 (28]. Endogenous states such as arousal or attention can also contribute to trial-to-trial variability and are often difficult to measure precisely [12 25 28 Additional sources of within-area variability, possibly arising from unobserved factors, may further complicate the individual- neuron dynamies. Although scaling the firing rates of trial with trial-wise gain constant, including observed behavioral covariates in the model, may help mitigate effects due to trial-to-trial variability, they may be too simple or require careful hands-on model engineering.
individual-neuron dynamies often exhibit substantial trial-to-trial variability [12 13
Deep artificial neural networks, which can capture complicated nonlinear relationships 29 have become popular for modeling large-scale, multi-area neural activity [3] 32
33
34
35
36
37
38 39 40). They are particularly attractive for predictive purposes, but lack the straightforward interpretability of GL.Ms. In this article, we introduce hybrid framework we call GL.M Transforrormrrr
combines GL.M-based coupling to identify cross-area interactions, with a Transformer- variational autoencoder ( (VAH) to capture individual-neuron dynamics that do not require interpretation
(so these dynamics act as high- nigh-dimensional nuisance parameters, in the terminology of statistics). Following previous neural applications of deep learning, we incorporate dimensionality reduction
32 40 in the spirit of Gaussian Process Factor Analysis (GPFA) 41, where a small number of latent factors drive the activity of tens or hundreds of neurons. We demonstrate that our method can accurately recover coupling effects in synthetic datasets, generated from both statistical GLM neurons and mechanistic exponential integrate- and-fire (EIF) models, while alternative methods do not. We then apply our GLM- TrMM-rrsffrmmer to the Allen Institute's s Visual Coding dataset [42), where identifies strong directional coupling from the primary visual area V1 to the lateral medial (LM) and anterolateral (AL) visual areas, in agreement with prior studies. Together, these results illustrate the utility of combining deep artificial neural networks with interpretable statistical models in identifying cross-area interactions from large-scale neural recordings.
31
2 Related work
Some recent advances, including LFADS 31) NDT (33, NDT2 34 MtM 33 and Meta- Meta- Dynamicall State Space Models [43]. aim to learn meaningful latent representations at each time bin that generate observed activity. In contrast, Deep Random Splines (DRS) B77 aims to learn trial-level latent representations, which capture trial-to-trial variability. While these models offer powerful representations and predictions, they do not easily identify interactions across neural populations.
Several approaches based on state-space models have accommodated neural interactions, including mp-srSLDS [4), MR-SDS [40], while STNDT [39] treats neurons tokens in Transformer architecture. In these models, interactions occur in latent space, and in MR-SDS and STNDT, interactions are parameterized by deep networks. They again do not provide direct interpretation of the way activity in one population affects activity in another.
3 GLM-Transformer
The FLLM-Tranfformr decomposes drivers of each neural population's 's activity into three components: setividual-neuroo. dynamics, the population's 's self-history effects, and cross- cross-population coupling effects. Individual-neuron dynamics replace the baseline term or stimulus effect in the traditional GL. framework in order to capture not only standard experimental effects, but also trial-to-trial variability in background dynamics arising from other sources, including subtle changes in behavior and endogenous states. Because background fluctuations are often correlated across
Overview

---

A
B
Transfomer enooder, 2/11 deooder
Indivitual-reurron dynamics
Area
population
tral-wise latent (F1)
scli-history effects
time
nonlincarity
animal a. trial
low-rank coupling,
coupling
Area
fring rate À
(pop 2
populetion
low-rank coupling
crose-area coupling (pop
Area
population
Illustration of data and Toar-rrasformerr model architecture. (A) Neural activity is recorded across multiple brain areas, with data from each area organized as three-dimensional army of neurons, time, and trials. (B) Schematic of TLMM-ransformer For each trial and animal a, Transformer encoder processes spike trains from a target population to produce trial-wise latent variable z,,,, which is decoded into smooth individual-neuron dynamics. The predicted firing rate is computed by combining individual-neuron dynamics with self-history effects from the population itself and coupling terms from other populations, followed by nonlinearity.
Figure 1:
areas [5 19
28 explicitly modeling these shared dynamics helps prevent misattribution
of such variability to cross-area coupling and improves robustness of coupling identification. The self-history component accounts for recurrent dynamics within each population and each neuron's own post-spike effect (e.g., its refractory period). The cross- prossppopulation coupling component aims to capture dependencies between neural populations while accounting for the other components in the model.
GLM-Transformer supports training across multiple animals recording sessions. In animal a and neural population brain area let y E N denote the spike count of neuron n at time bin in trial r. Here, j 1, 2. P. where P is the total number of populations. We assume spike counts follow Poisson distribution, with the mean corresponding to the underlying firing rate at that time point: yre. Poisson( A,,,,)). Throughout the paper, we use a bin size of ms and rarely observe more than one spike per bin. To emphasize the smooth nature of firing rates over time, we use the notations x (t) and ,,, interchangeably.
Inspired by Deep Random Splines [37]. model individual-neuron dynamics with trial-wise
dz-dimensional latent variable z.v. E dRd that captures trial-to-trial variability. The log-intensity of neuron n in population is modeled as
log (1) sad (2
Σ
Crn (1)
where the individual-neure intensity function (2,",, .) depends on both time t and the trial-wise latent variable .; ; it is modeled basis functions whose coefficients are produced by Transformer-basee encoder-decoder.
liwdiimesioonal smooth function constructed from B- B-spline
The coupling term e4.2-t (t) captures the effect of spikes from population i on neuron n in population j through low-rank structure, implemented as a linear combination of convolutions between spike trains from population and learned temporal filters parameterized by raised cosine basis functions (6. Raised cosine basis functions concentrate density at short time lags and provide an inductive bias that better captures rapid post-spike and coupling effects compared with the smoother and more uniform B- B-spline basis.
The final term, heyd (t), represents the neuron's 's own post-spike effect and is computed by convolving the single-neuron spike train with post-spike filter. Together, h (f) and the within-populationn component of the coupling term (when = j) constitute the self-history effect for population j. To avoid non-identifiability.
penalize roughness in the background dynamics that they vary more slowly than the coupling effects. In the following sections, we describe each component of the model in detail.
Although GL. Transformer models multiple sources of variability, its primary goal is to identify directed interactions between neural populations. The trial-specific latent variables and individual-

---

neuron dynamics control for background fluctuations that could obscure these interactions, enabling more accurate and interpretable estimates of cross-area coupling.
Encoder The goal of the encoder network is to estimate the trial-wise latent variable z for population . based on the neural activity observed in that trial. Since z unobserved, we use a variational inference framework. We place standard Gaussian prior over the latent, z ~ (0, I), and upproximate its posterior with a Gaussian distribution whose mean and variance are predicted by Transformer-basee encoder, conditioned on the spike trains from the trial:
g (2..0)
Y F,n. N (Ho(YP))) ( diag(o (X.:))
where Yo is the full spike train matrix of population j in trial r, with neurons and time spanning the two axes; Mg and o are the outputs of the Transformer- encoder Φ, representing the mean and variance of the approximate posterior.
We use 10 or 20 ms time bins as tokens for the Transformer (which are wider than the 2 ms bins used elsewhere in the model). These bin widths were selected through hyperparameter tuning. which showed that performance was stable within this range and degraded slightly for smaller or larger bins. Because the set of recorded neurons differs across animals, we apply an animal-specific linear transformation to map the spike count vector in each time bin to shared dmodde-ddmmene -iimensiaal embedding space. To provide temporal context, we include elementwise summation with positional encodings. The resulting sequence is passed through stack of Transformer layers. Consistent with prior work, we found no benefit in using more than two Transformer layers or attention heads [33]. Finally, we obtain summary embedding via average pooling over all tokens, which is then passed through a linear layer to produce the mean vector and variance matrix (which is diagonal) of the approximate posterior on the trial-wise latent vector for each population.
The decoder uses the trial-wise latent z(-J to produce representation of individual -ecurnn dynamics. As in previous generalized factor analysis models [BI 32 40, 41), we assume that the individual-neuron dynamics of population are governed by set of underlying low-dimensional individual-neuron dynamics factors, fimila. (z,-) t), each varying continuously over time for a given trial. Here, "ind" indicates individual-neuron dynamics, and the index / = 1, 2...... I_ind enumerates the T,ind factors. These factors are B-spline functions, where the trial-wise latent variable 2, is mapped to the B-spline coefficients via a feedforward network (a multilayer perceptron, MLP) with one hidden layer. To encourage smoothness, we apply an & penalty on the second derivative of each factor function. In practice, we typically use 10 20 B-spline basis functions. The 1-th individual-neuron dynamics factor function is computed as
Decoder
f aind (t) FB-splines ([
where FMLP((Z):)) E L maps the trial-wise latent to a matrix of B-spline coefficients for all factors, and [-) denotes the extraction of the coefficients corresponding to the 1-th factor. F'B-spiines then applies these specific coefficients to the basis functions. To get the individual- dynamics intensity function for each neuron, we apply a linear readout from these factors:
Lind
f .)
Σ wind
gindiad (zf.) t),
(4)
where winnda.d are learnable readout weights that linearly combine the factors to produce the final individual-neuron dynamics intensity for neuron in population j.
Low-rank population-to-opplaaii coupling effeets To reduce the number of coupling param- eters, which could be quadratic in the number of neurons, impose low-rank structure the coupling effects between neural populations. Specifically, the influence of source population on target population j is mediated through set of now-dimensiona coupling factors:
No.1
a.mplinng..i+))
Σ
aitt in
(g a.iit x T.
xya,
(t),
(5)

---

where No. is the number of neurons in source population i for animal and - 1, 1...... 2. coupling indexes the c,coupling coupling factors. The temporal filter o. is learnable and parameterized using smooth raised cosine basis functions [6), and the convolution * is applied the spike train r. of neuron n. The sending weight ai,n ai modulates the contribution of neuron n to the 1-th coupling factor. We have found it sufficient, in practice, to use only one factor for each coupling (i.e.,
= 1). The coupling effect on the target population is then
I. Lonepline
(t)=
W
Wroopling...-59))
f fcoupling...--1
(6)
where Woopplinga is the receiving weight that maps coupling factor / to neuron in the target population j. We add a small & penalty (1 x 10-9) on both the sending and receiving weights during training.
Self-history effects Self-history effects are modeled by combining within cottin-population coupling and neuron-specific post-spike filters. The within- aithin-population coupling from population i to itself handled using the same low-rank structure described for cross-populte poosslpopulttion interactions (i + i). Additionally, each neuron has its own post-spike filter that captures history- -ependent modulation, such as refractory periods. These post-spike filters are also parameterized using a set of smooth raised cosine basis in time delay window of 10 ms 6 45 train.
and are convolved with the neuron's 's own spike
Training Each dataset is split into training (70%), test (10%), and validation (20%) sets. We used a four-stage training procedure that progressively incorporates model components: (1) train only the trial-invariant
selvidual-neuron dynamics, without other components, using maximum likelihood; (2) enable the VAE and train the trial-varying individual- dynamics with evidence lower bound (ELBO); (3) include population coupling effects and continue training with ELBO; (4) include neuron-specifte post-spike filters and fine-tune the full model. This staged approach improves convergence and avoids poor local optima (compared with training all components together). We used Bayesian optimization for hyperparameter tuning. Full details about training and hyperparameter tuning are provided in Appendix E All experiments were run on single NVIDIA GeForce RTX 4090 GPU (24 GB) using PyTorch 2 with CUDA 12.1. Training on the Allen Institute dataset from ten animals took approximately 15-20 hours.
Identifiability If the individual-neuron dynamics (z,))) .) are too flexible, the model can become non-identifiable: the dynamic components may absorb variation that should be attributed to coupling effects. To avoid this, we made a strong assumption that coupling effects vary more rapidly than the background fluctuations arising from changes in behavior or endogenous states. This assumption is motivated by empirical observations that slower fluctuations in neural activity tend to have higher spectral power [12 and often driven by behavioral and arousal-related factors evolving over hundreds of milliseconds, whereas cross-area interactions typically occur at much shorter latencies [46). This is why we regularized the individual- dynamics and used low- dimensional triel-wise latent variable .,, Also. because cross-area coupling often involves only a small subset of neurons within a population [47 46 48), the loading patterns of individual-neuron dynamics factors and coupling factors typically differ, which further helps to disentangle these two components. When GLM- Transformer is reconfigured using a large number of B-spline basis functions, no smoothness penalty, a nigh-dimensiioal trial-wise latent variable "バ- and an increased number of factors for the individual- neuron dynamics term, the troublesome absorption of coupling effects by individual-neuro dynamics can occur. See Supplementary Figure ST
Summary of model components To clarify the methodological contribution, we summarize how
integrates and extends prior approaches. The individual- dynamics component builds on ideas from Deep Random Splines (DRS) [37 and Neural Data Transformer (NDT) (33. It combines the latent-per-tria structure and spline-based decoder from DRS with the flexible Transformer encoder from NDT, allowing trial-specific background dynamics to be modeled with both flexibility and structure. The low-rank GL.--tyll coupling component represents new integration: although inspired by reduced-rank regression, it extends this idea point process GL.Ms for estimating directed cross-area interactions in an interpretable way. The low-rank structure is
5

---

Individual-neuron dynamics loading matrik
Individual-neuren dynamics factor
Individual-neuron dynamics factor
factor factor
Sorted meeund Index
Fitted
0.0
0.0
0..5-
250 Time (ms)
500
250
500
Time (ms)
Weight
Weight
Coupling loading matrix sending
Coupling Filter
Coupling Factor
receiviing
1.0
Softed nouron
Sonted neuron
0.5-
0.0
Time
Weight
Weight
Figure Recovery of ground truth components in synthetic data generated from GLM spiking neurons. Gray: ground truth; blue: fitted. Top row: ground truth versus fitted individual- dynamics factors of an example trial (left, middle) and their loading matrix (right). Bottom row: ground truth versus fitted coupling filter (left), coupling factor of example trial (middle), and sending/receiving weights of the coupling term (right).
crucial it reduces model complexity and enable neural populations while preserving the interpretability coupling estimates. The post-spike history term follows established GLM formulations. Taken together, these elements form unified framework that combines the interpretability of GLMs with the representational power of deep latent models, providing a principled and scalable approach for disentangling background dynamics from true cross-area coupling.
scale efficiently large
4 Experiments and results
4.1 Ablation study
We provide an ablation study (Table in the Appendix that evaluates the contributions of key components. These results confirm that each design choice improves model performance on both synthetic and real data.
4.2 Validation on synthetic data from GLM neuron models
We first evaluated whether GLM- Transformer can recover the ground truth structure when data are generated from a known GLM model. We simulated spiking activity from two neural populations, each containing 50 neurons. The individual-neure
dynamics (i.e., baseline firing rates) were generated using two latent factors per population, drawn from Gaussian processes with temporal correlation structure defined by 300-ms time constant. These factors were linearly combined using different loading weights each neuron, mimicking structured trial-to-trial variability from unobserved sources. Ccros-poppulation coupling was introduced from the first population to the second using the same low-rank structure assumed by GLM-Transfoon Transformer. Nonzero sending and receiving weights were assigned only the first 10 neurons in each population, reflecting experimental findings that only a subset of neurons participate cross-area interactions 47 46 48), Spike trains were generated via inhomogeneous Poisson process, with conditional log- inggintensity given by the sum of the individual-neuroo dynamics, cross- coossloopulttion coupling, and a neuron- -ppenncc post-spike history term. The data set contained 2,000 trials, each 500 ms long, with average firing rate of 26 spikes/sec.
Fitting TrL--frraffrmmrr to this synthetic dataset, it successfully recovered the ground truth. As shown in Figure?? 2 the estimated individual-neuron dynamics and, more importantly, the cross- population coupling match the true underlying components. The fitted post-spike filters also align

---

B
Source area population
Exteraal Current (nv)
No coupling
Weak coupling
Strong caupling
Indiviuale dyuran Hogg relee)
input
O
0.2
O0
250
500
250
500
250
500
Time (ms)
Time tms)
Time (ms)
Synaptic Current (nt)
Coupling Blgg firing rate)
o
Mapay
Target area population
0.0
0.0
250
500
250
500
250
500
tril3
Time (ms)
Time (ms)
Time (ms)
Figure 3: Comparison of ground truth synaptic currents with fitted components from GLM- Transformer in EIF neuron simulations. (A) Schematic of the simulation setup: external input drives both source and target populations, with synaptic connections from source target neurons. (B) Each column shows results from one example trial under a different coupling condition: no coupling (left), weak coupling (middle), and strong coupling (right). Top external input (gray, left axis) and fitted individual-neurgy dynamics (blue, right axis), averaged across all neurons in the target population. Bottom row: ground truth synaptic input from source target (gray, left axis) and fitted coupling effect (blue, right axis), averaged over the 10 target neurons receiving cross poossppopulttoon input. Fitted coupling effects closely match the true synaptic currents across conditions.
with the ground truth (Supplementary Figure S2 These results validate both the accuracy of the T.MM-Tranffrmmrr and the training procedure.
4.3 Validation simulated data from mechanistic neuron models
To test the robustness of ELM-Transformer
more realistic data, simulated spike trains from two
interconnected neural populations using an exponential integrate-and-fin (EIF) model 449 50 51]
(see Appendix [F for details). External inputs had dual-peaked temporal shape, mimicking responses to drifting gratings in mouse visual cortex from the Allen Institute Visual Coding dataset, with peaks around 50 and 250 ms after stimulus onset (46) Trial-varying gain modulations, sampled from a Gaussian process, were added to mimic behavioral fluctuations. Each population contained 50 neurons, with excitatory synapses from 10 neurons in the first population to 10 in the second, and log-no log-normally distributed synaptic strengths [52]. Recurrent connections were also included within each population. We varied the mean cross-population synaptic strength to assess whether could recover the underlying coupling.
Figure shows that Transfororr produces interpretable components that qualitatively match the underlying biophysical ground truth across varying coupling strengths. Although the fitted individual-neuron dynamics do not exactly match the external input due to enforced smoothness, they capture the main temporal features, and more importantly the fitted coupling effect matches the ground truth synaptic input in all three conditions. To further evaluate scalability, simulated is still able
larger network with four interacting EIF populations and we found GLM-Transfoo recover the directed connections (Supplementary Figure SS}. These results demonstrate the reliability of the which is that the fitted components can reflect underlying currents or voltage signals 10.2.1].
while highlighting an advantage of GLM-based approaches,
[50] Section
4.4 Robustness to shared background dynamics
key motivation for ILMM-raanforrmrr is to avoid spurious coupling estimates in the presence of background gain fluctuations that are shared across neural populations. To evaluate this, we simulated two EIF neuron populations as in the previous experiment, but introduced correlated nonstationary gain modulation functions to the two populations. Importantly, no synaptic coupling was present between them. All neurons were simulated using the same EIF model described earlier. Hyperparameters for the GLM- TMM-Tranffrmer were also the same in the previous figure.
7

---

External input with share background
Source area population
B
O
O
Contribuziin of couplinat00
Tranfommrr
GLM Transfomar trialvarying baseline
OO
Reduced regression
Singla-neuron GLM
20-
Single-neuron trial tribl-spciffc gain
mDLAG
00
0.0
0.2
0.4
Correlation background dynamics
tial
Target area papulation
Figure 4: Robustness to shared background dynamics. (A) Schematic of the simulation setup: two neural populations receive correlated external input but have direct synaptic connections between them. (B) Percentage of variance in predicted firing rates attributed to coupling effects, as a function of background gain correlation. Although true coupling exists, most baseline methods spuriously attribute correlated inputs to coupling. In contrast, GLM- Transformer maintains low estimates across conditions, due to its trrainvarying individual-neurve dynamics component. The x-axis represents correlation in the underlying Gaussian process gain modulation, rather than the total input current.
We compared GLM- Transformer several existing methods, including reduced-rane regression (RRR) [48], delayed latents across multiple groups (mDLAG) [47], and single-neuron GLMs (with and without trial-wise constant gain modulation) [6] we also included variant of GLM-Transforme without trial-varying individual-neuron dynamics. Figure shows the proportion of variance in predicted firing rates attributed to coupling effects as the correlation in shared nonstationary gain increases. While all other methods incorrectly attributed shared nonstationary gain to coupling, the TLM--ransformer consistently maintained low estimates of spurious coupling across all conditions. Notably, the trial-varying individual-neuron dynamics component plays critical role, removing it from our model (see the orange curve Figure leads to overestimated coupling effects, similar to the other methods.
4.5 Identifying cross-area interactions in the Allen Institute Visual Coding dataset
We applied Tran-foafforrmr to the Allen Institute Visual Coding dataset, which contains simul- taneous Neuropixels recordings from six visual areas in awake mice, including the primary visual cortex (VI), the lateral medial area (LM), and the anterolateral area (AL). During each trial, mice passively viewed visual stimuli including Gabors, drifting gratings, and natural movies [42). The model was trained on available data from ten animals, comprising 121,679 trials of 400 ms each.
After hyperparameter tuning, used two factors for individual-neuron dynamics, one factor for each coupling term, and chose d, (see Appendix for details).
TLM--ransformer identified strong coupling from VI to LM and from VI to AL (Supplementary Figure S4). consistent with the known visual hierarchy 28 and prior findings from the same dataset [46] Figure |3 ppeeeenns detailed results for LM's individual-neuron dynamics and the coupling from V1 LM in one animal previously analyzed by Chen al. 46). As shown in Figure [5A, LM's individual-neuron dynamics exhibit substantial trial-to-trial variability. Notably, the neurons with the strongest coupling weights, both in VI (sending) and in LM (receiving), substantially overlap with the "cross-pop" neurons identified by Chen et al. 44 (Figure DD, neurons in red). Chen et al. showed that this subset of neurons plays a central role in cross-area interactions, as the population firing peak times in VI and LM were highly correlated when computed using only these neurons. In contrast, we observe less alignment between these "cross-pop" neurons and the loading weights of LM's individual- dynamics factors (Figure 3)), suggesting that the two terms capture different patterns.
To further examine the trial-wise latent . we visualized LM's latents using t-SNE (Figure 5). The embedding shows clear organization with respect to locomotion and pupil diameter, a proxy for arousal. Since these behavioral variables are major sources of background dynamics across the brain [13
15 [6), this alignment suggests that the inferred trial-wise latent captures meaningful trial-to-trial variability.
8

---

individual-nnuron term factor 2
B
Loading matrix factor
2? of LM
factor
factor 2
Safted meurrn indee.
-0.5
Seroed neuron
10-
-1.0
trial
-1.5-
trfal
-3.0
brial
200
400
200
400
0.0
Time (ms)
Time (ms)
Weight
Weight
V1
coupling
Loading matrix
coupling filter
coupling Tactor
sending
receiving
Sorted index
-10
-5
0
5
10
Locomotion
Pupil diameter
Stationary
Min pupil
200
0.5
0.0 0.5
Running
Max pupil
Time lag (ms)
Time (ms]
Weight
Weight
Figure 5: Individual-neuron dynamics and coupling effects inferred by TLM-Transformer in the Allen Institute Visual Coding dataset. (A) Trial-specific individual-neuron dynamics in LM captured by two smooth latent factors, shown across three example trials. (B) Loading weights from the individual-nenee previously identified by Chen et al. [46] as participating in cross-area coordination. (C) Coupling filter (left) from VI LM and the corresponding coupling factor for three example trials (right). We only have one coupling factor from an area to another. (D) Sending and receiving weights of the coupling factor from VI to LM. Red squares again indicate the "cross-pop" neurons, which match the neurons with highest sending/receiving weight in G..--Tranfformer. (E) t-SNE embedding of the trial-wise latent variable z of LM colored by behavioral state (locomotion: stationary vs. running) and sized by pupil diameter. The learned latents capture both locomotion and pupil diameter effectively.
dynamics factors to LM neurons. Red squares mark "cross-pop" neurons
The alignment between "cross-pop" neurons and those with large coupling weights, along with the organization of trial-wise latents by locomotion and arousal, provides indirect evidence for the reliability of HLMM-ransformert two key components: coupling and individual- dynamics. These results suggest that Traa-rrrsorrm can jointly capture interpretable cross-area interactions and rich trial-to-trial variability.
5 Discussion
The SILM--Trasffrrme combines GLM-style coupling components with flexible Transformer model for within-area dynamics. thereby enables interpretable estimation of directed interactions across populations of spike trains in the presence of trrainvarying dynamics, as observed in experimental data [13 14 15 16). The method is also scalable in practice, supporting joint training across multiple animals, sessions, and datasets, which can improve generalization and robustness.
Despite its advantages, the GLM- Transformer has several limitations that could be explored in future work. One issue is the dependence of the Transformer training. Although our multi-stage training procedure mitigates some of the variability across alternative settings, and we observe consistent performance. results across broad range of hyperparameters, there are no theoretical guarantees easy to imagine substantial effort devoted to training under wide variety of cireumstances based on the ever-increasing availability of multiple spike train data. In addition, CL.M-Transformer more data demanding than traditional GLMs, typically requiring over a thousand trials to converge, which is about order of magnitude more data than standard GLMs. This limitation could be mitigated by initializing the model from a pretrained network when only a small dataset is available.
We have also made several simplifying assumptions. For example, our implementation assumes fixed trial length, which limits applicability datasets with variable temporal structures. practical workaround for datasets with a few discrete trial lengths would be to learn separate decoders per trial type. A more general extension would involve replacing the spline-based decoder with a GPT-style architecture capable of handling variable-length sequences, which could broaden the scope of tasks addressable by .rTanssoorrmrr.. addition, we currently use the same number of factors for the individual-
neiuro neuron and coupling components, choice that could be relaxed in future work.

---

Finally, a key scientific challenge is to develop methods for statistical inference about effects of interest models like GLM--rrasformes where repeated retraining would be computationally expensive. Recent research on post-hoc debiasing and uncertainty-aware estimation 53 suggests these strategies might be made applicable to GiaM-Tranfformer. Once rigorous inference about cross-area interaction effects becomes available, models following the general approach deployed GLM-Transformer could become useful many studies involving large-scale neural recordings.
54
55
Code availability
GitHub repository with code to reproduce all experiments is available at: Qi TainGGM--Transforr
com/
References
[1] Joshua Siegle, Xiaoxuan Jía, Séverine Durand, Sam Gale, Corbett Bennett, Nile Graddis, Greggory Heller, Tamina K Ramirez, Hannah Choi, Jennifer A Luviano, et al. Survey of spiking the mouse visual system reveals functional hierarchy. Nature, 592((852):86 92, 2021.
[2] James Jun, Nicholas A Steinmetz, Joshua Siegle, Daniel J Denman, Marius Bauza, Brian Barbarits, Albert K Lee, Costas A Anastassiou, Alexandru Andrei, Çağatay Aydın, al. Fully integrated silicon probes for high-density recording neural activity. Nature, 551(7679):232- 236, 2017
[3] Nicholas Steinmetz, Cagatay Aydin, Anna Lebedeva, Michael Okun, Marius Pachitariu, Marius Bauza, Maxime Beau, Jai Bhagat, Claudia Böhm, Martijn Broux, et al. Neuropix- els 2.0: A miniaturized high-density probe for stable, long-term brain recordings. Science, 3e266599) 2021.
[4] Timothy A Machado, Isaac V Kauvar, and Karl Deisseroth. Multiregion neuronal activity: the forest and the trees. Nature Reviews Neuroscience, 23(1 1):683-7444 2022.
[5] Wilson Truccolo, Uri T Eden, Matthew R Fellows, John P Donoghue, and Emery N Brown. A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects. Journal of neurophysiology, 93(2):1074-1089) 2005.
[6] Jonathan W Pillow, Jonathon Shlens, Liam Paninski, Alexander Sher, Alan M Litke, EJ Chichilnisky, and Eero P Simoncelli. Spatio-temporal correlations and visual signalling in complete neuronal population. Nature, 454(7207):995 9999, 2008
[7] Jacob L Yates, Il Memming Park, Leor N Katz, Jonathan W Pillow, and Alexander Huk. Func- tional dissection of signal and noise in mt and lip during decision-making Nature neuroscience, 20(9): 1292, 2017.
[8] Eric Hart and Alexander C Huk. Recurrent circuit dynamics underlie persistent activity in the macaque frontoparietal network. Elife, 9:e52460, 2020.
[9] Stephen L Keeley, David M Zoltowski, Mikio C Aoi, and Jonathan Pillow. Modeling statistical dependencies in multi-region spike train data. Current Opinion in Neurobiology, 65:194-202, 2020.
[10] Adam Kohn, Anna Jasper, João D Semedo, Evren Gokcen, Christian K Machens, and M Yu By- ron. Principles of corticocortical communication: proposed schemes and design considerations. Trends in Neurosciences, 43(9):725-737. 2020.
[11] Robert E Kass and Valérie Ventura. A spike-train probability model. Neural computation, 1713-177,), 2001.
[12] Robert E Kass, Heejong Bong, Motolani Olarinre, Qi Xin, and Konrad N Urban. Identifi- cation of interacting neural populations: methods and statistical considerations. Journal of Neurophysiology, 130(3):475445) 2023.
10

---

[13] International Brain Lab, Brandon Benson, Julius Benson, Daniel Birman, Niccolo Bonacchi, Matteo Carandini, Joana A Catarino, Gaelle A Chapuis, Anne Churchland, Yang Dan, et al. A brain-wide map of neural activity during complex behaviour. bioRxiv, pages 2023-07, 2023..
[14] Simon Musall, Matthew T Kaufman, Ashley L Juavinett, Steven Gluf, and Anne K Churchland. Single-trial neural dynamics are dominated by richly varied movements. Nature neuroscience, 22(10):1677-1686 2019.
[15] Carsen Stringer, Marius Pachitariu, Nicholas Steinmetz, Charu Bai Reddy, Matteo Carandini, and Kenneth D Harris. Spontaneous behaviors drive multidimensional, brainwide activity. Science, 364(6443))caavv7811 2019.
[16] Dario L Ringach. Spontaneous and driven cortical activity: implications for computation. Current opinion in neurobiology, 19(4):439-444 2009.
[17] Cristopher M Niell and Michael P Stryker. Modulation of visual responses by behavioral state in mouse visual cortex. Neuron, 65(4):472-479, 2010.
[18] Georg B Keller, Tobias Bonhoeffer, and Mark Hübener. Sensorimotor mismatch signals in primary visual cortex of the behaving mouse. Neuron, 74(5):809-815, 2012.
[19] Jérôme Lecoq, Joan Savall, Dejan Vučinić, Benjamin F Grewe, Hyun Kim, Jin Zhong Li, Lacey J Kitch, and Mark Schnitzer. Visualizing mammalian brain area interactions by dual-axis two-photon calcium imaging. Nature neuroscience, 17(1))1825-1829, 2014.
[20] Mark Andermann, Aaron M Kerlin, Demetris K Roumis, Lindsey L Glickfeld, and Clay Reid. Functional specialization of mouse higher visual cortical areas. Neuron, 1225)1102-- 2011.
[21] Ash Ayaz, Aman B Saleem, Marieke Schölvinck, and Matteo Carandini. Locomotion controls spatial integration in mouse visual cortex. Current Biology, 23(10):890-89) 894, 2013.
[22] Corbett Bennett, Sergio Arroyo, and Shaul Hestrin. Subthreshold mechanisms underlying state-dependent modulation of visual responses. Neuron, 30(2):350-337,200 2013.
[23] Pierre-Olivier Polack, Jonathan Friedman, and Peyman Golshani. Cellular mechanisms of brain state-dependent gain modulation visual cortex. Nature neuroscience, 16(9):133--139, 2013.
[24] Yu Fu, Jason M Tucciarone, J Sebastian Espinosa, Nengyin Sheng, Daniel P Darcy, Roger A Nicoll, Z Josh Huang. and Michael Stryker. cortical circuit for gain control by behavioral state. Cell, 156(6):1 (39-)) 152, 2014.
[25] Martin Vinck, Renata Batista-Brito Ulf Knoblich, and Jessica A Cardin. Arousal and locomotion make distinct contributions to cortical activity patterns and visual encoding. Neuron, 86(3):740-
754, 2015.
[26] Maria Dadarlat and Michael Stryker. Locomotion enhances neural encoding of visual stimuli in mouse vl. Journal of Neuroscience, 37(14)33764-37 3775, 2017.
[27] Amelia J. Christensen and Jonathan W. Pillow. Reduced neural activity but improved coding in rodent higher-order visual cortex during locomotion. Nature Communications, 13(1):1676, 2022.
[28) Lindsey L. Glickfeld and Shawn R Olsen. Higher-order areas of the mouse visual cortex. Anmual review of vision science, 3:251-273, 2017.
[29] Manuel Brenner, Elias Weber, Georgia Koppe, and Daniel Durstewitz. Learning interpretable hierarchical dynamical systems models from time series data. arXiv preprint 4rXiv:2410.08114, 2024.
[30] Matthieu Kirchmeyer, Yuan Yin, Jérémie Donå, Nicolas Baskiotis, Alain Rakotomamonjy, and Patrick Gallinari. Generalizing to new physical systems context-informed dynamics model. In International Conference on Machine Learning, pages 1123-11300 PMLR, 2022.
11

---

[31] Chethan Pandarinath, Daniel J O' Shea, Jasmine Collins, Rafal Jozefowicz, Sergey D Stavisky, Jonathan C Kao, Eric M Trautmann, Matthew Kaufman. Stephen Ryu, Leigh R Hochberg. et al. Inferring single-trial neural population dynamics using sequential auto-encoders. Nature methods, 15(1)):805 815, 2018.
[32] Andrew R Sedler and Chethan Pandarinath. Ifads-torch: A modular and extensible imple- mentation of latent factor analysis via dynamical systems. arXiv preprint arXiv:2309.01230 2023.
[33] Joel Ye and Chethan Pandarinath. Representation learning for neural population activity with neural data transformers. arXiv preprint arXiv:2 10::::22000.10 2021.
[34] Joel Ye, Jennifer Collinger, Leila Wehbe, and Robert Gaunt. Neural data transformer 2: multi- context pretraining for neural spiking activity. Advances in Neural Information Processing Systems,36:8035222 36:80322-80374, 2023.
[35] Mehdi Azabou, Vinam Arora, Venkataramana Ganesh, Ximeng Mao. Santosh Nachimuthu, Michael Mendelson, Blake Richards, Matthew Perich, Guillaume Lajoie, and Eva Dyer. A unified, scalable framework for neural population decoding. Advances in Neural Information Processing Systems, 36:44937-44956, 2023.
[36] Yizi Zhang, Yanchen Wang. Donato Jiménez-Benetó, Zixuan Wang, Mehdi Azabou, Blake Richards, Renee Tung, Olivier Winter, Eva Dyer, Liam Paninski, al. Towards a" universal translator" for neural dynamics at single-cell, single-spike resolution. Advances in Neural Information Processing Systems, 77880495-80521, 2024.
[37] Gabriel Loaiza-Ganem, Sean Perkins, Karen Schroeder, Mark Churchland, and John Cun- ningham. Deep random splines for point process intensity estimation of neural population data. Advances in Neural Information Processing Systems, 32, 2019.
[38] Ran Liu, Mehdi Azabou, Max Dabagia, Jingyun Xiao, and Eva Dyer. Seeing the forest and the tree: Building representations of both individual and collective dynamics with transformers. Advances in neural information processing systems, 35:2377-2391. 2022.
[39] Trung Le and Eli Shlizerman. Stndt: Modeling neural population activity with spatiotemporal transformers. Advances in Neural Information Processing Systems, 012117926-7999 2022.
[40] Orren Carniol-Tambour David M Zoltowski, E Mika Diamanti, Lucas Pinto, Carlos D Brody. David W Tank, and Jonathan W Pillow. Modeling detendent communication between brain regions with switching nonlinear dynamical systems. In The Twelfth International Conference on Learning Representations, 2024.
[41] Byron M Yu, John P Cunningham, Gopal Santhanam, Stephen Ryu, Krishna V Shenoy, and Maneesh Sahani. Gaussian prrcces factor analysis for low now-dimensional single-trial analysis of neural population activity. Advances in neural information processing systems, 21, 2008.
[42] Allen Institute MindScope Program. Allen brain observatory - neuropixels visual coding (dataset). 2019.
[43] Ayesha Vermani, Josue Nassar, Hyungju Jeon, Matthew Dowling, and Il Memming Park. Meta-dynamical state space models for integrative neural data analysis. arXiv preprint arXiv:2410.05454, 2024.
[44] Joshua Glaser, Matthew Whiteway, John P Cunningham, Liam Paninski, and Scott Linderman. Recurrent switching dynamical systems models for multiple interacting neural populations. Advances in neural information processing systems, 33: 33:14867-14878. 2020.
[45] Yu Chen, Qi Xin, Valérie Ventura, and Robert E Kass. Stability of point process spiking neuron models. Journal of computational neuroscience, 46(1): (661)) 2019.
[46] Yu Chen, Hannah Douglas, Bryan J Medina, Motolani Olarinre, Joshua H Siegle, and Robert E Kass. Population burst propagation across interacting areas of the brain. Journal of Neurophysi- ology, 128(6) 2022.
12

---

[47] Evren Gokcen, Anna Jasper, Alison Xu, Adam Kohn, Christian K Machens, and Byron M Yu. Uncovering motifs of concurrent signaling across multiple neuronal populations. Advances in Neural Information Processing Systems, 0e::4711-344722 2023.
[48] João D Semedo, Amin Zandvakili, Christian K Machens, M Yu Byron, and Adam Kohn. Cortical areas interact through a communication subspace. Neuron, 102(():249922)) 259, 2019.
[49] Nicolas Fourcaud- Troommééé David Hansel, Carl Van Vreeswijk, and Nicolas Brunel. How spike generation mechanisms determine the neuronal response to fluctuating inputs. Journal of neuroscience, 23(37): 1 1628-1 16640, 2003.
[50] Wulfram Gerstner, Werner M Kistler, Richard Naud, and Liam Paninski, Neuronal dynamics: From single neurons networks and models of cognition. Cambridge University Press, 2014.
[51] Chengcheng Huang. Douglas Ruff, Ryan Pyle, Robert Rosenbaum, Marlene Cohen, and Brent Doiron. Circuit models of low-dimensional shared variability in cortical networks. Neuron, 101(2):337-348, 2019.
[52] Sen Song, Per Jesper Sjöström, Markus Reigl, Sacha Nelson, and Dmitri B Chklovskii, Highly nonrandom features of synaptic connectivity in local cortical circuits. PLoS biology, 3(3):e68, 2005.
[53] Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins. Double/debiased machine learning for treatment and structural parameters, 2018.
[54] Edward H Kennedy.
Semiparametric doubly robust targeted double machine learning: review. Handbook of Statistical Methods for Precision Medicine, pages 207-236, 2024.
[55] Jakob Gawlikowski, Cedrique Rovile Njieutcheu Tassi, Mohsin Ali, Jongseok Lee, Matthias Humt, Jianxiang Feng, Anna Kruspe, Rudolph Triebel, Peter Jung, Ribana Roscher, et al. survey of uncertainty deep neural networks. Artificial Intelligence Review, 56(Suppl 1):1513-1589, 2023.
13

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
The abstract and introduction clearly describe the key contribution of our method, GLM-Transformer which combines GLM-style coupling with Transformer-basee variational autoencoder for modeling triainvaryingg individual andiiduull-neuron -nuron dynamics. The claims are supported by experiments on both synthetic data and real data in section Experiments and results.
Justification:
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We discuss the limitations of our method in the Discussion section.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g.. independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or speech-to-text system might not be used reliably provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't acknowledged
the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
14

---

Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
Justification:
We don't have theoretical result.
Guidelines:
The answer NA means that the paper does include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: We provide complete descriptions of the model architecture, training pro- cedure, datasets, and evaluation metrics the Methods, Experiments and Results, and Supplementary Materials. All experimental settings and hyperparameters, including ablation configurations and data splits, are described in detail to support reproducibility of our main results.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often good way to accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
15

---

(d)
We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We use publicly available data from the Allen Institute Visual Coding dataset, and provide GitHub repository containing all code needed reproduce our experi- ments.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. cc/
for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/ /nips.ccc/pullc//uuui
Cooe/obbeeeeeeeee wwtmonssssooni ssionPolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if' applicable).
Providing as much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: We provide complete details of the experimental setup, including data splits, training procedures, hyperparameter choices, and optimization settings. These are described in the Methods and Experiments and Results sections, with additional information in the Supplementary Materials.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
16

---

Justification: For all results that support the main claims (as opposed to illustrative ex- amples), we report confidence intervals. These include Figure 4. where we compare Gra--Tranerr with other methods, and Supplementary Table $1, which presents an ablation study.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We report compute resources in the main text under the "Training" section.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification:
Our research uses only publicly available datasets and does not involve human subjects, personally identifiable information, or sensitive content. We have reviewed the
NeurlPS Code of Ethics and confirm that our work conforms all guidelines.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
17

---

The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: This work focuses methodological advances in modeling neural spike train data and does not involve applications with direct societal impact.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g.. disinformation, generating fake profiles, surveillance), fairness considerations (e.g.. deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied partieular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our work does not involve releasing models or datasets that pose high risk for misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
18

---

Answer: [Yes]
Justification: use the publicly available Allen Institute Visual Coding dataset, which we properly cite in the paper.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should provided.
assets are released, the license, copyright information, and terms use in package should be provided. For popular datasets, paperswi thcode. com/ datasete has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. this information is not available online, the authors are encouraged to reach out to the asset's 's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
We provide an anonymous GitHub repository containing the full implementa- tion of GLM-Transformer including code for model training, evaluation, and reproduction all experiments. The repository includes documentation and example scripts support ease of use by other researchers.
Justification:
Guidelines:
The answer NA means that the paper does release assets. Researchers should communicate the details of the dataset/code/model part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: This paper does not involve crowdsourcing research with human subjects. All analyses were performed on publicly available neural recordings from non-hume non-human animals.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
19

---

15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification: This study does not involve research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
Declaration LLM usage
Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing. editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
This work does not involve use of large language models (LLMs) as part of its core methods. Any LLM use was limited to standard writing editing assistance and did not affect the scientific content or methodology.
Justification:
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy for what should or should be described.
wwwpeurips.cc/ooofeeeeeeeeeeeeee / /LLM)
20

---

A
Ablation Study
Table S1: Ablation study results. Values test loss in terms of negative oog-likelihod (lower is better), with standard error in parentheses.
Dataset
Model Variant
GLM Simulation
EIF Simulation
Allen Visual Coding
GLM- Transformer (Full)
3.008 10" (2.6 10")
1.513 10" (4.9 10
4.218
10 (3.8 10
Bidirectional RNN encoder
3.032 10 (5.3 x 10")
1,515 x 10" (1.0 x 10")
4.225
10 (2.7 x 10")
No trial-wise latent
3.035 x 107 (4.4 x 104)
1.521 x 107 (7.3 x 10%)
4.264 x 107 (3.3 x 10")
No coupling term
3.065 107 (1.2 10%)
1.537 107 (7.3 10%)
4.263
107 (2.6 x 10")
No post-spike history
3.109 10 (4.1 10%)
1.603
107 (2.8 10f)
4.294
107 (1.0 10
To assess the contribution of key components GLm-Transformrr, conducted series of ablation experiments, with results summarized in Table ST Each variant removes or alters one component of the full model, and we report the resulting test negative log- ligeeeilioood (NLL) on three datasets: synthetic GLM simulation, an EIF simulation, and the Allen Institute Visual Coding dataset. In all experiments, we use the mean of the approximate posterior as the latent without sampling.
Replacing the Transformer encoder with a bidirectional RNN yields slightly higher test loss on the synthetic dataset, with minimal differences the EIF and real datasets. This suggests that while the Transformer may better capture certain trial-specific features, both architectures can perform comparably in practice. Removing the trial-wise latent variable entirely leads to further degradation in performance across all datasets, confirming the importance of capturing trial-to-trial variability in the individual-neuroo dynamics component. Omitting the coupling term or the post-spike history filter significantly worsens model fit.
Overall, the ablation results validate the design of GLM- Transformer: each component contributes meaningfully to the model's performance, and the full model consistently achieves the best results across all settings.
B Supplementary figures
B.1 Non-identifiability example
Supplementary Figure ST illustrates a case of non-identifiability when regularization is removed and the individual-neuron dynamics component made overly ficxible. In this example, the model is configured with larger number of B-spline basis functions (increased from 10 to 50), no smoothness penalty (regularization weight reduced from 1000 to 0), a high-dimensional trial-wise latent variable 2 (dimension increased from 4 to 64), and more latent factors for the individual- dynamics term (increased from 2 to 4). As a result, the individual- neuron dynamics component absorbs temporal structure that should instead be attributed to coupling.
While the total predicted log firing rate closely matches the ground truth, the decomposition into individual-neuron dynamics and coupling is incorrect. This highlights the necessity of adopting modeling assumptions that make the two effects distinguishable. In our case, we impose smoothness on the individual-neure dynamics.
21

---

Individual-neuron dynamics
Coupling
-2.0
0.4
log firing rate
log firing rate
-2.5
0.2
-3.0
Ground truth
0.0
Fitted
-0.2
250
500
250
500
Time
Time
Supplementary Figure S1: Non-identifiability example from the GLM simulation dataset. Left: Ground truth (gray) and fitted (blue) individual-neuron dynamics. Right: Ground truth (gray) and fitted (blue) coupling term. Although the predicted firing rate somewhat correct, the decomposition incorrect: the overly flexible indivonual-neuron dynamics absorbs the coupling component.
B.2 Post-spike filter recovery
All neurons in the GLM simulation share the same ground truth self-history filter, which reflects refractory period. Figure S2 shows that GLM-Traae Trr-rfssfrrm recovers the shape and timing of this filter across neurons.
Estimated and ground truth self-history filter
-0.5-
-1.0
Ground truth Estimated
-2.-
10
Time
Supplementary Figure S2: Recovery the shared post-spike filter the GLM simulation. The estimated filter (blue) matches the ground truth (gray).
B.3 EIF simulation with four neural populations
We conducted an additional EIF simulation involving four neural populations to evaluate the scalability of GL.M-Transformer The ground-truth connectivity was designed such that population projected to populations and 4, and population 2 projected to population 3. G.M--Tranfoommer successfully recovered the true directed interactions, demonstrating ability to generalize to more than two neural populations.
22

---

Synaptic Current (nA)
from to
from 3 1
from
Coupling (log finin rate))
0.1
0.1
0.1
0.2
0.2
0.2
0.0
0.0
0.0
0.0
0.0
0.0
Synapeic Current (nA)
from
from 2
Coupling (log firing rate)
0.2
0.2
0.2
0.0
Synoptic Current (nA)
from 1 to
3
from 2 to
from 4 to 3
Coupling Blog firing rate)
0.1
0.2
0.1
0.0
0.0
0.0
0.0
0.0
0.0
200
400
Synaptic Current (nA)
from
from to
from 3 4
Time (ms)
0.1-
0.1
0.1
Coupling (log firing
0.2
0.2
0.2
0.0
0.0
0.0
0.0
0.0
0.0
200
400
200
400
200
400
Time (ms)
Time (ms)
Time (ms)
Supplementary Figure $3: EIF simulation with four neural populations, each consisting of 50 neurons. Each panel shows the estimated coupling filters (blue) and corresponding ground-truth synaptic currents (green) for all pairs of source and target populations. identifies the true connectivity structure (1-2, 2-+3, 1-)4).
successfully
B.4 Outputs of the six area Allen Institute Visual Coding dataset
We applied GLM-Transformer to six visual cortical areas in the Allen Institute Visual Coding dataset (V1, LM, AL, RL, AM, PM). Each subplot shows the baseline and estimated coupling effects from one source area to a target area. Consistent with known anatomical hierarchies, the model identified strong feedforward connections from VI to higher visual areas, especially V1-1L and -I-AA Therefore, we focus on V1-+LM in Figure 3 These results further suggest that GL GLM--ransformer can extract interpretable and biologically meaningful connectivity patterns from large-scale neural recordings.
23

---

individual-nerron dynamics
Coupling
PM to AM
V1
to AM
to AM
0.51
0.5
0.5
055
0.t
-4.0
0.0
0.0
0,0
0.0
PM
AM to PM
V1 to PM
LM to PM
to PM
RL to PM
-3.5
AM
V1
LM V1
to
0,5
4.0
0.5
0.5
05-
0.5
-4.0
0.0
0.0
0.-
0.0
AL
AM to AL
to AL
V1 to AL
LM to AL
RL to Al.
3355
0.0
0
0.0
RL
AM RL
PM RL
VI to
LM
AL to RL
3.5
-..0
200
200
200 400
400
200
Time (ms]
Time (ms)
Time Ims)
Time (ms)
Time (ms)
Time (ms)
Supplementary Figure S4: Outputs of the six areas in the Allen Institute Visual Coding dataset. Each row shows the individual-neuron dynamics (left column) and estimated coupling effects (remaining columns) for one target area. Strong coupling is observed from VI to LM and from V1 to AL, consistent with the known feedforward hierarchy of the mouse visual cortex.
C Regularization terms
We two forms of regularization to improve identifiability and interpretability.
Smoothing penalty. To encourage smooth temporal structure in the individual-neuron dynamics, apply le penalty the second-order finite difffrrences of each latent factor. Specifically, let f(t) E R²T denote a time-varying factor evaluated at I discrete time bins. The smoothing penalty is:
((t+))-2f()++f)) 1))²
Ysmooth
where Ysmooth is tunable regularization weight. This penalizes curvature in the factor trajectory, promoting smoothness over time.
Sparsity penalty. To encourage sparsity in the coupling structure, apply é penalty the sending and receiving weights of the coupling factors.
D Shared-across-animee and animal-specific components
To encourage the model to learn fundamental features that generalize across animals, we promote parameter sharing wherever feasible. For the individual-neuron dynamics component, the first and last processing steps are kept animal-specific, as they directly operate on neurons, and each animal involves a different set of neurons. However, all intermediate steps-inc steps-including the Transformer encoder and MLP decoder-are shared across animals and brain areas.

---

In contrast, the coupling component involves only two processing steps, and none of its parameters shared across animals. This design choice allows the model capture animal-specific variations in functional interactions across brain areas.
Hyperparameter tuning
For the Allen Institute Visual Coding dataset, we performed hyperparameter tuning for GLM- Transformer using Bayesian optimization with the hyperopt library and the Tree-structured Parzen Estimator (TPE) algorithm. The objective was to minimize the test negative log-likelihood.
Fach trial followed the four-stage training procedure described in the main text, and evaluation used the mean of the approximate posterior (i.c., without sampling). Due to computational constraints, hyperparameter optimization conducted on subset of the data consisting of two animals (IDs: 757216464 and 715093703), and the resulting optimal configuration was used for all real-data experiments.
Search space. The following parameters were tuned:
B-spline basis: number of basis functions in {10, 30, 50)
Transformer encoder: number layers in (1, 2. 4), hidden size dmosel {128, 256, 512),
feedforward dimension in (256, 512, 1024], dropout in (0.0. 0.2, 0.4), number of attention heads in (1, 2), and token time bin length (ms) in (10, 20}
Trial-wise latent variable 2,:;; dimensionality d2 in (2, 4, 8)
Decoder: number of latent factors in 2, and MLP hidden layer size in (dz, 2dz, 4d2
Coupling: number of basis functions in |3, 5), window length (maximum delay allowed) in 130, 50, 70), coupling subspace dimensionality [1, 2]
Penalties: spline smoothing in (100, 1000, 10000), coupling subgroup penalty in (le-5, le-4, 1c-3, 1e-2)
Optimization: Adam optimizer with learning rates for each component in (1e-4, 1e-3, 1e-2), batch size in (32, 64, 128)
Final selected hyperparameters. The best-performing configuration was:
B-spline basis: 10 basis functions
Transformer encoder: layer, dmodel = 128, feedforward dimension 512, dropout = 0.0, attention head, token time bin length 20 ms
Trial-wise latent variable z.:²: d2 =
Decoder: 2 individual- dynamics factors, MLP hidden layer size 4
Coupling: 3 basis functions, window length (maximum delay allowed) 50, subspace dimension
Self-history filters: basis functions, window length (maximum delay allowed) 10
Penalties: spline smoothing = 1000, coupling subgroup penalty = le-5
Optimization: Adam optimizer with transformer encoder learning rate le-4, decoder and coupling learning rates = le-2, and batch size = 64.
For the GLM simulation dataset, we use the same set of hyperparameters as those selected for the smaller spline smoothing penalty
Allen Institute data. For the EIF simulation dataset, (= 100) and a shorter coupling effect window (35 ms); all other hyperparameters remain the same.
We also explored the robustness of the model by jittering few hyperparameters around the final selected values. The results remained consistent: the inputs from V1 LM and from V1 to AL were almost always the strongest coupling effects.
25

---

F Exponential integrate-and-fire (EIF) neuron model
We used the exponential integrate-and-fire (EIF) neuron model
49
generate synthetic 50
51
datasets. Unlike the simpler leaky integrate-and-fire (LIF) model, the EIF model more accurately captures the exponential rise in membrane potential near the firing threshold and is considered more biologically realistic [49]. Our simulated spiking network consists of two neural populations, source and target, each containing 50 neurons. The membrane potential V (t) of neuron n evolves according to:
Cm -g ET) ATe + (1) p (1) (t). di
dVn
spike is emitted when Vn,(t) reaches the threshold Vih. The neuron then enters a refractory period Teer after which its membrane potential is reset to Ve-. We use the same parameters as in Huang al. [5]] Tm Cm//L = 15 ms, EL = -60 mV, Vi -50m mV, Vi 10 AT 2mV,
V2 65 and Tiet = 1 ms.
The input terms T (1), [eoise (I), and (1) represent synaptic input, background noise, and external stimulus, respectively. /noiss (t) is modeled as Gaussian white noise with amplitude selected to yield a baseline firing rate of 10 Hz in the absence of synaptic external inputs. The synaptie input current given by:
(t) Cm
N Σ 3
Jn,m
Σ
n)(
t
(8)
where Jnmm the synaptic weight from neuron m to neuron n, and t denotes the ktth spike time of neuron m. The postsynaptic current kernel 17(t) is defined as:
-l/ra e//rr,
t 0, t < 0,
711) = Td
0,
with тa = 5 ms and ग. = ms. Synaptic weights are sampled from a log-normal distribution 52 In Figure 3 the mean synaptic weight from 10 source neurons to 10 target neurons is 0.2. Recurrent connections within each population have mean weight of 0.01.
The external input (t) consists of two components. The first mimics visual stimuli from the Allen Institute Visual Coding dataset using dual-peaked temporal profile with peaks approximately 50 ms and 250 ms after stimulus onset [46]. The second component introduces trial-to-trial variability through gain modulation. A smooth gain factor is drawn per trial from a Gaussian process (with time constant of 300 ms and amplitude 0.05) and projected to each neuron using a loading vector whose entries are independently sampled from Uniform 00.5, 1.0). In Figure 3 the Gaussian process factors of the two populations are uncorrelated, while they are correlated in Figure 4
Simulations were performed with 0.2 time step, and spike trains were then binned into 2 ms intervals for use in
and other models. A total of 2,000 trials were simulated.
G Animals used in the Allen Institute Visual Coding dataset
All real-data experiments (e.g., Figure 5 and Supplementary Figure S4 were conducted using Neuropixels recordings from ten animals in the Allen Institute Visual Coding dataset. The animal IDs are: 757216464, 79891 1424, 715093703, 719161530, 721 737581020 739448407, 742951821, 743475441, and 744228101.
We used all available trial types. Trials shorter than 400 ms were concatenated with the following trial, since these stimuli (e.g., Gabors, flashes) were presented without intervening blackout periods. For trials longer than 400 ms (e.g., drifting gratings, which last 500 ms), we used only the first 400 ms. For stimulus types that span long continuous durations (e.g., natural movies), we divided the full sequence (c.g., 10 minutes) into non-overlapping 400 ms segments.
26


Reviews:
[{"review_id": "7P5jKak3iS", "reviewer": "Reviewer_pXQU", "summary": "The authors present a hybrid approach to modeling neural population dynamics that accounts for 1) individual neuronal latent dynamics, 2) cross-area coupling, 3) post-spike filters and spike history using a combination of GLMs and transformers. The authors outline the construction of their model and convincingly validate it on both real and simulated (i.e. ground truth) data.", "questions": "I have 3 questions: \n- Why did the authors use a combination of splines, raised cosine basis functions, and “penalized roughness” to estimation the 3 components of their model. Why use 3 different things? Could you not have used B-splines for all 3? It seems to make the model overly-complex this way.\n- You used 20-30 ms time bins for the transformer. Why?\n- Could the authors comment on how data hungry their approach might be compared to end alternative, non-transformer method?", "limitations": "Yes", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The paper was exceptionally easy to read, well organized, and very polished. The approach is novel and demonstrates excellent performance both in terms of ground-truth recovery  and interpretability. I think that the method could be very valuable to analysis of high-density recordings. I have no major criticisms but I do have some questions which I detail below.", "was_revised": false, "final_justification": "", "rebuttal": "We thank the Reviewer pXQU for the helpful comments and for recognizing the value of our work. We address your questions below, and ***we will update the revision to provide more clarification***.\n\n> Q1) Why use a combination of splines, raised cosine basis functions, and “penalized roughness”?\n\nWe agree this could benefit from additional explanation. In principle, several types of basis functions could be used for all components of our model, as they are often roughly equivalent in representational power. However, we chose different basis functions for specific reasons grounded in modeling convention and performance.\n\nWe use raised cosine basis functions for coupling because they are denser at short time lags, which allows for rapid changes immediately following a spike. This aligns well with known properties of post-spike and coupling effects and follows established practices in GLM modeling [1].\n\nFor individual-neuron dynamics, we use uniform B-spline basis functions, which do not have the inductive biases built into the raised cosine basis. \n\nFinally, we apply penalized roughness to encourage smoothness in the estimated individual-neuron dynamics. While it is possible to achieve smoothness using a small number of B-spline basis functions, we believe the penalty approach offers more systematic, automatic, and controllable regularization, as discussed in the spline-based modeling literature [2,3].\n\n> Q2) Why use 20-30 ms time bins for the transformer\n\nThank you for bringing this up! We tried shorter or longer time bin widths; the performance (log likelihood) drops slightly if using time bins bigger than 20ms or smaller than 10ms. The model is relatively insensitive within the 10–20 ms range, and we found no substantive differences between these settings. A possible explanation is that larger time bins lose temporal information, while very small bins require more data to train and may not carry too much additional information.\n\nBelow is a table showing the relative negative log-likelihood (NLL) of the model across different time bin widths for each dataset. Values are offset by the best (minimum) NLL for each dataset, so lower is better, and 0 indicates the best performance.\n\nAlthough the 20 ms bin width appears to perform better than 10 ms in this table, the results here were obtained on a smaller dataset ( 2000, 2000, and 5915 trials for the three datasets, respectively) due to limited time during the rebuttal phase. In contrast, our formal hyperparameter tuning (performed on a larger real dataset)  selected 10 ms as the optimal bin width, which is what we used in the final model.\n\n| Time bin width (ms)  | NLL on GLM synthetic dataset  | NLL on EIF synthetic dataset | NLL on Allen Institute dataset\n| ------------- |:-------------:| ------------- |:-------------:|\n| 2      | 2.43e+04     | 3.39e+04      | 6.39e+04     |\n| 5      | 1.09e+04     | 3.07e+04      | 6.30e+04     |\n| 10      | 1.43e+02     | 9.56e+01      | 1.42e+04     |\n| 20      | 0     | 0      | 0     |\n| 40      | 2.21e+03     | 6.03e+04      | 1.35e+04     |\n\n> Q3) how data hungry compared to non-transformer method\n\nThis is a great point and we will add a comment to the results in the revised manuscript. We evaluated model performance as a function of training set size using our EIF synthetic dataset, comparing the GLM-Transformer with a simpler non-Transformer GLM variant. Both models eventually reach a performance plateau, but the simpler GLM converges quickly with as few as 100 training trials. In contrast, the GLM-Transformer generally requires around 1500 trials to reach a higher-performing plateau.\n\nThis higher data demand reflects the flexibility of the Transformer-based component and its capacity to model complex trial-to-trial variability. One possible future direction when encountering smaller dataset is to pretrain a foundation model across multiple neural datasets, then fine-tuning on the smaller one.\n\n### References\n\n[1] Pillow J W, Shlens J, Paninski L, et al. Spatio-temporal correlations and visual signalling in a complete neuronal population[J]. Nature, 2008, 454(7207): 995-999.\n\n[2] Hastie T, Tibshirani R, Friedman J H, et al. The elements of statistical learning: data mining, inference, and prediction[M]. New York: springer, 2009.\n\n[3] Hastie, T., Tibshirani, R., & Friedman, J. (2009). *The Elements of Statistical Learning: Data Mining, Inference, and Prediction* (2nd ed.). Springer."}, {"review_id": "7HOFANg6V7", "reviewer": "Reviewer_BKwJ", "summary": "This paper introduces GLM-Transformer, a new method for identifying how different brain areas communicate with each other using neural spike data. The main challenge is that neural activity varies significantly from trial to trial due to factors like behavioral changes, which can mask true interactions between brain regions. The authors' approach combines a Transformer-based neural network to capture these trial-to-trial variations in individual neuron activity with a generalized linear model framework to identify directed connections between brain areas. They tested their method on both synthetic data and real recordings from mouse visual cortex, showing that it can accurately recover known connections while remaining robust to confounding background activity that trips up other methods. When applied to the Allen Institute Visual Coding dataset, GLM-Transformer successfully identified feedforward pathways from V1 to higher visual areas, consistent with established brain anatomy. The method offers a way to better understand brain connectivity by separating true inter-area communication from within-area variability.", "questions": "1. Could you explain the mechanistic role of trial-varying individual-neuron dynamics in achieving robustness to shared background fluctuations?\n\n2. Could you design a statistical framework to assess the significance of detected coupling effects and provide confidence intervals for coupling strength estimates?\n\n3. Could you design a simulation experiment to evaluate the model's performance and computational scalability across varying numbers of brain areas (>2)?", "limitations": "The authors discussed the limitations in Section 5, and I agree with their statements.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "**Strengths**:\n\n1. This paper addresses a challenge in computational neuroscience by combining the interpretability of GLMs with the flexibility of Transformer-based VAEs. This hybrid design is creative, allowing the model to capture complex trial-to-trial dynamics while maintaining interpretable cross-area coupling terms.\n\n**Weaknesses**:\n\n1. All experiments are conducted on relatively small-scale systems (2-3 brain areas), with Figure 5 demonstrating interactions only between V1 and LM. This limited scope raises serious questions about the method's scalability to realistic neuroscience applications, which typically involve dozens of brain areas and hundreds of neurons.", "was_revised": true, "final_justification": "My concerns have been resolved based on the new results and the authors' clarifications.", "rebuttal": "We thank Reviewer BKwJ for recognizing that “this hybrid design is creative” and for the thoughtful feedback. We agree it will help strengthen the work. We address your questions below, and ***we will revise the manuscript to reflect these valuable suggestions***.\n\n> Weakness 1) All experiments are conducted on relatively small-scale systems (2-3 brain areas)\n\nWe apologize for not making it clearer in the original submission. In our real dataset experiments, the model was applied to **six brain areas**: V1, LM, AL, RL, AM, and PM, each with many tens of neurons.\n\nThe table below shows the **mean coupling contributions** (in log firing rate, averaged across trials, neurons, and time) between areas. The V1 -> LM coupling is the largest among well-sampled regions, consistent with the known feedforward hierarchy in mouse visual cortex. This is why we focused on V1 -> LM coupling in Figure 5. We will also add a supplementary figure to show this. \n\n|From ↓ To → | V1 | LM | AL | RL | AM | PM\n| ------------- |:-------------:|:-------------:|:-------------:| ------------- |:-------------:|:-------------:|\n| V1      | \\     | 0.153      | 0.075     | 0.029      | -0.006      |-0.023     |\n| LM      | 0.006     | \\      | 0.000     | 0.000      | 0.108      |0.000     |\n| AL      | 0.000     | 0.000      | \\     | 0.000      | 0.036      |0.187     |\n| RL      | 0.024     | 0.000      | 0.000     | \\      | 0.000      |0.000     |\n| AM      | 0.013     | 0.108     | 0.000     | 0.000      | \\      |0.000     |\n| PM      | 0.026     | 0.000      | 0.000     | 0.000      | 0.000      |\\     |\n\nWe also added a new simulated dataset with **four populations** to further support model scalability. Please see our response to Q3 below.\n\n\n\n> Q1) mechanistic role of trial-varying individual-neuron dynamics in achieving robustness to shared background fluctuations\n\nSorry we may not have explained this very well. In practice, each neuron in each trial exhibits background fluctuations, and they are often correlated across areas [4,5,6]. However, Without modeling these, a method might mistake shared background variability for inter-area coupling.\n\nOur trial-specific individual-neuron dynamics component explicitly models this variability, thereby isolating shared background dynamics from true cross-area interactions. This decomposition prevents the model from erroneously attributing correlated activity to coupling. We will clarify this point in the revised text.\n\n> Q2) design a statistical framework to assess the significance of detected coupling effects and provide confidence intervals for coupling strength estimates\n\nWe agree that developing such inference procedures is a very important next step. While it is beyond the scope of the current work, we are actively exploring approaches inspired by recent work on post-hoc inference in deep learning, such as double/debiasing machine learning [7]. These ideas could be extended to estimate uncertainty around coupling effects in GLM-Transformer, and we view this as a promising direction for future work.\n\n> Q3) design a simulation experiment to evaluate the model's performance and computational scalability across varying numbers of brain areas (>2)\n\nWe appreciate this suggestion and have conducted a new simulation experiment with **four neural populations** to directly address this point.\n\nWe created a new EIF dataset with populations A, B, C, and D, each containing 50 neurons, to mimic the activity patterns of V1, LM, AL, and RL in the real dataset. The ground truth connectivity structure was designed as follows: A to B, B to C, and A to D.\n\nThe table below shows the mean contributions between areas (in log firing rate, averaged across trials, neurons, and time bins). As shown, GLM Transformer successfully recovers the ground truth connections:\n\n|From ↓ To → | A | B | C | D |\n| ------------- |:-------------:|:-------------:|:-------------:| ------------- |\n| A     | \\     | 0.219      | 0.000     | 0.108      |\n| B     | 0.000     | \\      | 0.473     | 0.000      |\n| C     | 0.000     | 0.000      | \\     | 0.000      |\n| D     | 0.000     | 0.000      | 0.000     | \\      |\n\nThis larger scale dataset demonstrates that GLM-Transformer can scale to more complex multi area networks, and we plan to replace the original two population simulation (Figure 2) with this updated result in the revised manuscript.\n\n### References\n\n[4] Musall S, Kaufman M T, Juavinett A L, et al. Single-trial neural dynamics are dominated by richly varied movements[J]. Nature neuroscience, 2019, 22(10): 1677-1686.\n\n[5] Stringer C, Pachitariu M, Steinmetz N, et al. Spontaneous behaviors drive multidimensional, brainwide activity[J]. Science, 2019, 364(6437): eaav7893.\n\n[6] International Brain Laboratory, Benson B, Benson J, et al. A brain-wide map of neural activity during complex behaviour[J]. Biorxiv, 2023: 2023.07. 04.547681.\n\n[7] Kennedy E H. Semiparametric doubly robust targeted double machine learning: a review[J]. Handbook of statistical methods for precision medicine, 2024: 207-236."}, {"review_id": "y8ZOARaMHx", "reviewer": "Reviewer_pEdw", "summary": "This paper introduces GLM-Transformer, a hybrid model that combines a generalized linear model (GLM) for interpretable cross-population coupling with a Transformer-based variational autoencoder to capture individual-neuron dynamics and trial-to-trial variability. The authors estimate a low-dimensional latent vector per trial and population,m which maps onto B-spline-based background firing rate trajectories. This is combined with self-history and coupling terms that are modeled via low-rank temporal filters. The model is trained in four stages, iteratively increasing complexity, with an ablation study showing that every component meaningfully improves performance. The authors validate their approch on synthetic, mechanistic, and real spike train data, where it recovers known interaction structure.", "questions": "- assumption of fixed trial length: I initially found it a bit hard to understand how time-dependence arises from the fixed trial latent factors via the splines. It seems the approach assumes the splines are defined over a fixed temporal support. Could you clarify how restrictive this assumption is, e.g., how would the model handle variable-length trials or generalization beyond the training interval? You mention the fixed length briefly in the Discussion (“we also were able to assume a fixed trial length...”), but this assumption appears central and could benefit from more explicit discussion directly  in the Methods section.\n\n- You mention that “in all experiments, we use the mean of the approximate posterior as the latent without sampling.” Out of curiosity: how large is the latent variance in practice? What would happen if one trained the model deterministically without a VAE setup? This would help clarify whether the probabilistic treatment adds substantial value (could also be added to table S1 if time allows)\n\n- Table S1:The absolute differences in test loss are very small (e.g., 0.0567 vs. 0.0568 for Transformer vs. RNN encoder). Are these differences practically meaningful? Do they persist across longer trials or different datasets?\n\n- Encoder: What adjustments, if any, were necessary to adapt the Transformer to spike data? Did you mainly take the settings from the literature? Are there specific preprocessing or normalization steps for the token sequences? While the hyperparameter search for the Transformer is in the supplement, no other details are given as far as I could see.\nSecond, how does computational load compare to the RNN variant, particularly as trial length increases? What is the computational load vs. the RNN/runtime? Especially with respect to trial length?", "limitations": "Clear limitations section.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 2, "strengths_and_weaknesses": "*Strengths*:\n\nThis is a well-written and logically structured paper that presents a principled approach to disentangling background activity and cross-population coupling in neural spike train data. The introduction is strong, and the method is clearly motivated. The combination of a GLM framework with a (Transformer-based) variational autoencoder is sensible and timely, and the authors have done a good job of situating their work within the relevant literature on spike train modeling. The four-stage training procedure and the corresponding ablation study in Table S1 help clarify the contributions of individual components. The results section is well-organized.\n\n*Weaknesses*:\n\nModeling assumptions: While much of this is unavoidable when combining multiple components into sophisticated architectures, there are some points in the methods section where modeling assumptions or design decisions are introduced without enough justification. In particular:\n- \"to avoid this, we make the seemingly reasonable assumption that coupling \neffects vary more rapidly than the background fluctuations arising from changes in behavior or endogenous state\". It’s fair that we need to introduce some inductive biases to assure identifiability, but this assumption is not very intuitive to me. In practice, behavioral or arousal-related state changes (e.g., due to attention or movement) can also occur on fast time scales? What would be a scenario in which this breaks down?\n\n- “We use 10 or 20 ms time bins as tokens for the Transformer (which are wider than the 2-ms bins used elsewhere in the model)” Why? I understand that this may be driven by attention scaling (e.g., quadratic cost), but an explicit justification would be useful.\n-“Throughout the paper, we use a bin size of 2 ms and rarely observe\nmore than one spike per bin” Have you considered alternative count models (e.g., zero-inflated Poisson or negative binomial)?\n- The posterior over trial-wise latent variables assumes a diagonal covariance matrix. While this is common, some discussion of whether this is a modeling assumption or a practical constraint (e.g., for stability or training speed) could be helpful.\n\nClarity of contribution:\nWhile the method is well-motivated, I found the paper could more clearly articulate what is actually novel in terms of methodology. Many of the model components, e.g., trial-wise latent variables, B-spline decoders (as in Deep Random Splines [36]), smooth temporal filters via raised cosine bases [6], and Transformer encoders adapted to spike data [31], are taken from prior work. That's entirely valid, and combining these elements in a novel  architecture can itself be a contribution. However, a short summary paragraph at the end of the Methods section highlighting which components are standard, which are adapted, and which are new would be really helpful.\n\nFraming and overall contribution:\nRelating to the clarity of contribution, one point that could also benefit from clarification is the overall goal of the paper. The title emphasizes “identifying interactions across brain areas”, which would suggest that the coupling terms are the object of interest. However, much of the model’s complexity and interpretability comes from the other components (the individualized response functions, trial-specific latent variables and self-history terms). It’s not always clear how central they are to the paper’s scientific aims. Are the latents z meant to be interpreted (e.g., in relation to behavior, as in Figure 5E), or are they primarily there to recover accurate coupling terms? A clearer statement early in the paper, and perhaps reiterated at the end of the Methods section, about what is meant to be interpreted versus abstracted away would be helpful.\n\nSuggestions:\nThe related work section is concise and well-organized but somewhat narrow in scope. It focuses primarily on two points: (1) that Deep Random Splines (DRS) uniquely capture trial-level variability, and (2) that prior models do not provide interpretable population-level coupling. However, trial-specific modeling is an increasingly active area, with several recent approaches proposing hierarchical or meta-dynamical structures that extract trial-level latents in a dynamics context, e.g. \nMeta-Dynamical State Space Models for Integrative Neural Data Analysis (Vermani et al, ICLR 2025)\nLearning Interpretable Hierarchical Dynamical Systems Models from Time Series Data (Brenner et al., ICLR 2025) \nGeneralizing to New Physical Systems via Context-Informed Dynamics Model (Kirchmayer et al., ICML 2022).\n\nAlso, the claim that prior latent space models \"do not provide direct interpretation\" of cross-population effects could be clarified: why can such interpretability not also be achieved through structured latent recurrences or constrained connectivity (e.g., low-rank or sparse matrices). Finally, some recent work also focuses on directly inferring latent factors from spike trains that may also be worth referencing: e.g. \nLatent Diffusion for Neural Spiking Data (Kapoor et al., NeurIPS 2024)\nIntegrating Multimodal Data for Joint Generative Modeling of Complex Dynamics (Brenner et al., ICML 2024)", "was_revised": true, "final_justification": "As highlighted in my final comment, a few limitations remain. First, the Transformer encoder, despite its prominence in the title, yields virtually no performance gain over a bidirectional RNN in your ablation (Table S1), while the suggestion that it will outperform on other datasets feels a bit speculative. Second, the fixed-length-trial assumption is quite substantial, since in many real-world settings trials vary substantially. Finally, as mentioned by other referees, the methodological novelty beyond the low-rank GLM integration feels incremental. I usually don't like raising the novelty claim, since it is hard to argue against, and the authors rightly acknowledge that they're building on considerable prior work, but for that reason I’d expect the other parts of the paper (e.g., more compelling empirical gains or broader applicability) to compensate.\n\nThat being said, overall, the GLM-Transformer represents an interesting integration and yields solid empirical results, but these remaining concerns keep me on the fence. I will therefore retain my current borderline-accept rating.", "rebuttal": "We thank Reviewer pEdw for the detailed and constructive review and for describing our work as “sensible and timely.” We greatly appreciate that you read deeply with the main text and supplementary materials. Your feedback has already helped us identify several ways to improve the work. We address your comments below, and ***we will revise the manuscript accordingly***.\n\n> Weakness 1) Modeling assumptions: assumption to assure identifiability\n\nThis is good criticithsm. We overstated that point, and we changed \"we made a seamingly reasonable assumption\" to \"we made a strong assumption\" in our revision. \n\nOur reasoning is based on the fact that slower fluctuations in neural activity tend to have higher power in the spectrum [8]. By regularizing the individual-neuron dynamics to be smoother, we are effectively adjusting for background fluctuations that dominate lower-frequency bands. \n\nAlso, in the Allen Institute Visual Coding dataset, we believe the major sources of trial-to-trial background fluctuation are behavioral factors such as locomotion, arousal, and reaction time. These tend to evolve on the scale of hundreds of milliseconds, whereas cross-area interactions typically occur at much faster time scales [9]. So, we believe our assumption holds in this setting and in many similar experimental paradigms. That said, we acknowledge that in other situations (e.g., fast-changing internal states), the assumption may be violated.\n\nLastly, we want to emphasize that we view functional coupling as capturing conditional dependencies, not causality. Even if some fast-changing unobserved confounders exist, the estimated coupling that is conditioned on the slow-changing background component, can still reflect useful and interpretable patterns of interaction across areas.\n\n> Weakness 1) Modeling assumptions: explicit justification for 10ms or 20ms time bin\n\nAnother reviewer also brought this up. We tried shorter or longer time bin widths; the performance (log likelihood) drops slightly if using time bins bigger than 20ms or smaller than 10ms. The model is relatively insensitive within the 10–20 ms range, and we found no substantive differences between these settings. A possible explanation is that larger time bins lose temporal information, while very small bins require more data to train and may not carry too much additional information.\n\nBelow is a table showing the relative negative log-likelihood (NLL) of the model across different time bin widths for each dataset. Values are offset by the best (minimum) NLL for each dataset, so lower is better, and 0 indicates the best performance.\n\nAlthough the 20 ms bin width appears to perform better than 10 ms in this table, the results here were obtained on a smaller dataset ( 2000, 2000, and 5915 trials for the three datasets, respectively) due to limited time during the rebuttal phase. In contrast, our formal hyperparameter tuning (performed on a larger real dataset)  selected 10 ms as the optimal bin width, which is what we used in the final model.\n\n| Time bin width (ms)  | NLL on GLM synthetic dataset  | NLL on EIF synthetic dataset | NLL on Allen Institute dataset\n| ------------- |:-------------:| ------------- |:-------------:|\n| 2      | 2.43e+04     | 3.39e+04      | 6.39e+04     |\n| 5      | 1.09e+04     | 3.07e+04      | 6.30e+04     |\n| 10      | 1.43e+02     | 9.56e+01      | 1.42e+04     |\n| 20      | 0     | 0      | 0     |\n| 40      | 2.21e+03     | 6.03e+04      | 1.35e+04     |\n\n\n> Weakness 1) Modeling assumptions: posterior assumes a diagonal covariance matrix\n\nThis is primarily a practical constraint. Modeling a full covariance matrix would require predicting and regularizing a positive semi-definite matrix with quadratic complexity, which would substantially increase model and training complexity. While possible, we found diagonal posteriors sufficiently expressive in practice, and the added cost was not justified for our current goals.\n\n> Weakness 2) Clarity of contribution\n\nWe will absolutely add a summary paragraph at the end of the Methods section to clearly delineate which model components are standard, adapted, or novel. Due to length limitations in this rebuttal, we outline the key points below:\n\n- The individual-neuron dynamics component adapts ideas from Deep Random Splines (DRS) [10] and Neural Data Transformer (NDT) [11]: we use DRS’s latent-per-trial structure and B-spline decoder, NDT’s Transformer architecture, and dimension reduction to couple multiple neurons through shared latent structure.\n- The low-rank GLM-style coupling component is, to our knowledge, a novel integration. While inspired by reduced-rank regression, we are not aware of prior work combining it with spike-based GLMs to model directed interactions.\n- The post-spike history term is standard in GLM literature.\n- The overall design, which explicitly separates trial-varying background dynamics from coupling effects by leveraging the complementary strengths of GLMs and deep learning models, is also novel and central to our contribution.\n\n> Weakness 3) Framing and overall contribution\n\nThank you for pointing this out. Our primary scientific goal is to estimate cross-area coupling. The individual-neuron latent component is treated as a nuisance term, which is important for isolating coupling, but not the focus of interpretation. In Figure 5E, we show that the latent aligns with behavior to help validate that it is capturing meaningful trial-to-trial structure. We will revise both the Introduction and Methods to better clarify this point.\n\n> Weakness 4) related work section\n\nThank you for the helpful suggestions! We will incorporate the referenced recent works on trial-wise latent modeling.\n\n> Weakness 5) prior latent space models\n\nYes, we agree that interpretability can also be achieved through structured latent connectivity, and there are many valuable models that take this approach. We think these are all important models for understanding population dynamics.\n\nThat said, our GLM-style coupling term captures interactions in a different way. Rather than modeling communication from latent variables to latent variables, it directly estimates the effect of a spike in one neuron on the firing rate of another neuron at different time lags. This more explicit spike-to-rate formulation offers interpretable coupling filters that are especially intepretable for studying directed interactions between brain areas. We see this as a complementary perspective to latent-to-latent coupling approaches.\n\n> Q1) assumption of fixed trial length\n\nWe agree this is a central assumption and currently a limitation of our model. For datasets with a few discrete trial lengths, one workaround would be to learn separate decoders per trial length. A more general and exciting direction would be to use a GPT-style decoder that supports variable-length outputs.\n\n\n> Q2) latent variance\n\nThank you for catching this. In our ablation experiments, we disabled sampling to reduce variance in training and isolate the impact of other components. However, when sampling is turned off entirely, the posterior variance is no longer trained. In the real-data experiments, we do include sampling. The average standard deviation of the latent posterior is 0.57, averaged across dimensions and trials, indicating meaningful variability is captured.\n\n> Q3) Table S1\n\nWe apologize for the confusion and thank you for pointing this out. The values in Table S1 are the negative log likelihood averaged over hundreds of thousands of time bins, so they are much smaller than the actual total negative log likelihood, which would require multiplying these numbers by the total number of time bins. This averaging compresses the scale and can make differences look deceptively small. We will change all the metrics to total negative log likelihood in the revision.\n\nAs for the Transformer vs. RNN encoder comparison, the difference in negative log likelihood on the Allen dataset is within the confidence interval, so we don’t claim a strong advantage in this setting. However, we believe the benefit of the Transformer may become clearer in settings with longer trials or more complex neural dynamics, where its ability to model long-range dependencies might provide more value.\n\n> Q4) Encoder details\n\nWe did not apply any special preprocessing or normalization to the token sequences; we largely followed prior practices from Neural Data Transformer [11]. In terms of runtime, we did not observe a major difference between the Transformer and the RNN encoder in our experiments. This is likely because the number of time bins (i.e., tokens) per trial is relatively small. In practice, the main bottleneck was converting sparse spike times to dense spike trains during data loading, which dominated the overall compute time and reduced the impact of encoder choice.\n\n\n### References\n\n[8] Kass R E, Bong H, Olarinre M, et al. Identification of interacting neural populations: methods and statistical considerations[J]. Journal of Neurophysiology, 2023, 130(3): 475-496.\n\n[9] Chen Y, Douglas H, Medina B J, et al. Population burst propagation across interacting areas of the brain[J]. Journal of Neurophysiology, 2022, 128(6): 1578-1592.\n\n[10] Loaiza-Ganem G, Perkins S, Schroeder K, et al. Deep random splines for point process intensity estimation of neural population data[J]. Advances in Neural Information Processing Systems, 2019, 32.\n\n[11] Ye J, Pandarinath C. Representation learning for neural population activity with neural data transformers[J]. arXiv preprint arXiv:2108.01210, 2021."}, {"review_id": "PVpCVXyqDr", "reviewer": "Reviewer_wbH6", "summary": "This paper proposes a GLM-Transformer for multi-population and multi-trial neural spike data. The GLM-Transformer generalizes the widely used GLM with a Transformer-based VAE, which decomposes neural activity into three components: individual-neuron dynamics, self-history effects and cross-population coupling. This decomposition makes the model flexible but keeps the interpretability as in GLM. After validating by simulations, they further applied it to the Allend Institute Visual Coding dataset and identified feed forward pathways consistent with established visual hierarchies.", "questions": "The paper is clearly written. Besides the question in the weakness section, I don’t have other major questions.", "limitations": "The paper has discussed the limitations in the Discussion section: dependence of the Transformer on training, the assumption and inference.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "**Strengths**:\n- The work is well-motivated by challenges in analyzing real neural spike train data.\n- The writing is clear, structured, and easy to follow.\n- The model extends the GLM using a Transformer-based VAE in latent space, enhancing flexibility while preserving interpretability through a biologically meaningful decomposition.\n\n**Weaknesses**:\n- The GLM-Transformer combines GLM with Transformer-based VAE in a straightforward way, which may raise questions regarding novelty.\n- Although the method focuses on the interpretability and inference. It would be better to compare it with other methods. Will simple GLM based method be enough to gain insights/ biological interpretation?", "was_revised": true, "final_justification": "I am between borderline accept & reject, as reasons provided in the response. Will let AC make final decision.", "rebuttal": "We thank reviewer wbH6 for their helpful review and for recognizing our work is \"well-motivated by challenges in analyzing real neural spike train data\". We address your questions below. \n\n> Weakness) The GLM-Transformer combines GLM with Transformer-based VAE in a straightforward way, which may raise questions regarding novelty\n\nWe appreciate this comment and would like to clarify our perspective. We believe novel methods often build on existing techniques. The GLM-Transformer (a) decomposes coupling identification from background dynamics and (b) takes advantage of the Transformer to do what it is good at: flexibly fitting complicated temporal processes where interpretation is not needed. This is a novel conception and we spent considerable time and effort to demonstrate its viability.\n\n> Weakness) Although the method focuses on the interpretability and inference. It would be better to compare it with other methods. Will simple GLM based method be enough to gain insights/ biological interpretation\n\nWe fully agree about the importance of comparing to simpler GLM-based baselines and we did include such a comparison. As shown in Figure 3, we directly compared GLM-Transformer against simple GLMs and other alternatives, in the presence of shared background dynamics. These comparisons demonstrate that while standard GLMs can provide interpretable components, they lead to biased coupling when background dynamics are present. GLM-Transformer remains robust in this setting by modeling trial-varying individual-neuron dynamics explicitly, and we discuss this result in the main text."}]

Author rebuttal:

Meta-review: This paper proposes GLM-Transformer, a hybrid framework that combines a Transformer-based VAE to capture trial-to-trial, neuron-specific variability with a GLM-style coupling component to identify interpretable cross-area interactions. The approach directly addresses a key confound in multi-area neural data analysis: background fluctuations that can obscure true directed connectivity. Through simulations, mechanistic models, and application to the Allen Institute Visual Coding dataset, the method is shown to recover known hierarchical pathways while remaining robust to shared variability.

Reviewers consistently found the paper well-written, technically sound, and relevant, highlighting the creative integration of GLMs and deep latent models. Strengths noted include the principled separation of nuisance dynamics from coupling, the biological interpretability of the filters, and the robustness demonstrated across datasets. The authors’ rebuttal clarified several concerns, adding analyses on larger populations, new simulations with four areas, and expanded discussion of assumptions and related work. These additions resolved most critical questions, and two reviewers raised their scores after rebuttal.

The main reservations concern the degree of methodological novelty, the limited advantage of the Transformer encoder over simpler RNNs in ablations, and the fixed-trial-length assumption, which could restrict applicability. Still, the consensus is that the paper makes a meaningful contribution by framing and solving an important challenge with a practical, interpretable, and scalable approach. On balance, I recommend acceptance.
Decision: Accept (poster)