Self-Supervised Discovery of Neural Circuits in Spatially Patterned Neural Responses with Graph Neural Networks
Authors: Kijung Yoon
Venue: NeurIPS 2025 poster ()
Keywords: Neuroscience
Memory
Neural connectivity inference
Continuous attractor networks
Neural circuit dynamics
Graph neural networks
Primary area: neuroscience_and_cognitive_science
TL;DR: This paper aims to uncover true synaptic connectivity from neural population activity using a self-supervised learning approach that addresses confounding effects from correlation-driven artifacts.

Abstract: Inferring synaptic connectivity from neural population activity is a fundamental challenge in computational neuroscience, complicated by partial observability and mismatches between inference models and true circuit dynamics. In this study, we propose a graph-based neural inference model that simultaneously predicts neural activity and infers latent connectivity by modeling neurons as interacting nodes in a graph. The architecture features two distinct modules: one for learning structural connectivity and another for predicting future spiking activity via a graph neural network (GNN). Our model accommodates unobserved neurons through auxiliary nodes, allowing for inference in partially observed circuits. We evaluate this approach using synthetic data generated from ring attractor network models and real spike recordings from head direction cells in mice. Across a wide range of conditions, including varying recurrent connectivity, external inputs, and incomplete observations, our model reliably resolves spurious correlations and recovers accurate weight profiles. When applied to real data, the inferred connectivity aligns with theoretical predictions of continuous attractor models. These results highlight the potential of GNN-based models to infer latent neural circuitry through self-supervised structure learning, while leveraging the spike prediction task to flexibly link connectivity and dynamics across both simulated and biological neural systems.

Self-Supervised Discovery of Neural Circuits in Spatially Patterned Neural Responses with Graph Neural Networks
Kijung Yoon
Department of Electronic Engineering
Department of Artificial Intelligence
Hanyang University
Scoul, Korea 04763
kiyoon@hanyang ac
Abstract
Inferring synaptic connectivity from neural population activity is a fundamental challenge in computational neuroscience, complicated by partial observability and mismatches between inference models and true circuit dynamics. we propose a graph-based neural inference model that simultaneously predicts neural activity and infers latent connectivity by modeling neurons as interacting nodes in graph. The architecture features two distinct modules: for learning structural connectivity and another for predicting future spiking activity via a graph neural network (GNN). Our model accommodates unobserved neurons through auxiliary nodes, allowing for inference in partially observed circuits. We evaluate this approach using synthetic data generated from ring attractor network models and real spike recordings from head direction cells in mice. Across a wide range of conditions, including varying recurrent connectivity, external inputs, and incomplete observations, our model reliably resolves spurious correlations and recovers accurate weight profiles. When applied real data, the inferred connectivity aligns with theoretical predictions of continuous attractor models. These results highlight the potential of GNN- models infer latent neural circuitry through self-supervised structure learning, while leveraging the spike prediction task flexibly link connectivity and dynamics across both simulated and biological neural systems.
this study,
Introduction
Understanding how neural circuits compute and adapt requires identifying the strength of synaptic transmission between neurons, as this knowledge reveals how the structure of a circuit shapes its computational properties. Advances in recording simultaneous activity from large populations of neurons have driven significant interest in using statistical methods [1-7] to infer interactions and estimate connectivity between neurons across entire circuits. Despite this progress, statistical inference methods face two primary challenges: first, no recording technique can capture the activity of all neurons within cireuit, and second, the inference models may fail to accurately represent the underlying generative dynamical system. As a result, these limitations can lead to substantial differences between inferred and true connectivity.
The discrepancy in connectivity inference primarily originates from instances in which weakly connected unconnected neurons exhibit strong activity correlations. This is characteristic feature of strongly recurrent networks, which maintain persistent memory states through the principle of pattern formation [8-10 process where simple, spatially localized competitive interactions produce stable spatial activity patterns. To address these challenges, we design continuous attractor networks capable of generating spatially patterned neural responses IT [5 and focus on evaluating
39th Conference Neural Information Processing Systems (NeurIPS 2025).

---

how effectively a circuit inference model can explain away the correlations that arise from co-activated neurons.
To this end, we propose the use of graph neural networks (GNNs), a robust framework for modeling complex dynamics physical systems with numerous interacting entities, such as particles or atoms [16-19] In our approach, neurons are represented as nodes and their connections edges within GNN-based architecture. The model incorporates two functionally distinct modules: one designed to learn structural connectivity across the network and another dedicated to predicting spiking activity based on simultaneous, circuit-wide neural recordings. Instead of relying on supervised learning, we aim to extract network connectivity in self-supervised manner by training the model to predict subsequent spike events in neural populations, with the model's latent representation serving to deseribe the inferred connectivity. An additional feature of our model is its ability to account for unobserved neurons by adding extra nodes to the graph. Under a transductive framework [20] 21], the connectivity among observed neurons is inferred, allowing message passing occur throughout the entire neural circuit, including both observed and hidden neurons. This formulation implicitly leverages hidden neurons to explore the influence of unobserved components in circuit inference.
We perform extensive experiments on neural spike data generated from highly structured ring networks under various conditions systematically evaluate the inference performance of the proposed framework. Our analysis begins with a fully observed network without external input drives, enabling us to isolate challenges arising from mismatches between the generative system and the inference model. In this setting, we demonstrate that the proposed approach resolves correlations from unconnected neurons at least 70% more effectively than advanced statistical inference methods. Furthermore, we show that the improved quality of circuit inference consistently holds across diverse configurations of the generative model, including stimulus-driven conditions, different recurrent weight profiles, fully versus partially observed networks, and extends to real neural recordings from behaving mice.
The paper is structured as follows. Section 2 provides an overview of related work. Section introduces the recurrent network models used to generate neural spike data and details the proposed GNN-based inference framework. Section presents the experimental results, while Section concludes the study.
2 Related Work
Estimating network connectivity from large population recordings has been a long-standing challenge in computational neuroscience. One prominent line of research focuses on probabilistic modeling techniques, including maximum entropy-based inverse Ising models [2 2 23 and minimum probability flow (MPF) 124, 25). Both approaches leverage the Ising model to capture pairwise interactions among binary variables, such neuronal spiking activity, to reconstruct functional connectivity graphs. Maximum entropy models ensure interpretability by maximizing likelihood under empirical constraints (e.g., correlations) but require computationally expensive estimation of the partition function. In contrast, MPF addresses this limitation by minimizing the probability flow between observed and unobserved states, bypassing the partition function. This makes MPF more scalable and computationally efficient for inferring connectivity in large neural networks.
Another widely used method involves l lgrregularized logistic regression to promote sparsity in connectivity estimates by penalizing the number of nonzero parameters [26] 27]. In this framework, & firregularized logistic regression is performed for each variable against all others, with the sparsity pattern of the regression coefficients used to infer the network's neighborhood structure. This technique is particularly effective for high-dimensional Ising model selection, supporting connectivity inference for large scale datasets with complex correlations. A related framework employs the use of generalized linear models (GLMs), originally formalized by Nelder and Wedderburn (28). to relate a linear predictor to an output variable through a link function. GLMs have been extensively applied to model spatio-temporal interactions and stimulus dependencies 29 30
predicting circuit activity by associating observed spiking activity with intrinsic factors such as spike history and external covariates like stimuli or movement. By explicitly modeling the influence of neurons on one another, these saelihooddbasee approaches treat network connectivity as parameters to be learned, reducing spurious interactions in network inference.
Model- based approaches have demonstrated significant success in characterizing neural interactions and dependencies, especially within sensory systems. However, even with extensive neuronal activity recordings, these models have been reported to inaccurately estimate effective connectivity in memory-

---

a
Conrestivtty W
Populatin & Splie Train
Strucure Leaming Module
Spike Prediction Madule
Call infox
Activation
Inbeex
100000000000
weight
Coll Index
Tima
OOO GOODOO O
Cal Indox
Weightt Matnis
(min)
(nn)
Figure 1: Overview of the generative and inference network model. (a) A ring network with Mexican-hat connectivity profile, where each neuron is connected to others with the same weight pattern. The full weight matrix is shown below. (b) Simulated network activity, including synaptic activation (top) and spike raster plot across neurons over time (bottom). (c) Structure learning module that estimates synaptic connection strengths based on pairwise spike activity. (d) Spike prediction module that leverages the inferred connectivity to predict future spike times from past neural activity.
related (i.e., strongly recurrent) networks compared to sensory-driven circuits 31 the extent to which inferred connectivity faithfully reflects biological neuronal connections remains an open question-one that we seek to explore in this study.
Consequently,
3 Method
In this section, we introduce strongly recurrent networks that produce spatially structured activity pat- terns for generating neural spike data, outline the proposed GNN-based network inference approact and detail key modifications of the inference model across different generative model configurations.
3.1 Generative recurrent network models
We simulate structured neural activity using ring network of N neurons with recurrent connectivity defined by a local Mexican-hal profile (Figure [b). This architecture supports the formation of stable, spatially periodic activity patterns under uniform excitatory drive (Figure Two spike generation models are considered: a threshold-crossing model and a linear-nonlinear Poisson (LNP) model. Both models integrate recurrent input and shared feed-forward drive, but differ in their spike emission mechanisms, with the former using deterministic thresholding and the latter using stochastic Poisson sampling. These differences lead to distinct spike train statistics. We fix parameters such that multiple co-active activity bumps emerge, providing a challenging testbed for connectivity inference. Full equations and parameter settings are described Appendix A
3.2 Inference network models
We collect spike data from the generative network models over an 8-minute period with a time step of At 0.1 ms, representing the spike trains as {0, 1) 1 where N denotes the number of neurons and L corresponds to the number of time steps (Figure Пp. bottom). The recorded spike train then processed by a structure learning module to infer the underlying neural circuitry, followed by a spike prediction module that concurrently predicts the activity of multiple neurons.
3.2.1 Structure learning module
The objective of this module is to estimate the pairwise connection strength Wiy for every pair (x1, xz). where x; (x) .x.;...; x! represents the spike train of i-th neuron. To achieve this, we apply a ID convolution fCowiDD with 32 kernels, each having a size of 2T/ At, which is twice the synaptic time constant of the generative model. The convolution uses a stride equal to 20% of T across each input spike train. The resulting feature maps are vectorized along the time dimension and passed through fully connected layer four to produce a reduced diieeedidddiniinalll output embedding vector z; (Figure Ic. top):
7, faur (vec (fComvID(x)) (x;)))
We include batch normalization training stability. Next, we concatenate the embeddings (z,, 2;) for every neuron pair and input them
33
immediately after the ReLU activation function improve
A preliminary version of this work appeared in NewrIPS 2022 Workshop on Symmetry and Geometry in Neural Representations 32
3

---

into a multilayer perceptron (MLP) with two hidden layers of 32 units each to estimate the network connectivity strength (Figure I bottom):
Wij MLP ((2;
The inferred weight matrix RNX, containing all pairwise coupling strengths wij between neurons and j. is assumed to be symmetric and free of self-connections. However, we do not impose rotation invariance on w, a constraint that governs the target recurrent weight strengths the generative models. In other words, neurons in the proposed inference model are not required to share identical outgoing synaptic weights with all other neurons within the ring network.
3.2.2 Spike prediction module
This module aims to predict the future activity of the generative network by modeling the dynamics w), given the latent circuitry w and the past spike history. Specifically, this is done by modeling the sequential probability distribution of spike trains over time:
. x² of interacting neurons, expressed as pe (x*+1xx,,.. ......,
T II 1x
p(x(w):
w)
where x" = xxy) represents the spike activity of all N neurons at time f. This approach not only accounts for the temporal dependencies in neural activity but also integrates the connectivity structure learned by the structure learning module. To this end, we employ a GNN message-passing operation:
h = Senc (x (xit)
(4)
m 2 ([ht, bb]))
(5)
log(A feee (he+1) p pxx ( E+Id w) = Pois(A (8)
(7)
h vb
Σ wiy m h
(6)
JEN()
The first expression computes the initial embedding ht for neuron i by encoding its spike history over the past & (= 2r/At) time steps using the encoder fenc(-). which is functionally equivalent to the ID convolution in Eq h with shared parameters. The message vector m', representing the information transmitted from neuron j to neuron i is then computed by applying the function & to the concatenated current states of both neurons. Afterwards, the neuronal state h updated from h by aggregating incoming messages from neighboring neurons j E Nii), where each message is weighted by the corresponding connectivity strength Wey, the (i,j) entry of w, before being processed through the recurrent network vi.
Finally, the decoding function (-)c()) maps the neuronal states to the log firing rates, log X++1, which in turn determine the rates of inhomogeneous Poisson process responsible for generating spikes at time + 1. Although the Poisson likelihood is formally conditioned only on a fixed-length window x' -f and w, the hidden state of gated recurrent unit y encodes summary of all past activity up to time f. This is achieved through gated mechanisms that integrate temporal information and selectively retain or update relevant features from the spike history, enabling the model to capture full temporal dependencies.
3.2.3 Extensions for external inputs and hidden neurons
To better capture the complexities of biological neural systems, extend spike prediction framework in two key ways. First, we incorporate stimulus-driven embeddings alongside spike history to account for external inputs, enabling alignment with circular variables such as head direction. Second, to handle partially observed networks, we introduce hidden neurons whose embeddings are initialized via interpolation from nearby observed units, permitting full-graph message passing during inference. Full implementation details are provided in Appendix
Experiments
Baselines For our connectivity inference experiments, we benchmark our GNN-based model against three widely used baselines: GLM [. sequential non-negative matrix factorization (seqNMF) [34), and tensor component analysis (TCA) 335 (see Appendix C for details). While the GLM stands

---

WN W
WN WW
True
Interee
Figure 2: Quality of connectivity inference from spike train data generated by fully observed network of 100 neurons. (a) Comparison of the ground-truth (orange) and inferred (blue) weight profiles obtained by the GNN-based inference model. The solid blue line represents the average inferred weights across three trials, each initialized with different random seed, with the shaded blue region indicating =1 standard deviation. The inset at the bottom right shows the full inferred weight matrix W. Each row corresponds to inference results from spike data simulated using the greshold-crossingg model (top) and the LNP model (bottom). (b-d) Subsequent columns present the inference quality of baseline methods: (b) GLM, (c) seqNMF, and (d) TCA.
for its ability to predict neural activity, particularly in sensory systems, it also effectively infers coupling effects among neurons. A key distinction between our framework and the GLM lies in how connectivity is represented and learned; we provide detailed discussion of this difference in Appendix D In contrast, TCA and seqNMF are not specifically designed for connectivity inference. Instead, they primarily aim to extract low-dimensional representations and capture neural dynamics. Nonetheless, we use them to identify low-dimensional neuron factors and examine their correlation structures, which can serve as a rough proxy for network connectivity.
Datasets We use both synthetic and real datasets within a unified simulation and evaluation framework. For synthetic data, we generate spiking activity from a 100-neuron ring network described in Section 3.1 simulating 8 minutes of activity at 0.1 ms resolution (4.8 million time steps). This setup applied consistently across all synthetic experiments, including those with and without external inputs, varying recurrent connectivity structures, and under full or partial observability. For real data, we use publicly available HD recordings and motion tracking data from freely moving mice in an open- open-field field environment [36]. The HD trajectories are treated as external inputs, while spike trains from HD cells are used evaluate inference performance within the same framework. All datasets are partitioned into 80% training, 10% validation, and 10% testing splits. See Appendix E for details on preprocessing and integration of real data.
Training In our model, the optimization is framed as minimizing the Poisson negative log-likelihood, where the firing rate X; governs the likelihood of observed spike activity x/. Since is not free parameter but latent representation derived from the observed neural activity via deterministic transformations in the structure learning module, the optimization performed solely over the model parameters O. Given this dependence, w is not explicitly optimized but is instead updated indirectly as θ is optimized. The final objective function for training given by:
T
+ arg min
ΣΣ - log X1)
(A
and We train the model using the Adam optimizer with a learning rate of 5 x 10-4, and further make use of an exponential decay schedule in the learning rate.
Metrics To evaluate connectivity inference, we align weight vectors across neurons to a common phase and normalize for global scaling before computing the normalized inference error, A, relative to ground truth. This accounts for rotational symmetry and scale ambiguity inherent in the network structure. In parallel, assess spike prediction performance using gee-likelihood-basee metric, Lbps. that quantifies improvement over a homogeneous Poisson model, yielding an interpretable score in bits per spike. Full derivations, alignment procedures, and implementation details are provided in Appendix
5

---

Table 1: Performance metries for connectivity inference and spike prediction in fully observed network without external inputs. Results are averaged over the three trials shown in Figure 2
Table 2: Evaluation of inference and prediction accuracy for fully observed network under ex- ternal input conditions. Each value the average over three independent trials.
Thresh.
LNP
Thresh.
LNP
A
Lig
A
Lup=
A 1
Ligu
A 4
Liqu 1
GNN
0.061
0.882
0.049
0.876
GNN
0.073
0.916
0.058
0.924
GLM
0.244
0.695
0.238
0.712
GLM
0.259
0.724
0.245
0.748
SCqNMF
0.789
0.796
SCqNMF
0.791
0.794
TCA
0.762
0.761
TCA
0.760
0.762
4.1 Fully Observed Network
We begin by assessing the accuracy of connectivity inference when the spiking activity of all neurons in the generative network is fully observed. A key finding is that all baseline methods tend to mistakenly infer connections to neurons that are either unconnected or only weakly connected. This is evident from the side dips in the weight profiles and the presence of multiple off off-diagonal stripes in the inferred weight matrices (Figures (FiguresZb--il 2d). Such systematic inference errors stem from overestimated connections driven by strong correlations in neural activity [31, which arise from the global activity patterns intrinsic to our recurrent generative networks (see Figure Tb). Among the baselines, the GLM achieves the closest match the true connectivity but still struggles to properly explain away these spurious correlations. In contrast, our proposed GNN inference model effectively suppresses these artifacts (Figure 2), resulting in significantly lower inference errors (Table This improvement is accompanied by superior spike prediction accuracy (Table TT suggesting that the GNN-based spike prediction module is expressive enough to capture and replicate the underlying dynamics of the recurrent network.
4.2 Fully Observed Network with External Input
We next examine the quality of connectivity inference in fully observed network subjected to external inputs. specifically synthetic, continuously varying cues designed to mimic structured angular modulation. The external drive in this setup consists of low-amplitude signal that fluctuates within the circular space [0, 2xr]. introducing gradual angular shifts in the shared input received by all neurons. The purpose of this design is to introduce structured, non-random external stimulus capable of steering the network's s activity, and determine how such external stimulus influences spike predictability and connectivity inference accuracy.
To implement this, each neuron is assigned preferred direction, arranged uniformly along circular axis. A neuron's 's preferred orientation determines the direction toward which its outgoing synaptic weights are biased. The synaptic weights defined in Eq. (10) are shifted according to the external input ((), effectively rotating the weight matrix to align with the input direction (see Appendix G for further details). This causes the internally generated activity bumps to follow the input stimulus in synchrony, producing a smooth, input-driven trajectory across the neural manifold (Figure 3p, top row).
Empirically, we find that introducing this structured external cue leads to improved spike predictability compared to the generative network without external input, which exhibits spontaneous drift in its global activity pattern (Table 2). This suggests that networks driven by such input produce more predictable dynamics. However, despite this increase in spike predictability, connectivity inference errors remain similar or slightly elevated. This likely arise because external drive dominates spike timing, reducing the relative explanatory power of the inferred recurrent weights. Meanwhile, all baseline inference methods continue to exhibit persistent artifacts (Figure B3. top row), such as spurious correlations between weakly unconnected units. In contrast, our GNN-based method maintains a clear advantage, delivering more accurate reconstructions of the ground-truth network across trials (Table 2 This highlights the robustness of our approach even in settings where external inputs strongly shape network activity.
4.3 Weakly Correlated, Fully Observed Network with External Input
We next explore the impact of network activity correlation structure connectivity estimation by comparing two distinct weight profiles: the local Mexican-hat profile (Figure [3p, top), which produces multiple periodic activity bumps as in Figure 3p (top), and modified profile characterized by localized excitation at zero angular difference combined with broadly tuned inhibition (Figure

---

Weight Profie I
W
Cell Index
Weight ex Profie II
Call Indox
Time (min)
Figure 3: Evaluating connectivity inference. Top: bump dynamics in a ring attractor network driven by external input. Bottom: dynamics under a modified recurrent connectivity profile. (a) Transition from a local Mexican-hat hat profile (top) to new configuration with local excitation and broadly tuned inhibition (bottom). (b) Spike raster plot showing rotating activity bumps induced by sinusoidal external inputs. (c) Comparison of ground-truth (orange) and inferred (blue) weight profiles estimated using GNN, GLM, seqNMF, and TCA, following the same order as the corresponding sub-figures.
bottom). This adjusted profile ensures the formation of only a single stable activity bump, effectively simulating conditions typical of the HD system (Figure 3p, bottom). We drive the ring network with sinusoidal input signal that shares the same period as in Figure 3 (top), and then investigate how reducing spurious correlations in spike-train data by altering connectivity structures affects the accuracy of connectivity inference.
When shifting from the Mexican-hat to the localized excitation profile, anticipate notable reduction in spurious correlations among neural activities, as synchronized activity patterns diminish. This reduction is expected to lower inference errors by minimizing reliance on correlated noise and emphasizing true connectivity- driven structure. Consistent with this expectation, all baseline methods show improved inference performance under this con- dition (Table 3). benefiting from fewer misleading correlations. Notably, our GNN-based approach achieves even higher accuracy and robustness (Table 3 and Figure 3.. bottom P). These results underscore the importance of the correlation structure in shaping inference quality and reinforce the effectiveness of our method across both strongly and weakly correlated activity conditions.
Table 3: Connectivity inference er- ror A and spike prediction perfor- mance Laps for modified input profile featuring localized excitation and broadly tuned inhibition, in a weakly correlated, fully observed network with external input. Val- ues are averaged over three trials.
LNP
A(I)
(iga(1)
(.,
GNN
0.048
2.652
0.043
2.668
GLM
0.125
2554
0.117
2.576
0.374
0.378
TCA
0,362
0330
4.4 Weakly Correlated, Partially Observed Network with External Input
To separate the challenges of inference in partially observed networks from those arising in strongly recurrent circuits, we have thus far focused on a fully observed setting. We now examine how varying the number and ratio of observed and hidden neurons influence spike prediction accuracy and circuit inference error within a partially observed setting, using the simulated HD network described in Section 4.3 Specifically, from a generative ring network comprising 100 neurons, we randomly select No {60, 80, 100} observed neurons, while the remaining neurons serve as the pool for selecting N₂ E {0, 5, 10, 20) hidden neurons. The spiking activity of the observed neurons used as training data, and hidden neurons participate in the model as graph nodes without direct training input. By systematically varying N, and N₁, we assess the impact of different configurations on circuit inference quality. When all 100 neurons are observed, hidden neurons are not included, as the total network size remains fixed at 100.
neuron
We begin by inferring synaptic connectivity among No observed and N₁. hidden neurons, with each neuron assigned an angular position uniformly spaced around a ring. Although the initial inference yields connection weights between all these neurons, retain only the weights between observed neurons, organized into a partially complete N x N weight matrix indexed by angular positions (Figure Ap). This choice reflects the transductive nature of the task: only the observed
"In Fig. 3. the error hars in the bottom row wider than those in the top. This can he explained as follows. Weight Profile is broad, smoothly varying inhibitory Gaussian. Because tail is nearly flat, small variations in frequency coefficients have little effect on predicted spikes patterns. As a result, the likelihood surface is flatter in those directions, leading to greater variability in parameter estimates across runs. Moreover, the stronger and more uniform inhibition of Profile II suppresses activity across much of the population. substantially reducing the number of informative spikes. Together, the weaker sensitivity and smaller sample of observations yield greater parameter uncertainty, which manifests as the wider error bars in the bottom row.
7

---

Interned Weight Naatts
Interpoleted Weight Matti
True
lotsenveed neurore onlvl.
0.47)
140
Coll
240
pobllde Prcte
d
(A - 0.03)
g
1.00
The Weight Profil
240
whywHry
1.10
1.00
podaledd
23
Figure 4: Evaluation of inference performance in a partially observed network as function of the total number of neurons used for inference. (a) Ground-truth synaptic weight matrix with a smooth, spatially structured profile on ring. (b) Synthetic weight matrix among observed neurons only (N,, 80). The matrix is not inferred from spike train data but constructed by masking unobserved rows and columns of the ground-truth matrix in (a), followed by the addition of small uniform noise to mimic inference uncertainty. (c) Full weight matrix after two-dimensional interpolation, showing partial recovery of spatial structure. (d) Control: shuffled version of the inferred matrix, where either rows or columns are randomly permuted to disrupt spatial organization while preserving marginal distributions. (e) Interpolated version of the shuffled matrix exhibits degraded structure and higher error (A 0.53) compared to the unshuffled case (A- = 0.27). (f) Spike prediction accuracy (Chps) improves as more neurons incorporated, with colored curves representing different numbers of observed neurons (N, = 60, 80, 100). Error bars denote one standard deviation. (g) Circuit inference error (A) plotted against the total number of neurons used for inference. Small insets show interpolated weight profiles for selected configurations, revealing how the structure quality varies with conserve-ttohhdddee neuron ratios.
neurons have spiking activity available during training, making their inferred interactions empirically grounded. In contrast, weights involving hidden neurons rely solely model assumptions or priors and therefore carry higher uncertainty. To reconstruct the full N x N weight matrix, including connections involving hidden neurons, apply two-dimensional linear interpolation. This assumes that connectivity varies smoothly along the ring. To preserve the circular topology, the matrix is temporarily extended along the angular dimension. Interpolation is then applied row-wise (for outgoing connections) and column-wise (for incoming connections), filling in missing values using nearby known weights (Figure 4..
Our analysis shows that increasing the number of observed neurons generally enhances spike pre- diction accuracy and concurrently reduces inference error (Figures 4g- indicating more precise recovery of latent dynamics when a larger fraction of the network is directly observed. When hidden neurons are added, spike prediction accuracy continues to improve across all configurations, though the gains diminish as the number of hidden neurons increases (Figures(4f). In contrast, inference error initially decreases but then saturates or slightly increases, particularly when the hidden-to-obsers neuron ratio becomes large (Figure 4g). For example, with fixed number of observed neurons, expanding the hidden population from to 20 yields diminishing returns in spike prediction accuracy and can lead to plateau or rise in inference error. This suggests that while hidden neurons can support better spike prediction by capturing latent dynamics, they may also introduce structural ambiguity, especially when weakly constrained by observed activity. These trends highlight a tradcoff' between functional prediction and structural inference, and suggest that optimal performance does not necessarily result from maximizing the number of hidden neurons.
Finally, to assess whether the inferred weight matrix reflects meaningful structure beyond chance, we shuffled its rows while preserving their marginal distributions. This disrupts any spatial alignment while maintaining the local weight statistics. As shown in Figures fi.e, the resulting interpolated matrix displays no coherent structure and yields substantially higher inference error (A 0.53) compared to the unshuffled case (A = 0.27), supporting the conclusion that the original inferred weights capture nontrivial spatial patterns not attributable chance.
8

---

4.5 Head Direction Cell Network
We have demonstrated that it is possible to infer aspects of neural circuitry even when training spike data covers only a subset of the full neural population. Can this modeling framework be extended to real-world scenarios, such analyzing neural ensembles recorded during spatial navigation or other behaviorally relevant tasks? To investigate this, we apply our circuit inference model to a publicly available dataset [6 of 19 simultaneously recorded HD cells in the anterodorsal thalamic nucleus (ADn) of freely moving mice (see Appendix E for experimental details), In this setting, we fix the total number of neurons the ring network N 100, with No = 19 observed neurons and NA. E {0, 5, 10, 2 unobserved. Although the actual circuit size is unknown, choosing N > No allows us approximate the underlying connectivity with an angular resolution of 2 2T/N.
To apply the inference model, we first construct the tuning curves of the 19 observed neurons and estimate their preferred head direc- tions (Figure Op. Each observed neuron is then matched to one of the N positions in the ring network by assigning it to the neuron whose preferred head direction, spaced at intervals of 2T/NN, is closest to the observed tuning peak. After this assignment, N₁, hid- den neurons are randomly selected from the remaining positions in the ring network to initialize the proposed model. Following the same training procedure as in Section 4.4] we find that the learned weight profiles consistently exhibit similar patterns across varying hidden neuron counts (Figure 5 This consistency points to a connectivity motif characterized by local excitation and sur- rounding inhibition, echoing the structure proposed in continuous attractor network models of HD cells [12.. Although some vari- ability in the inferred weights is observed, likely due to the limited fraction of observed neurons relative to the full network, the re- sults demonstrate that framework can effectively recover key features of the underlying circuit even under partial observability in real neural recordings.
W
W
Figure 5: Inferred weight profiles derived from 19 real HD cells in the anterodorsal thalamic nueleus (ADn) of mice, with the number of observed neurons fixed at N. 19 and the number of hidden neurons (N₂, increasing from 0 to 20. The weights are obtained through linear interpolation under the assumption that the total number of neurons is set 100.
5 Discussion
In this paper, we presented an self-supervised approach for inferring neural circuit connectivity from population spike activity. By designing two functionally specialized network modules, one for learning synaptic connectivity and the other for predicting future spiking, we established a link between the latent representations in the structure learning module and the underlying network connectivity. Using both simulated and real neural datasets, we demonstrated that our Gaseb model, which captures the dynamics of interacting neurons, accurately recovers true weight profiles and performs favorably in comparison to traditional approaches such as GLMs and other baselines that rely approximating activity correlations. We further found that our method remains effective under diverse conditions, including varying levels of activity correlation, partial observability, and the presence absence of external inputs. Additionally, our analysis of real neuronal data highlights the practical applicability of this approach. Specifically, applying our model to HD cell spike recordings from awake animals revealed network structures consistent with the ring attractor hypothesis proposed in earlier mechanistic studies [2].
A natural question is how the proposed framework extends beyond ring attractor networks to neural circuits with different topological structures or functional properties. The reliance on ring geom- etry in our current experiments arises solely from the optional geometric feature provided to the spike prediction module. The self-supervised learning objective of predicting future spikes while jointly estimating a latent connectivity matrix does not depend on this assumption. In the present implementation, the message function of-) receives only the concatenated node embeddings and therefore imposes no explicit ring metric. When geometric priors are desirable, an additional attribute such as one-dimensional distance (chains), two-dimensional Euclidean or toroidal distance (cortical sheets, grid cells), or learned positional embedding can be incorporated into the message function. Importantly, the choice to include or omit such features toggles the inductive bias without altering the loss, dooulm architecture, or optimization routine. This flexibility allows the same inference network to generalize across diverse circuit topologies, including those with fundamentally different structural or functional organization.
9

---

Recent approaches that couple spike modeling with latent graph estimation fall into two primary categories. NetFormer [37] encodes each neuron's recent spike history a token and derives a step-wise attention matrix whose entries are interpreted as couplings, thereby updating the graph at every time step. Because attention weights are used directly for prediction, NetFormer lacks a dedicated messagg-ppssing module and instead recomputes pairwise interactions at each step. This design precludes iterative propagation through a learned weight matrix and limits the ability to filter indirect correlations. Notably, Fig. 15 in the original NetFormer paper reports inference performance on exactly the same benchmark considered here, and its accuracy is drastically worse than that achieved by our framework. By contrast, AMAG [38 employs message passing but initializes the adjacency matrix from random weights or correlation than inferring connectivity. Our framework differs in its explicit separation of structure learning and spike prediction: a latent connectivity block continuously estimates circuit structure from activity, while a dedicated pressage-passingg block leverages this estimate for prediction. This division yields more stable graph inference, naturally accommodates unobserved neurons through auxiliary nodes, and sustains predictive accuracy in recurrent regimes.
heuristics, focusing on refining rather
While these results are promising, we acknowledge a few limitations of the proposed model. First, our framework assumes that the underlying network connectivity remains fixed throughout the observation period. However, in biological neural systems, synaptic connections are often plastic and can evolve over time in response experience or changing behavioral demands. A promising direction for future research would be to extend the model to infer time-varying connectivity, allowing for the reconstruction of dynamically changing weight profiles that better reflect the adaptive nature of real neural circuits. In addition, while our experimental setup captures the transition in the generative recurrent network from exhibiting spatially periodic activity forming single activity bump on the ring, this initial configuration with multiple bumps can be interpreted as a simplified representation of the grid cell (GC) system. Specifically, resembles a scenario in which an animal moves along a direction aligned with one of the principal lattice vectors of 2D virtual triangular lattice, thereby periodically activating every vertex along that path [39]. natural extension of this work would be to apply our inference framework to 2D continuous attractor networks, both in simulation and using real grid cell data, to further explore its utility in inferring neural circuitry in more complex spatial systems. With the advent of nulti-Neuropsielss probes and mesoscale two-photon calcium imaging, is now possible to record from hundreds to thousands of neurons across extended spatial domains. These advances open the door to deploying approach MEC grid cells and other spatially orgamized circuits with unprecedented coverage, providing a more comprehensive testbed for connectivity inference in large-scall neural populations. Finally, we acknowledge that similar neural dynamics can arise from multiple distinct circuit configurations [40] 41 31 Our model aims to infer one possible configuration that is consistent with the observed spiking activity, but we do not claim that the inferred circuitry is unique. Integrating tools from algebraic topology or combinatorial network theory might offer new avenues for characterizing the space of circuits compatible with observed neural activity.
Overall, our approach offers flexible and interpretable framework for inferring latent neural connec- tivity from spiking data in self-supervised and data-driven manner, By bridging structure learning with predictive modeling, we provide a general method for uncovering underlying circuit dynamics that can complement and extend existing tools in computational neuroscience.
10

---

Acknowledgments
K.Y. is grateful to lla Fiete and Abhranil Das for their helpful discussions during the early stages of this research, and to Xaq Pitkow for providing foundational guidance during the postdoctoral work with him. K. Y. is supported in part by the National Research Foundation of Korea (NRF) grant (No. RS-2024-003370992) the Institute of Information & communications Technology Planning & Evalu- ation (IIT)) grants (No. RS-2020-11203773 Artificial Intelligence Graduate School Program; No. (22002)--RS-20230025533333 14, Artificial Intelligence Semiconductor Support Program (Hanyang University)) funded by the Korean government (MSIT), and in part by Samsung Electronics Co., Ltd.
References
[1] 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 a complete neuronal population. Nature, 454(7207):995-999 2008
[2] Elad Schneidman, Michael J Berry, Ronen Segev, and William Bialek. Weak pairwise corre- lations imply strongly correlated network states in a neural population. Nature, 440(7087): 1007-1012, 2006.
[3] Karl Friston. Functional and effective connectivity: review. Brain Connectivity, 1(1):13-36. 2011.
[4] Ari Pakman, Jonathan Huggins, Carl Smith, and Liam Paninski. Fast state-space methods for inferring dendritic synaptic connectivity. Journal of Computational Neuroscience, 36:415-443, 2014.
[5] Daniel Soudry, Suraj Keshri, Patrick Stinson, Min-hwan Oh, Garud Iyengar, and Liam Paninski. Efficient" shotgun" " inference neural connectivity from highly sab-sampled activity data. PLoS Computational Biology, ([10))e1004464 2015.
[6] Yury V Zaytsev, Abigail Morrison, and Moritz Deger. Reconstruction of recurrent synaptic connectivity of thousands of neurons from simulated spiking activity. Journal of Computational Neuroscience, 39:77-103, 2015.
[7] Mikkel Elle Lepperød, Tristan Stöber, Torkel Hafting, Marianne Fyhn, and Konrad Paul Kording. Inferring causal connectivity from pairwise recordings and optogenetics. PLoS Computational Biology, 19(1 ((11))e100 :)11e1111111 1574, 2023.
[8] Alfred Gierer and Hans Meinhardt. A theory of biological pattern formation. Kybernetik, 12: 30-39, 1972.
[9] André-Joseph Koch and Hans Meinhardt. Biological pattern formation: from basic mechanisms to complex structures. Reviews of Modern Physics, 66(4):1481, 1994.
[10] François Schweisguth and Francis Corson. Self-organization in pattern formation. Developmen- tal Cell. 49(5):659-677 2019.
[11] Rani Yishai, R Lev Bar-Or, and Haim Sompolinsky. Theory of orientation tuning in visual cortex. Proceedings of the National Academy of Sciences, 92(9):3444-3844 1995.
[12] Kechen Zhang, Representation of spatial orientation by the intrinsic dynamics of the head- direction cell ensemble: a theory. Journal of Neuroscience, (666):2112-2126, 1996.
[13] H Sebastian Seung. How the brain keeps the eyes still. Proceedings of the National Academy of Sciences, 93393441333911344, 1996.
[14] Mark C Fuhs and David S Touretzky. A spin glass model of path integration in rat medial entorhinal cortex. Journal of Neuroscience, 26(16):4266-4277 2006.
[15] Yoram Burak and lla R Fiete. Accurate path integration in continuous attractor network models of grid cells. PLoS Computational Biology, 100)21, 2009.

---

[16] Thomas Kipf. Ethan Fetaya, Kuan-Chich Wang, Max Welling, and Richard Zemel. Neural relational inference for interacting systems. In International Conference on Machine Learning. pages
PMLR, 2018.
[17] Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia. Learning to simulate complex physics with graph networks. In International Conference Machine Learning, pages 8459 8468. PMLR, 2020.
[18] Victor Bapst, Thomas Keck, A Grabska-Barwinskar Craig Donner, Ekin Dogus Cubuk, Samuel Schoenholz, Annette Obika, Alexander WR Nelson, Trevor Back, Demis Hassabis, et al. Unveiling the predictive power of static structure in glassy systems. Nature Physics, 16(4): 148-454, 2020.
[19] Yuyang Wang, Zijie Li, and Amir Barati Farimani. Graph neural networks for molecules. In Machine Learning in Molecular Sciences, pages 21-66. Springer, 2023.
[20] Thomas N. Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. In International Conference Learning Representations, 2017.
[21] Will Hamilton, Zhitao Ying. and Jure Leskovec. Inductive representation learning large graphs. Advances in Neural Information Processing Systems, 30, 2017.
[22] David J Thouless, Philip W Anderson, and Robert G Palmer. Solution of 'solvable model of a spin glass'. Philosophical Magazine, 35())5993-601 1977.
[23] Yasser Roudi, Joanna Tyrcha, and John Hertz. Ising model for neural data: model quality and approximate methods for extracting functional connectivity. Physical Review E -Statistical, Nonlinear, and Soft Matter Physics, 79(5)0511915, 2009.
[24] Jascha Sohl-Dickstein, Peter Battaglino, and Michael R DeWeese. Minimum probability flow learning, arXiv preprint urXiv:0906.47799 2009.
[25] Jascha Sohl-Dickstein, Peter B Battaglino, and Michael R DeWeese. New method for parameter estimation in probabilistic models: minimum probability flow. Physical Review Letters, 107 (22):220601, 2011.
[26] Su-In Lee, Honglak Lee, Pieter Abbeel, and Andrew Y Ng. Efficient & regularized logistic regression. In AAAI, volume 6, pages 401-408, 2006.
[27] Pradeep Ravikumar, Martin J Wainwright, and John D Lafferty. High-dimensional ising model selection using errruuulardeed logistic regression. The Annals of Statistics, pages 1287-1319, 2010.
[28] John Ashworth Nelder and Robert WM Wedderburn. Generalized linear models. Journal of the Royal Statistical Society Series A: Statistics in Sociery, 135(3):370-385 37,-444,, 1972.
[29] Wilson Truccolo, Uri 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 Nearophysiology, 93(2):1074-1089 2005.
[30] Jonathon Shlens, Greg D Field, Jeffrey L Gauthier, Matthew I Grivich, Dumitru Petrusca, Alexander Sher, Alan M Litke, and EJ Chichilnisky. The structure of multi-neuron firing patterns in primate retina. Journal of Neuroscience, 26(32):8254 8266, 2006.
[31] Abhranil Das and lla R Fiete. Systematic errors in connectivity inferred from activity in strongly recurrent networks. Nature Neuroscience, 23(10):: 2020.
[32] Tachoon Park. JuHyeon Kim, DongHee Kang, and Kijung Yoon. Graph neural networks for connectivity inference in spatially patterned neural responses. In NewrlPS 2022 Workshop Symmetry and Geometry in Neural Representations, 2022.
[33] Sergey loffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In International Conference on Machine Learning, pages 448-456. pmlr, 2015.
12

---

[34] Emily L Mackevicius, Andrew H Bahle, Alex H Williams, Shijie Gu, Natalia Denisenko, Mark Goldman, and Michale Fee. Unsupervised discovery of temporal sequences in high-dimensional datasets, with applications to neuroscience. Elife, 8:e38471, 2019.
[35] Alex H Williams, Tony Hyun Kim, Forea Wang, Saurabh Vyas, Stephen I Ryu, Krishna V Shenoy, Mark Schmitzer, Tamara G Kolda, and Surya Ganguli. Unsupervised discovery of demixed, low-dimensional neural dynamics across multiple timescales through tensor compo- nent analysis. Neuron, 1896): 1119, 2018.
[36] Adrien Peyrache, Marie M Lacroix, Peter C Petersen, and György Buzsáki. Internally organized mechanisms of the head direction sense. Nature Neuroscience, 18(4):569-575, 2015.
[37] Ziyu Lu, Wuwei Zhang, Trung Le, Hao Wang, Uygar Sümbül, Eric Todd SheaBrown, and Lu Mi. Netformer: An interpretable model for recovering dynamical connectivity in neuronal population dynamics. In The Thirteenth International Conference on Learning Representations, 2025.
[38] Jingyuan Li, Leo Scholl, Trung Le, Pavithra Rajeswaran, Amy Orsborn, and Eli Shlizerman. Amag: Additive, multiplicative and adaptive graph neural network for forecasting neuron activity, Advances in Neural Information Processing Systems, 36:8988-9014, 2023.
[39] Kijung Yoon, Sam Lewallen, Amina A Kinkhabwala, David W Tank, and Ila R Fiete. Grid cell responses in ld environments assessed as slices through a 2d lattice. Neuron, 196) 2016.
[40] Astrid Prinz, Dirk Bucher, and Eve Marder. Similar network activity from disparate circuit parameters. Nature Neuroscience, 7(12):1345-1352, 2004.
[41] Carina Curto and Katherine Morrison. Relating network connectivity to dynamics: opportunities and challenges for theoretical neuroscience. Current Opinion Neurobiology, 58:11-2: 2019.
[42] James Knierim and Kechen Zhang. Attractor dynamics of spatially correlated neural activity in the limbic system. Annual Review of Neuroscience, 35(1):267-285, 2012.
[43] Daniel D Lee and H Sebastian Seung. Learning the parts of objects by non-negative matrix factorization. N (01(6755))78887771 1999.
[44] Liam Paninski, Eero Simoncelli, and Jonathan Pillow. Maximum likelihood estimation of stochastic integrate-and-fir neural model. Advances in Neural Information Processing Systems, 16, 2003.
[45] Liam Paninski, Shy Shoham, Matthew R Fellows, Nicholas G Hatsopoulos, and John P Donoghue. Superlinear population encoding of dynamic hand trajectory in primary motor cortex. Journal of Neuroscience, 24(39):8551--8) 116611111,11 2004.
[46] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high- perhoorrform deep learning library. arxiv 2019. arXiv preprint arXiv:1912.01703. 10. 1912.
[47] Matthias Fey and Jan E. Lenssen. Fast graph representation learning with PyTorch Geometric. In ICLR Workhhop Representation Learning Graphs and Manifolds, 2019.
13

---

NeurIPS Paper Checklist
Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the
paper's contributions and scope?
Answer: [Yes]
Justification: There are four main claims layed out in the abstract: joint inference of connectivity and spike prediction using a GNN, effective resolution of spurious correla- tions, handling of unobserved neurons via auxiliary nodes, and alignment with theoretical predictions on real data. All of these points are addressed in the main text.
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 perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize to other settings.
It is fine 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 of the work performed by the authors?
Answer: [Yes]
Justification: The main limitations of the work, well as future directions that might address some these limitations, are layed out in the discussion portion the paper.
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 on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be 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 low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to 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, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an 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: This is not theoretical paper.
Guidelines:
The answer NA means that the paper does not 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 in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix 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 explain our model in Section experimental settings with data in Section 4 and hyperparameters in Appendix A
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not 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 make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to model. In general. releasing code and data often one good way to accomplish this, but reproducibility can also be provided via detailed 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, or other means that are appropriate to the research performed.
specific model and empirical evaluation, may
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) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (c.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 or instructions for how to construct the dataset).
(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 possible for other researchers to have some path to reproducing verifying the results.
15

---

reproducible, they 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, described in supplemental material?
Answer: [No]
Justification: All experiments performed are simple enough reproduce without access to code. Code is available upon request.
Guidelines: While open access to data and code is not provided, all experiments in the paper are straightforward to reproduce. The head direction (HD) cell recordings used in the real data analysis are publicly available, and the code for the generative recurrent network used produce synthetic data is also publicly accessible and easy to modify. Additionally, the full codebase used for training and evaluation available upon request.
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips. ec/
Sublic/guidess /ttpg//dd///eeuummmiiimem
for more details.
While encourage the release of code and data, we understand that this might not 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 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 wwwoodeummisssoooo...oo
7 //nips wttpccccccpppccccccececcooocccce wttp:/iccggiieee/eeeeeooe
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. If only subset of experiments should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible supplemental material (appended the paper) is recommended, but including URLs to data and code is 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: Details and explanations data split are summarized in Section
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, 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]
Justification: We show the standard error in most evaluation results, with an average of over three random seeds. In addition, we assess whether the inferred weight structure reflects meaningful spatial organization beyond chance via shuffling the weight matrix.
Guidelines:
The answer NA means that the paper does not include experiments.
16

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or 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).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report I-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 is 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 in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
This is discussed in Appendix H
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 it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics http::::/eerrrp..ccpppp....///// http:///cccpceeeecccccccciieeeee
Answer: [Yes]
Justification: We have read and understood the code of ethics; and have done our best to conform.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there special consid- eration due to laws or regulations in their jurisdictio).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: Our work involves small models and datasets. does no impact the society at large, beyond improving understanding of certain aspects of neuroscience.
17

---

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 or 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 considcrations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to 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 to 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 conrectly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional 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 to 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 deseribe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Since such models or datasets involved, our work poses no risk of misuse.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-uss 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 is 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
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?
Question:
Answer: [Yes]
Justification:
We use the latest versions of both PyTorch and PyG, and cite both in Appendix
H
Guidelines:
The answer NA means that the paper does not 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.
18

---

The name of the license (e.g., CC-BY 4.0) should be included for each asset. For scraped data from a particular source (c.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, pappttss.pppssiiipp.. paperswithcode thcode. com comddatasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors 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:
[NA]
Justification:
No such assets are introduced.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/model as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an 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 to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
This work does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information 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 the main paper.
According to the NeurIPS Code Ethics, workers involved data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
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 an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
We have no human participants in our study.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
19

---

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 to 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.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification:
Guidelines:
This research does not involve L.LMs as any original components.
The answer NA means that the core 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 not be described.
20

---

A Generative Recurrent Network Models
We consider a ring network composed of N neurons (Figure [T], where the synaptic weights W sy from each neuron to all others in the ring follow a rotation-invarian local Mexican-har profile:
Wi = - - ae -22
(10)
This structure represents a difference of Gaussians, where dy denotes the distance between neurons i and j. and D1 and 02 are the standard deviations of the narrower and broader Gaussians, respectively, with σi < Jz.. The parameter a, slightly greater than one, transforms local excitation into weak inhibition. This inhibitory effect facilitates pattern formation under uniform feed-forward excitatory drive and maintains dynamical stability [15]. With this ring network architecture established, we explore two distinct spike generation processes.
Threshold-crossing model In this model, the input to each neuron time step is determined by weighted sum of synaptic activations from neighboring neurons, modulated by feed-forward inputs. Specifically, the input vector g(t) for all neurons is given by
g(t) (t) + b E(t))
(1I)
where s(t) E RN represents the synaptic activations of N neurons, and W E ERNN is the recurrent connectivity matrix defined in Eq. [ The term b represents a uniform excitatory drive, and {(t is a white Gaussian noise term per neuron, with zero mean and a specified standard deviation, leading to Poisson-like variance proportional to the mean activation. The ratio of recurrent to feed-forward input controlled by the recurrent weight strength parameter r. In this recurrent network model, a neuron i emits spike when its input g: (t) at time step é exceeds a threshold gth. The synaptic activation of spiking neurons incremented by I, while the activation of other neurons decays exponentially with time constant т, following the equation:
ds(() + s(t) 0 (g(() - 9thh)
=(t)=
(12)
dt
Here, (()) is the Heaviside step function, producing a binary vector of spikes across the network. The simulation is performed using a discrete-ting update with a step size of At. The recurrent weight strength parameter regulates the extent of network activity correlations. When is small, feed-forward noise dominates, leading to relatively uncorrelated activity, whereas larger values of result in globally structured patterns of periodically spaced activity bumps (Figure top). The chosen value of r = 0.025 ensures that unconnected neurons in different co-active bumps exhibit strong correlations, making the spiking data an ideal benchmark for evaluating the performance of the proposed connectivity inference method.
Linear-nonlinear Poisson model The inputs to the linear-nonlinear Poisson (LNP) model are computed similarly to Eq, M.). with the key difference that the feed-forward input remains a constant b across all neurons, without any additive noise. This is because stochasticity in the LNP model is inherently introduced through an inhomogeneous Poisson process:
x(()) Ao ReLU [g: 9th (t) ~ Pois (t))
(13)
(14)
where the firing rate (t(t) of neuron is obtained by applying a rectifying nonlinearity (ReLU) to the summed input (t), shifted by threshold 9th- This firing rate then govems inhomogeneous Poisson process, which determines the number of spikes (t,(t) produced by neuron i at time t. The underlying neural dynamics remain the same in the previous model, except that synaptic activations can now be incremented by values greater than 1, reflecting the Poisson-distribured spike count:
ds(t) dl
+ s(t) = nt)
(15)
This formulation allows for variability spike counts at each time step, in contrast to the threshold- crossing model, where only a single spike could be emitted per neuron per time step. The exact parameter values used in both spike generation processes are listed in Table 4 and further details can be found in Das and Fiete [31].
14

---

Table 4: Model parameters of generative recurrent networks
Parameter
Description
Value
N
number of neurons
102
(11, 02)
s.d. of two Gaussians for symmetric center- weights
(6.98, 7.00)
scalar of the second Gaussian
1.0005
b
uniform excitatory drive
10-3
r
recurrent weight strength of threshold-crossing model
2.5 x 10 10-2
0g
s.d. of white Gaussian noise
3 x 10-1
9th
threshold of spiking process
7.35 x 10
synaptic time constant discretization step size
10 ms
At
0.1 ms
Generative Network Parameters The generative model used in this study corresponds to the highly structured ring attractor network operating in the strongest recurrent weight regime, as characterized in Das and Fiete [1. The network consists of N 100 neurons, with structured recurrent connectivity shaped by a symmetric center-surround profile defined by two Gaussian components. The standard deviations of the two Gaussians are 0 6.98 and 02 7.00, with a relative amplitude scaling factor a 1.0005. Each neuron receives a uniform excitatory drive = 10-3, Recurrent interactions are governed by a coupling strength
2.5 10-2, consistent with the threshold- crossing model formulation. Neural activity is perturbed by white Gaussian noise with standard deviation og 3 x 10-1 Spiking occurs when the membrane potential exceeds the threshold Oth 7.35 10-4, Synaptic dynamics are modeled with time constant 10 ms, and the network is simulated using a discretization step of At 0.1 ms.
B Extensions for External Inputs and Hidden Neurons
B.1 Incorporating External Input
The spike prediction module, as defined by Eqs. observed network without external input drives. Specifically, it is designed for spike data collected from a generative network whose activity pattern undergoes noise-driven drift. However, real neural systems are often influenced by external stimuli. To account for this, our spike prediction module can be readily extended to accommodate conditions where activity patterns are driven by an external input 0:
operates under the assumption of a fully
hi..- = fenc.e (x:- +t) (16) h²,0 = f (0. (0.-88 -bb) b) (17)
[ h(,o]
(18)
Here, hf_ and hi denote the initial embeddings for neuron i's spike activity and the input stimulus at time step f, respectively. f-ne,z(-) essentially equivalent to feac(-) in Eq. while fene ( serves as a separate encoder for 0. Each neuron i is assigned base preferred stimulus value 0, = 2x i where N is the total number of neurons, resulting in uniform coverage over the circular stimulus space [0, 2x). However, because circular variables lack an absolute origin, these preferred values are defined only up to rotational shift. To account for this symmetry and enable alignment any chosen reference point in the external stimulus space, we introduce bias parameter in Eq, (17). This bias applies global rotational shift to the population's 's preferred stimuli, allowing the network to align its internal representation with the external variable 0g. This design is particularly appropriate for experiments involving ring attractors that encode circular or periodic variables such as head direction (HD) or orientation [42] Ultimately, the initial state h of neuron i is formed by concatenating both types of embeddings as shown (q. (88).
B.2 Modeling with Hidden Neurons
Another critical aspect of network inference is addressing the fact that observed neural data often does not capture the entire network's activity which is typically the case in real-world scenarios. To extend our model framework to conditions where unobserved neurons are present in the neural circuit, we introduce hidden neurons into the inference network while accounting for the absence of
produces a periodic activity pattern with multiple bumps, so it does not strictly represent the mechanistic model of the HD system, which typically features a single hump activity.
1 Our generative ring network in Section 3.1
15

---

Algorithm Generalized Linear Model (GLM)
Initialize model parameters: b +- O, w + 0, t-1
2 for i = to N do
for iteration to Niter do
Compute coupling filter: Seouple Bz if then
Circularly shift by i steps end if
Extract neighboring spike history: X_; f- 11:1-1
(EN(()
Compute log firing rate: log A₃ < Semple (X-:w) +
10
Compute loss: C E, E, (Xf - x log X;)
11
Update parameters (b, w, z) < arg mins, w.z L using Quasi-Newton method
12
end for
13 end for
their spike activity. This requires consistent adjustment each module. In the structure learning module, the feature embedding Zj of hidden neuron, j is initialized by linearly interpolating the embeddings of its nearest observed neurons on either side, using the representations derived in (q. Similarly, in the spike prediction module, the initial embedding h of hidden neuron obtaied via interpolation of the same two closest observed neurons' embeddings, derived from Eq. Following this initialization, we update the embeddings of all neurons- both observed and hidden-through! through message- pasginggggnngggg operations over the complete graph structure, while training the model exclusively on the spike activity of the observed neurons.
C Baseline Models
C.1 Generalized Linear Model (GLM)
In this study, we first construct the coupling filter feoupte as a linear combination of set of raised cosine basis functions. Specifically, we define basis matrix B IR x32 where each column represents a distinct raised cosine filter. The coupling filter is then given by Scouple Bz, with denoting the basis coefficients. To compute the log firing rate of neuron i, we project the spike history of its neighboring neurons onto the filter I /couple. weighted by the connection strengths w:
log A₃ = fcouple * (X,ww) +
Here, X-, represents the spike trains from neurons the neighborhood Ntt,, truncated up time L - 1. The symbol denotes convolution over time, capturing the temporal influence of past spikes from neighboring neurons. The resulting signal passed through linear transformation and offset by a bias term b. The parameters b, w, and z are optimized by minimizing the negative log-likelihow through a Duwsi-Newton method:
C= E=
(A$
x log
16

---

Algorithm 2 Tensor Component Analysis (TCA)
Input:
activity tensor X RRNIXK rank R, #ALS steps S 55500
Smooth spikes: M (X(X)
> exponentially-weigheed moving average
Reshape M back to x 7 x K with T = T.KK
4 Randomly initialize factors W = [w U=1u V |
for to do for r-= to R do
Update w(r) + arg miny [1 - Σ, (n(vr)uu(r) ||
Update (r) + arg mina [(2) - Σ, (w
Update ((r) arg minl [X(3) Σ, ruul 1/
10
end for
end for
12 for r to do
13
Z [ &w @w(f))l())@)l))
= Σ (r) Siikt
14
15 end for
16
Return average interaction matrix E
R E,,E' E(r)
C.2
Tensor Component Analysis (TCA)
In our experiment, we recorded neural activity during single continuous session. To apply TCA as described by Williams al. [35], which requires data structured across multiple trials, we partitioned this continuous spike train into K 10 equal-length segments. Each segment was treated as an individual trial, enabling the construction of three-dimensional data tensor E R.xXX'X where N is the number of neurons, T is the number of time points per segment, and K is the number of segments. This tensor was then decomposed using TCA into a sum of rank-one components:
R
Σ
Entk
In this decomposition, wi(t) represents the contribution of neuron n to component r, u, captures the
temporal dynamics within each segment for component r, and Ui (r) accounts for variations across
segments for component The factor matrices = [wn (r) U= [u]"), V= [ were estimated
using an alternating least squares (ALS) optimization procedure, which iteratively updates each
factor while keeping the others fixed to minimize the reconstruction error. To analyze the interactions
between neurons, for each component r, we computed the outer product:
Z (r) (r) w (r) (r)
(r) ())
This four-way way tensor (r) encapsulates the pairwise interactions between neurons modulated by segment and temporal dynamics. Summing over the segment and time dimensions yielded matrix E(r) that represents the average interaction pattern for component r:
T
K
k.t.t
Finally, averaging over all components provided the overall interaction matrix:
E
1
(r)
R
This matrix E offers a low-dimensional representation of the neural population's correlation structure, capturing both within-segment dynamics and across-segment variations.
17

---

Algorithm 3 Sequence Non-negative Matrix Factorization (seqNMF)
Input: Data matrix X TRNXT, number of components K; sequence length I,; regularization parameter A; smoothing matrix S E RTxT with S 1 i - jl L, else 0
Initialize: W E N and @ RRXT with non-negative values
3 for iteration = to Nie do
Compute reconstruction: X k== W A.k:: H
Compute cooss-oorellatin matrix: CR. Σ E Wi,kd" Xn,+4 for all k + k'
Update W and by minimizing:
L =
Σ kig
CRR S S.H Hgal
using multiplicative gradient descent
end for
for k == to K do
10
Compute component tensor: Zkk) W A,
W sh,1 & HR.
11
Compute interaction matrix: (k)
11=1 ((
12 end for
13 Compute average interaction matrix: E
Lk=1 k
E(4
Output: Factors W, H, and interaction matrix
C.3 Sequence Non-negative Matrix Factorization (seqNMF)
To extract repeated temporal motifs from high-dimensional neural recordings, apply Neq regularized form of convolutional NMF introduced by Mackevicius et al. (34). Given a neural activity matrix X IR NT,, where is the number of neurons and T is the number of time points, SeqNMF decomposes X into K components via a convolutional model:
K
W k. H
Here, We E N encodes the sequence templates of length L, and H E RxxT specifies when each template is active. The convolution operates over the temporal axis. The model is fit by minimizing the objective:
C- C-
x
||
+ W X S. S.H 11,
where S TxT a smoothing matrix with ones on a band of width L = 2 around the diagonal, and À 0.001 controls the strength of regularization term that penalizes redundancy across components. The factors W and are updated via multiplicative gradient descent 43]. After training, the interaction structure between neurons is estimated from each sequence component. For each k, we compute correlation tensor:
Z W W. & Ha:
Then, summing over time yields:
K
(()
(b) lidt
and
E
1
(()
K
k=1
This provides a low-dimensional representation of temporally organized pairwise neural interactions.
18

---

D Comparison with the GLM
The GLM is widely used approach for modeling neural spike trains. While it has proven effective for predicting neural activity and estimating coupling filters, its formulation differs fundamentally from our proposed framework in several respects.
In the GLM, the weight matrix w, which governs pairwise influences, is directly parameterized and optimized from random initialization via oog-likelhood maximization. There is no separate function or intermediate representation devoted inferring structure. In contrast, our model introduces dedicated structure learning module that explicitly extracts latent embeddings of neuron spike trains and predicts edge strengths using MLP. This step isolates the task of estimating connectivity and makes it learnable independently of the spike generation process, enabling the model to reason about connectivity patterns beyond what improves immediate likelihood fit.
Moreover, the GLM computes the log firing rate of each neuron as a linear projection of the temporally filtered spike history of its neighbors onto the coupling filters. This process lacks latent embedding space or dynamic aggregation step. Crucially, the influence from other neurons is not passed via abstract latent messages. but directly via convolved spike histories. This approach enforces an additive, filter-based influence model with limited capacity to capture nonlinear or recurrent interactions. In contrast, our spike prediction module performs message passing on top of learned embeddings. Each neuron's 's spike history is first embedded into a latent state, messages are computed between neuron pairs based on these embeddings, and then dynamically integrated using gated recurrent unit. This framework supports richer, nonlinear temporal dependencies and is known to be capable of modeling a wide class of dynamical systems.
In summary, the key distinction lies in the explicit modularity of our model and its reliance on latent, dynamically integrated representations. which enable capture and refine structural hypotheses based on rich temporal dependencies. This sets it apart from the GLM, which lacks both the architectural separation and the representational flexibility to perform such inference.
E
Head Direction Cell Benchmark
The head direction cell benchmark dataset, available from CRCNS (Collaborative Research in Computational Neuroscience). contains extracellular recordings from the antero-dorsal thalamic nucleus and post-subiculum of freely moving mice as they foraged for scattered food in 53 x 46em open arena, specifically designed for analyzing head direction (HD) cell dynamies [36]. Although not every recorded neuron exhibits HD properties, a session was chosen based on having the highest number of neurons meeting HD score threshold, resulting in simultaneously recorded HD cells (F. The HD scores were computed by binning head direction data into 3 intervals, computing firing rates within these bins, smoothing the results with a 14.5% boxcar filter, and calculating the Rayleigh vector length from the resultant tuning curves. The threshold for identifying HD cells was defined using the 99th percentile of a null distribution generated by by shuffling spike trains. To avoid verrrepreeentttionn of any particular angle on the ring network, sessions were further screened to minimize the KL divergence between each session's preferred direction distribution and a uniform distribution. These steps produce dataset well-suited for inferring underlying circuits of HD neurons.
F Evaluation Metrics
To evaluate weight matrix where each neuron's 's weight vector (i.e., each row) is expected exhibit rotational invariance, the first step is to align all rows to a common reference phase. This is necessary because the weight vectors can be interpreted as the same underlying pattern, each shifted by different rotation. Therefore, each row must be shifted so that all vectors share a consistent phase. Additionally, the inferred weights may differ from the true weights by arbitrary global scale factor-this this can result from the multiplicative interaction between connection strength and incoming messages in 66. Since such scaling should not affect the evaluation, we first compute the average aligned weight vector w across neurons. We then learn scale factor that minimizes the L {ddisttnnc between and the true weight profile (e.g., a representative row from W), and apply this same scaling to the entire inferred matrix w to obtain W. Finally, we assess the quality of the inferred
JW-W WE where weights using the normalized inference error B1), given by A [[]] the Frobenius norm.
|| denotes
19

---

Figure 6: Tuning curves of 19 simultancously recorded HD cells that passed the HD cell threshold criteria. Each polar plot represents the directional tuning of a single HD cell, showing firing rate as function of head direction. The cells displayed here include those with the most widely distributed preferred directions across all sessions and animals, highlighting the diversity of head direction representations within the recorded population. Red dots indicate the peak firing direction (preferred head direction) for each cell.
As separate measure from connectivity inference, we evaluate spike prediction performance using log-li log-likelihood-. based metric that compares our model to homogeneous Poisson process [4] 45 I This metric quantifies how much better the model predicts spike timing relative baseline that assumes a constant firing rate. A higher score indicates that the model more accurately captures the temporal structure of neural activity. Specifically, the metric defined as
N Den
[<x] log XE 1--M X1 (x) (og d. Xx]
Dbps
(19)
N
Σ
where X; =; ET of represents the empirical average firing rate of neuron By normalizing the aigliielihood improvement by the total number of spikes, this metric yields an interpretable value in bits per spike. The final score Lbps obtained by averaging across all neurons.
G Externally Driven Rotation of Activity Bumps in a Ring Network
To incorporate external cues-such head direction, denoted by into ring attractor model, we adapt the synaptic connectivity to dynamically reflect changes in the input. In the standard ring attractor, stationary activity bumps are maintained through symmetric recurrent connections. Here, we introduce mechanism in which the synaptic weight matrix is shifted in real time based on the external signal 0(t), causing the activity bumps to rotate along the ring in synchrony with the stimulus.
2ar forming uniform eircular Each neuron E {0, I,..., NN 1} is assigned a preferred angle Φ; N arrangement. This setup allows us to interpret neuron indices as discrete angular positions, facilitating a direct mapping between the continuous angular domain and its discrete neural implementation. The synaptic weight from neuron i to neuron j. adjusted by the input angle 6(t), is given by:
Wo (d(Ø, - Φ5- - -8)), d(a)) = mod(a + n, 2T) n
20

---

where d(-) computes the signed shortest angular distance. The base connectivity profile W6(6) is modeled using a Mexican-ha kernel:
8² Wo( - exp
82
a axep
2
20
In practice, the rotation induced by (t) is implemented by translating neuron indices via integer shifts: the continuous input is scaled a discrete offset ( N-g-9(t)
where g gain factor that controls how strongly the external input steers the bump position. The weights are then indexed as (j - (0)) mod to maintain circular continuity. This shifting scheme results in coherent translation of the activity bumps that tracks the external input ((t) while preserving the internal activity profile (see Fig The network thus transitions from a self-sustained attractor to one that is steerable, integrating structured external input with internally generated dynamics.
2x
H Details of Compute Resources
All experiments in this study required approximately 7 GPU days on NVIDIA V100 32GB GPUs. Both the generative recurrent networks and the inference models were implemented using the PyTorch [66] and PyG 47 libraries. While the inference model relatively lightweight to train, generating synthetic spike trains from recurrent dynamical models can take several days depending network size and simulation length. Nonetheless, overall compute demands remain modest and accessible.
21


Reviews:
[{"review_id": "gPCMvMdGGw", "reviewer": "Reviewer_Yi57", "summary": "This paper addresses the important challenge of estimating synaptic connectivity from neural population activity. The authors propose a GNN-based approach that simultaneously predicts neural activity and infers connectivity given past observed neural activity and stimuli as inputs. The proposed model has two components: (i) a structure learning module -- an MLP that estimates pairwise connectivity from the embeddings of a pair of neurons' activities, and (ii) a spike prediction module that is a GNN, where each node corresponds to a neuron, and uses message passing to compute node embeddings. The message passing operations of the GNN, specifically the aggregation step, use the pairwise connectivity computed by the structure learning module. These node embeddings are then decoded into log firing rates. The entire network is trained by minimizing the Poisson negative log-likelihood of the observed neural activity. The proposed approach is first evaluated on simulated data from ring attractor networks. The results demonstrate that the model outperforms baselines at both predicting neural activity and inferring the connectivity across a range of conditions, including varying recurrent connectivity, external inputs, and incomplete observations. The authors also show that this framework can effectively recover key features of the underlying circuit in real neural data, which is promising.", "questions": "- Does $f_{out}$ in eq 1 include a nonlinearity?\n\n- The results reported in tables 1, 2 and 3 are averaged over multiple trials. Could you clarify what the trials refer to? Are these multiple trials from the test split? Or do they correspond to seeds for training the network? \n\n- In Fig 3C, the error bars are wider in the bottom row compared to the top. Why is that the case?", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "This paper proposes a novel approach for inferring synaptic connectivity from neural population activity, an important challenge in computational neuroscience. The authors demonstrate that their method successfully predicts neural activity and connectivity in simulated networks, outperforming existing baselines. Importantly, they also show the method's potential by applying it to real neural recordings. They also propose interesting future directions, such as applying the model to grid cell data.\n\nThe authors also discuss a few limitations of their work: (i) the model currently assumes that the network connectivity is static, and (ii) different underlying circuit configurations can produce similar neural dynamics, thus the model can correctly predict neural activity while inferring connectivity that is different from the ground truth. Another interesting question is the model's scalability. The GNN in the spike prediction module assigns one node per neuron. How efficient would this approach be for larger systems with thousands of neurons or more, particularly under conditions of partial observability where only a few hundred neurons are recorded?\n\nOverall, the paper is very clearly presented. Minor suggestions for improvement would be to add color bars to the weight matrix heatmaps and add axis labels to the few panels missing them.", "was_revised": true, "final_justification": "Thank you to the authors for addressing my questions. I am happy to maintain my initial score of 5.", "rebuttal": "We thank the reviewer for the thoughtful feedback and address each point below. All clarifications will be incorporated into the revised manuscript.\n\n> Q1. Scalability under partial observability (thousands of neurons, only a few hundred recorded)\n\n- The dominant cost in our spike prediction module is the fully connected message‑passing step of the GNN, which scales as $O(N^2)$ with the number of nodes $N$. When only  ~10% of neurons are observed, the effective graph size drops by an order of magnitude, while the number of auxiliary nodes remains a small constant. Consequently, both computation and memory fall by roughly two orders of magnitude, making the method practical even for circuits containing many thousands of neurons.\n\n> Q2. Missing color bars and axis labels\n\n- Color bars and the previously missing axis labels have been added to all weight matrix heat maps and line plots.\n\n> Q3. Non‑linearity in $f_{out}$ (Eq. 1)\n\n- Yes. After the final linear projection we apply a ReLU activation. This has now been stated explicitly in Eq. 1.\n\n> Q4. Definition of “trials” in Tables 1–3\n\n- “Trials” correspond to independent training runs with different random seeds. For every condition we train three such models and report the mean across them. The captions of Tables 1–3 have been updated accordingly.\n\n> Q5. Wider error bars in the bottom row of Fig. 3C\n\n- Weight Profile II is a broad, smoothly varying inhibitory Gaussian. Its nearly flat tail means that small changes in higher‑frequency coefficients have little impact on predicted spikes; the likelihood surface is therefore much flatter along those directions, and parameter estimates vary more between runs. In addition, the stronger and uniform inhibition of Profile II suppresses firing over much of the population, sharply reducing the number of informative spikes. With both weaker sensitivity and fewer observations, parameter uncertainties are larger, producing the wider error bars seen in the bottom row of Fig. 3C."}, {"review_id": "yG4QOHAM0m", "reviewer": "Reviewer_uHFU", "summary": "This paper introduces a self-supervised approach for inferring neural circuit connectivity from population spike activity using graph neural networks (GNNs). The model consists of two main modules: one for learning structural connectivity and another for predicting future spiking activity. With such design, the model can be applied to infer structural connectivity even under partially observed circuits. The method is evaluated on both synthetic data and real spiking recordings.", "questions": "The current work focuses on ring attractor networks. How might the proposed GNN architecture and self-supervised learning objective need to be adapted or extended to infer connectivity in neural circuits with fundamentally different topological structures or functional properties?\n\nFor the real HD neural, how the parameter $N$, $N_o$, $N_h$ are decided, especially is there any general instruction for different real neural datasets?\n\nHow the proposed method can be generalized to large scale neural activities? What is the computation cost, and inference cost,", "limitations": "Yes.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 3, "originality": 2, "strengths_and_weaknesses": "## Strength:\n\nIt is interesting that the introduced method can be applied to partially observed neural circuit. The task is challenge but is existing for the neural recordings. \n\nAuthors provide thorough discussion on different circuit scenario from fully observed network, weakly corrected fully observed, weakly correlated and partially observed.\n\n## Weakness \nConnectivity pattern with local Mexican-hat profile is used for synthetic data generation, it would be good to discuss how the performance change with different connectivity patters. \n\nRelated work discussion is not completed, for example how is the method comparable to the \n\n[1] NetFormer: https://openreview.net/pdf?id=bcTjW5kS4W\n\n[2] AMAG: https://proceedings.neurips.cc/paper_files/paper/2023/file/1c70ba3591d0694a535089e1c25888d7-Paper-Conference.pdf", "was_revised": true, "final_justification": "Authors answered my major questions, and explained the difference between proposed method and previous methods, so I raised my score to borderline accept.", "rebuttal": "Thank you for the reviewer’s valuable feedback. We respond to each point in detail below, and will ensure all clarifications are included in the revised manuscript.\n\n> Q1. Topology generalization\n\n- The reliance of our framework on ring geometry arises solely from the optional geometric feature provided to the spike prediction module. The core self‑supervised objective, predicting future spikes while co‑learning a latent connectivity matrix, remains unchanged.\n- In the current implementation, the message function $\\phi(\\cdot)$ receives only the concatenated node embeddings and therefore imposes no explicit ring metric. When geometric priors are desirable, an additional attribute such as one-dimensional distance (chains), two-dimensional Euclidean or toroidal distance (cortical sheets, grid cells), or a learned positional embedding can be used in the message function. Injecting or omitting this attribute toggles the inductive bias without altering the loss, the dual-module architecture, or the optimization routine, allowing the same inference network to accommodate a wide range of circuit topologies.\n\n> Q2. Choosing $N,N_o,N_h$ for real HD data\n\n- For the head‑direction (HD) dataset, the number of observed nodes $N_o$ is dictated by the recording and tuning‑curve selection, yielding 19 simultaneously recorded HD cells. We then choose the total network size $N$ so that the angular bin width $2\\pi/N$ is finer than the narrowest observed tuning curve; with $N=100$ this criterion is satisfied for both the synthetic and real data. To keep the synthetic and experimental settings comparable, we match the values of the parameters as closely as possible by utilizing the same number of hidden neurons in a network of the same size. The same heuristic transfers to other datasets: fix $N_o$ by the recording, pick $N$ to meet the desired spatial or functional resolution, and sweep a small set of $N_h$ values to provide sufficient latent capacity without over-parameterizing the model.\n\n> Q3. Scaling to large populations\n\n- Fully connected message passing scales as $O(N^2)$ in memory and computation, which becomes prohibitive for very large populations. Three complementary modifications might alleviate this burden. First, replacing the dense graph by a sparse adjacency built from k‑nearest neighbours or learned top‑k attention masks reduces the complexity to $O(kN)$ with $k\\ll N$. Second, grouping unrecorded neurons into coarse auxiliary sectors can further cuts the number of edges. Third, expressing the connectivity matrix in low‑rank form $W=UV^\\top$ stores two tall matrices instead of the full tensor. Together these changes have the potential to keep both training and inference costs manageable while preserving the self‑supervised objective and predictive accuracy.\n\n> Q4. Relation to NetFormer and AMAG\n\n- We thank the reviewer for pointing out two closely related studies and have expanded Section 2 accordingly.\n\n- *Paragraph to be inserted in Section 2:* Recent work that jointly models spike sequences and latent graphs splits into two main lines. NetFormer treats each neuron’s recent spike history as a token and derives a step‑wise attention matrix whose entries are interpreted as couplings, enabling it to update the graph anew every time step [1]. Because the attention is used directly to make the next prediction, NetFormer lacks the dedicated message‑passing neural network that, in our framework, forms the spike‑prediction module; consequently, it forgoes iterative propagation through a weight matrix and instead recomputes interactions from scratch at each step, which limits its ability to filter indirect correlations on benchmarks with strong recurrent connections (see Fig. 15 in the original paper). AMAG does employ message passing but begins with a randomly initialized (or correlation‑based) adjacency matrix [2]; as such, it does not attempt to infer latent circuit structure directly from spikes, but rather tunes an existing graph for better forecasting. In contrast, our approach comprises two explicit and functionally specialized modules: a structure learning block that estimates a continuous connectivity matrix from population activity, and a separate message‑passing block that uses this matrix to predict future spikes. This division enables stable connectivity estimation, accommodates unobserved neurons via auxiliary nodes, and maintains robust predictive accuracy in circuits with strong recurrence.\n\n  [1] Lu, Ziyu, et al. \"NetFormer: An interpretable model for recovering dynamical connectivity in neuronal population dynamics.\" *The Thirteenth International Conference on Learning Representations*. (2025)\n\n  [2] Li, Jingyuan, et al. \"Amag: Additive, multiplicative and adaptive graph neural network for forecasting neuron activity.\" *Advances in Neural Information Processing Systems*. (2023)"}, {"review_id": "g0x1d8LKrd", "reviewer": "Reviewer_UCAq", "summary": "This paper addresses the interesting question of how to infer synaptic connectivity from neural population activity. The authors develop a graph neural network (GNN)-based model that predicts both latent connectivity and neural activity by designing two functionally specialized network modules. As a result, they demonstrate consistently strong performance across diverse conditions on both synthetic and real datasets. Overall, although I am not an expert in GNNs, I find this paper easy to read and well supported by extensive results (i.e., various types of datasets and generative conditions).", "questions": "- In real-world recordings, neural activity is easily affected by various sources of noise. It would be helpful for the authors to further explain or validate the robustness of the proposed method under such conditions.\n\n- In Tables 1 and 2, what is the magnitude of the standard deviation for the listed models?", "limitations": "- In real-world recordings, neural activity is easily affected by various sources of noise. It would be helpful for the authors to further explain or validate the robustness of the proposed method under such conditions.\n\n- As partially discussed by the authors, one main limitation is the assumption that latent connectivity is fixed and follows a one-to-one mapping (i.e., each neuron pair has only one connection). It would be beneficial to discuss how the proposed method could be adapted to scenarios involving synaptic plasticity or more complex mapping relationships.", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 4, "strengths_and_weaknesses": "+ The writing is clear and well supported by extensive results.\n\n+ The paper addresses interesting research questions and demonstrates strong performance.\n\n- I did not identify any significant weaknesses.", "was_revised": true, "final_justification": "The authors have addressed  my concerns.  I am inclined to recommend acceptance.", "rebuttal": "We are grateful for the reviewer’s thoughtful input. Please find our detailed responses below. The suggested clarifications will be integrated into the revised version.\n\n> Q1. Robustness to biologically realistic noise\n\n- Neural recordings are affected by several well‑known noise sources: spike‑sorting errors that delete or mis‑assign individual events, and population‑wide modulations driven by unobserved inputs. Our current implementation already mitigates parts of this spectrum without additional mechanisms. First, the likelihood is Poisson and evaluated over time bins that pool many spikes; random mis‑sorts therefore act as small fluctuations whose effect on the log‑likelihood is down‑weighted relative to systematic co‑fluctuations across the population. Second, the message‑passing network aggregates information over multiple neighbours (and/or layers possibly), so an isolated spike error rarely propagates far through the graph. Third, the auxiliary nodes provide a low‑dimensional buffer for shared drive: correlations that are common to many observed neurons tend to be absorbed by observed‑to‑auxiliary edges rather than inflating observed‑to‑observed weights. While these ingredients offer a degree of robustness, we acknowledge that extreme noise could still degrade performance. Incorporating explicit noise‑handling components such as spike‑dropout augmentation is therefore a natural extension and will be explored in future work.\n\n> Q2. Magnitude of standard deviations in Tables 1 & 2\n\n- Each entry in Tables 1 and 2 is the mean across three independent training seeds; the corresponding sample standard deviations are listed below.\n\n  | **Table 1: no external input** | Δ ↓ (Thresh.)     | 𝓛_bps ↑ (Thresh.) | Δ ↓ (LNP)         | 𝓛_bps ↑ (LNP)     |\n  | ------------------------------ | ----------------- | ----------------- | ----------------- | ----------------- |\n  | **GNN**                        | **0.061 ± 0.006** | **0.882 ± 0.004** | **0.049 ± 0.005** | **0.876 ± 0.004** |\n  | GLM                            | 0.244 ± 0.012     | 0.695 ± 0.009     | 0.238 ± 0.013     | 0.712 ± 0.010     |\n  | seqNMF                         | 0.789 ± 0.021     | –                 | 0.796 ± 0.023     | –                 |\n  | TCA                            | 0.762 ± 0.019     | –                 | 0.761 ± 0.020     | –                 |\n\n  \n\n  | **Table 2: with external input** | Δ ↓ (Thresh.)     | 𝓛_bps ↑ (Thresh.) | Δ ↓ (LNP)         | 𝓛_bps ↑ (LNP)     |\n  | -------------------------------- | ----------------- | ----------------- | ----------------- | ----------------- |\n  | **GNN**                          | **0.073 ± 0.007** | **0.916 ± 0.005** | **0.058 ± 0.006** | **0.924 ± 0.005** |\n  | GLM                              | 0.259 ± 0.014     | 0.724 ± 0.010     | 0.245 ± 0.015     | 0.748 ± 0.011     |\n  | seqNMF                           | 0.791 ± 0.022     | –                 | 0.794 ± 0.023     | –                 |\n  | TCA                              | 0.760 ± 0.020     | –                 | 0.762 ± 0.021     | –                 |\n\n  \n\n  The standard deviations are roughly an order of magnitude smaller than the performance gaps between methods, so all improvements remain statistically significant. These revisions have been incorporated into the manuscript tables.\n\n> Q3. Extending to synaptic plasticity or many‑to‑many mappings\n\n- Our current experiments assume that the latent connectivity is stationary and represented by a single weight for each neuron pair, but the framework is readily adaptable to richer scenarios, a direction we are actively pursuing. For synaptic plasticity, the idea is to let the structure learning module update its connectivity map gradually as a learned function of recent spike history, so the graph drifts over time while remaining smooth on the short horizons needed for prediction. Because the spike prediction block always consumes whatever graph the structure module outputs at the current step, no architectural change is required beyond adding a lightweight plasticity head that proposes small weight updates at each iteration.\n- A second extension handles circuits in which multiple interaction types coexist. Rather than storing a single adjacency, we maintain a small set of edge‑type matrices and learn how strongly each contributes to the effective coupling at prediction time. The spike prediction module then aggregates messages across these channels, allowing it to capture multiplexed or many‑to‑many relationships. Implementing and validating these dynamic‑ and multi‑graph versions is part of our ongoing research agenda and will be outlined in the revised Discussion section."}, {"review_id": "2M7GNO5y5f", "reviewer": "Reviewer_qw7z", "summary": "The authors contribute to the question of computational neuroscience of how connectivity can be extracted from spiking data. They focus on two datasets. The first one is a simplified variant of the ring-model that was previously used by Das and Fiete (as cited). The second data set consists of electrophysiological recordings from head direction cells.\n\nThe  paper's main contribution is to address a well-known and acknowledged issue on connectivity estimation: the overestimation of connection in the presence of  strong correlation in activity (Das and Fiete 2020).\n\nThe proposed method of the authors uses a message-passing approach of graph neural networks and gated recurrent units. It works well on the artificial data (ring model of Das and Fiete).", "questions": "1. The authors show that on the ring data,  their approach indeed 'explains away the correlations that arise from co-activated neurons'. It would be great if they could explain how their equations (3) -(6) achieve this.\n\n2. The improvement over GLM in the weakly correlated case is not  impressive. Moreover, I had (may be wrongly) the impression that the authors injected inputs by changing the recurrent weights. This unusual procedure may put the GLM method at a disadvantage. How does the performance of the proposed method compare to GLM in a more classic setup when the external inputs are injected directly into neurons, for fixed weight matrix?\n\n3. As an additional control for the GLM fit: If you initialized the GLM method with the correct weight matrix (and use the procedure under point 2), would  the GLM fitting procedure (as a minimal condition) stay at the solution or still move away from the correct solution during fitting. ?\n\n4. How would the procedure work on more realistic spiking neural networks. In general, the model of a ring-attractor is of rather limited generality because of its rotation symmetry. No real data ever will have these symmetries. Can you test the method also on some more realistic spiking neural network?\n\n5. Line 129, what is the origin of the symmetry of the weight matrix? In Fig 5, is N_0=21 or 19 correct?\n\n6. How does a simulation of the final extracted model in the ABSENCE OF INPUT look like? Can you ensure that there is no run-away activity? Do bumps form spontaneously?\n\n7.  I wonder if the authors could test on any grid cell recordings, which is closer to the multiple bump setup than the chosen head direction data set.\n\n8. The choice of comparison points/alternative methods is somewhat arbitrary and limited.\nAlternative methods: The authors may want to also consider  work of Felipe Gerhard [1,2],  Guillaume Bellec [3], or Schimel [4[\n\n[1] F. Gerhard et al., (2013) Successful Reconstruction of a Physiological Circuit with Known Connectivity from Spiking Activity Alone; PLOS Comput. Biol. https://doi.org/10.1371/journal.pcbi.1003138\n\n[2] F. Gerhard et al. On the stability and dynamics of stochastic spiking neuron models: Nonlinear Hawkes process and point process GLMs, PLoS computational biology 13 e1005390 (2017)\n\n[3] G. Bellec et al. (2021)  Fitting summary statistics of neural data with a differentiable spiking network simulator. Neural Information Processing Systems 34 (NeurIPS 2021)\n\n[4] M Schimel et al., (2021) iLQR-VAE: control-based learning of input-driven dynamics with applications to neural data], bioRxiv, 2021.10. 07.463540", "limitations": "no issues", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 3, "originality": 4, "strengths_and_weaknesses": "The artificial data set was designed to contain strong correlations that not directly reflected in the weights. Based on the artificial data, the method proposed in the manuscript does seem to work. The problem it addresses is arguably serious and essential and hence an important contribution. This is clearly a strength of the paper and the title makes clear that this is also the main focus.\n\nUnfortunately, the case with hidden neurons does not seem to work very well which reduces the importance of the approach.\n\nThe results on the real experimental dataset  (head direction cells) do not offer any new insights. Moreover, I am not sure if the real dataset  is suitable to back up the main claims.\n\nThe specific artificial data set has inherent weaknesses.  In contrast to a real spiking model (e.g., nonlinear Hawkes process with strong negative self-feedback or integrate-and-fire models with reset), the neuron model used in  the artificial data has no self-feedback and is made to run in the high-noise limit (inhomogeneous Poisson process). As described, neither the  threshold-crossing model nor the LNP model have any refractoriness or reset. The authors say that they follow Das and Fiete (cited), and if this is what they do the data there is artificially made to have a mean interspike interval of 16 ms, with a broad Poisson-like interval distribution.", "was_revised": true, "final_justification": "I appreciate the good discussion with the authors. I feel that my concerns  regarding the limitations continue to be justified and I do not change my overall rating.", "rebuttal": "We appreciate the reviewer’s insightful comments and address each point below. Revisions reflecting these clarifications will be made accordingly.\n\n> Q1. Mechanism behind “explaining away” in Eqs. (3)–(6)\n\n- We thank the reviewer for requesting a clearer connection between Eqs. (3)–(6) and the model’s capacity to “explain away’’ correlations. The key mechanism is that the weighted message‑passing update enforces *dynamical consistency* between the currently estimated graph and the next‑step spike prediction, while the separation into structure learning and spike prediction modules converts every prediction error into a weight‑specific adjustment.\n- Concretely, each hidden state update in Eq. (6) aggregates messages of the form $w_{ij}\\mathbf m_{ij}^{t}$. If a weight $w_{ij}$ is spuriously large because neurons $i$ and $j$ happen to co‑fire, the term $w_{ij}\\mathbf m_{ij}^{t}$ over‑drives $\\mathbf h_i^{t+1}$, increasing the prediction loss at time $t+1$. The gradient of this loss flows through the multiplier $w_{ij}$ and therefore reduces that edge. Weights whose messages improve the prediction receive positive updates. The two‑module architecture sharpens this corrective signal: the structure head proposes the connectivity matrix, the prediction head immediately evaluates it by predicting spikes one step ahead, and the shared objective aligns both tasks. Because no auxiliary pathways exist in the fully observed setting, shared variance must be routed through accurate pairwise weights, which discourages edges that merely reflect co‑activation and retains those that convey causal influence.\n- We will incorporate this explanation at the end of Section 3.2.2 to clarify how Eqs. (3)–(6) drive the observed “explain‑away’’ behavior.\n\n> Q2. Fixed-weight external input\n\n- We thank the reviewer for raising the possibility that our weight shift implementation of head direction input might bias the comparison with the GLM baseline. We chose this formulation primarily for clarity: translating a weight kernel by an integer offset makes the bump’s motion and the connectivity profile easy to visualise. Importantly, shifting the weight matrix by $\\Delta \\theta$ is algebraically equivalent to holding the weights fixed and injecting an external input proportional to $W_0(\\phi_i-\\theta(t))$, up to an additive constant that is absorbed by each neuron’s bias term. A GLM equipped with separate bias and stimulus regressors therefore operates on the same effective synaptic drive and is not intrinsically disadvantaged.\n- To verify this point, we repeated the experiment in a conventional fixed‑weight setup. The recurrent matrix remains $W_0$, and each neuron receives an input $b\\left(1+\\boldsymbol{\\xi}(t)\\right)+\\alpha W_0(\\phi_i-\\theta(t))$, where the gain $\\alpha$ is chosen so that the bump amplitude matches the weight shift case. Under this stimulus the ranking of methods is unchanged: our model still produced lower weight recovery error and slightly higher spike prediction log‑likelihood, while the absolute gap to the GLM in the weakly correlated regime remains modest. These control results will be included in the supplementary material.\n\n> Q3. GLM initialized at the true weights\n\n- As a further control, we fit the GLM starting from the ground truth connectivity. When we initialize the GLM with the exact connectivity and keep that matrix fixed, the optimizer remains within a small neighbourhood of the correct solution ($\\Delta\\approx 0.024$) and achieves a minor log-likelihood gain, indicating that the fitting procedure is stable when the graph is \"known\". Allowing the weights to update causes the GLM to drift away from the true matrix even with strong $\\ell_2$ regularization.\n\n> Q4. Generality beyond ring attractor data\n\n- We agree that the generative recurrent networks possess strong rotational symmetry and therefore constitute only one of many possible circuit motifs. We would like to emphasize, however, that our *inference* architecture does **not** rely on that symmetry. As stated in Sec. 3.2.1, the weight matrix $\\mathbf{w}$ we infer is merely required to be symmetric and non-self-coupled; no rotation-invariance is imposed. All geometric structure resides in the *generative* benchmark, chosen because it offers ground truth connectivity and well-studied baselines. Consequently, the structure learning module can adopt or discard inductive biases (e.g., Toeplitz, low-rank, block-sparse) on a per-dataset basis without altering the optimization or loss.\n\n> Q5. Symmetric weight assumption and Fig. 5 typo\n\n- We thank the reviewer for noting these two points. The assumption $w_{ij}=w_{ji}$ was introduced to reduce parameter count and therefore sample complexity. It is not required by our optimization or loss. To confirm this, we repeated the fully observed experiments using a fully asymmetric weight matrix $\\mathbf{w}$. Weight inference error and spike prediction likelihood changed by less than 3%, indicating that the inference quality is essentially unaffected when the symmetry prior is removed.\n- The correct number of observed head-direction cells is $N_o=19$. The figure label 21 was a typo and will be corrected in the revised manuscript.\n\n> Q6. Behavior of the extracted network with the input drive removed\n\n- We evaluated the inferred connectivity by replacing the generative model’s ground truth weights with the learned matrix $\\mathbf{w}^\\ast$ and re-simulating the threshold-crossing dynamics after setting the external input to zero. Under this condition, $\\mathbf{w}^\\ast$ is effectively close to the Mexican-hat profile, and the no-input simulation behaves as a standard ring attractor: firing rates remain bounded (no run-away excitation) and multiple bumps form within tens of milliseconds and then drifts neutrally around the ring, confirming stability.\n\n> Q7. Testing on grid-cell data\n\n- We agree that medial entorhinal grid cells constitute a natural next test bed because their underlying dynamics correspond to a two-dimensional, multi-bump attractor. Extending our framework in that direction is technically straightforward. In our current implementation, the message function $\\phi(\\cdot)$ receives only the concatenated node embeddings $[\\mathbf{h}_i,\\mathbf{h}_j]$. It does not use an explicit one-dimensional ring distance. Adapting to grid cells therefore requires augmenting these learned embeddings with a two-dimensional positional code and leaving the remainder of the architecture and training objective unchanged. The main practical challenge is data availability: publicly released grid cell datasets typically contain only tens of simultaneously recorded units drawn from multiple lattice modules, which provides limited coverage of the continuous sheet and complicates quantitative validation of inferred connectivity. We are assembling a hybrid benchmark that augments these recordings with simulated hidden units so that ground truth weights are known for the combined population, and we plan to report a full evaluation in future work.\n\n> Q8. Baseline diversity concern\n\n- Thank you for the pointers. Gerhard et al. (2013, 2017) build on the same point-process framework we already employ, adding spike history terms, BIC-based filter selection, and a quasi-renewal stability test. These features refine and indirectly extend our GLM baseline rather than introducing a fundamentally different model. Bellec et al. and Schimel et al. move beyond GLMs toward latent or hidden-unit dynamics; incorporating them responsibly requires new hyper-parameter searches and biological constraints that cannot be completed within the rebuttal period.\n- In the revision, we will cite these works and outline a planned comparison with Bellec’s differentiable RSNN and a connectivity-constrained iLQR-VAE in follow-up experiments. We appreciate the reviewer’s suggestions, which will sharpen both our discussion and future roadmap."}]

Author rebuttal:

Meta-review: This paper proposes a self-supervised framework that jointly predicts spikes and infers latent synaptic connectivity with a graph neural network, using a dedicated structure-learning block coupled to a message-passing predictor. The model handles partial observability via auxiliary nodes and is evaluated on synthetic ring-attractor circuits and real head-direction recordings, showing improved spike prediction and more faithful recovery of weight profiles relative to standard baselines.

The paper is clearly written and technically sound, with a simple architecture that addresses correlation-driven artifacts in connectivity inference. Reviewers raised fair questions about scalability, generalization beyond ring geometries, the distinction from GLMs, and robustness under partial observability; the authors provided clarifications, added controls, and strengthened the discussion without introducing new risks or unresolved contradictions. With all reviews leaning to accept, the AC also agrees that the manuscript meets the bar for clarity, substance, and empirical support.
Decision: Accept (poster)