Scalable inference of functional neural connectivity at submillisecond timescales
Authors: Arina Medvedeva
Edoardo Balzani
Alex H Williams
Stephen L Keeley
Venue: NeurIPS 2025 poster ()
Keywords: Generalized linear model
functional connectivity
poisson point process
stochastic optimization
Primary area: neuroscience_and_cognitive_science

Abstract: The Poisson Generalized Linear Model (GLM) is a foundational tool for analyzing neural spike train data.
However, standard implementations rely on discretizing spike times into binned count data, limiting temporal resolution and scalability. Here, we develop stochastic optimization methods and polynomial approximations to the continuous-time analog of these models, and show them to be advantageous over their discrete-time counterparts. Further, we propose using a set of exponentially scaled Laguerre polynomials as an orthogonal temporal basis, which improves filter identification and yields closed-form integral solutions under the polynomial approximation. Applied to both synthetic and real spike-time data from rodent hippocampus, our methods demonstrate superior accuracy and scalability compared to traditional binned GLMs, enabling functional connectivity inference in large-scale neural recordings that are temporally precise on the order of synaptic dynamical timescales. We provide open-source implementations of both MC and PA estimators, optimized for GPU acceleration, to facilitate adoption in the neuroscience community.

Scalable inference of functional neural connectivity at submillisecond timescales
Arina Medvedeva
Edoardo Balzani Flatiron Institute
Flatiron Institute New York, NY
New York, NY
org
ebalzani www.miiflllttiooie aalllanniflattirooonsttttteee
Alex H Williams Flatiron Institute, New York University
Stephen L Keeley
Fordham University New York, NY
New York, NY orlliimmeffaiiiroiisst.bbu.. org
edu
Abstract
The Poisson Generalized Linear Model (GLM) is a foundational tool for analyzing neural spike train data. However, standard implementations rely on discretizing spike times into binned count data, limiting temporal resolution and scalability. Here, we develop Monte Carlo (MC) methods and polynomial approximations (PA) to the continuous-time analog of these models, and show them to be advan- tageous over their discrete-time counterparts. Further, we propose using a of exponentially scaled Laguerre polynomials as orthogonal temporal basis, which improves filter identification and yields closed-form integral solutions under the polynomial approximation. Applied to both synthetic and real spike-time data from rodent hippocampus, our methods demonstrate superior accuracy and scalability compared to traditional binned GLMs, enabling functional connectivity inference in large-scale neural recordings that are temporally precise on the order of synap- tic dynamical timescales and in agreement with known anatomical properties of hippocampal subregions. We provide open-source implementations of both MC and PA estimators, optimized for GPU acceleration, facilitate adoption in the neuroscience community
Introduction
As recording technologies in neuroscience advance, there is growing need to improve the scalability
of statistical methods for analyzing neural spiking activity. A key challenge in understanding neural computation lies in accurately estimating functional connectivity- statistical dependencies between neurons that reflect synaptic interactions. The Poisson Generalized Linear Model (GLM) is powerful tool for this purpose, capable of inferring both stimulus encoding properties and coupling between spiking units. However, the standard implementation of the GLM requires binning the timeseries data into large design matrix, X, of discrete spike counts. The time resolution of this binning is often coarse (~~ to 10 ms) TH5 compared to the timescale of synaptic dynamics, which rise and fall submillisecond timescales 6 This means conventional GLM implementations fail
to capture synaptic coupling filters on biophysically realistic scale [ [33 5 9. Moreover, as the bin size decreases. X grows size, posing significant computational and memory storage challenges.
'The Poisson point process GL. code poisson-process glm. git
available at http:://gittubb.commaaarrr/
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

We find that for most modern neural datasets, storing X in memory is infeasible, requiring users to batch X, which renders inference unstable even with state-of-the optimizers
Here, we propose methods that avoid these issues by considering the limit of infinitely small time bins, in which case the model becomes a Poisson point process (see e.g. Chapter 19 of [T]). Although point process models have been explored by the neuroscience community [1H17] most prior work either develops theoretical tools for continuous-tin
models without presenting fitting procedures
(e.g., convexity of the log-likelihood [ or error bounds 12), D, or explores related model classes therefore, 13 14.. In our setting of interest, a point process model is able to capture fine-scale spike time correlations between co-recorded neurons, which can be indicative of monosynaptic connections 6 77. Furthermore,
[5 16], , or uses numerical integration methods that do not scale to large datasets
we limit our benchmark comparison to discrete-time GLM implementations [
4
inputs to the model can be represented as a sequence of spike times instead of a large design matrix. However, to fit the point process model, we must numerically approximate an analytically intractable integral that appears in the likelihood function. We provide two approaches to deal with this integral: 1) Monte Carlo sampling- sampling-based approach (MC) and 2) second-order polynomial approximation, inspired by prior work [4] 18 19]. Both methods demonstrate improvements in accuracy over conventional approaches while maintaining computational tractability. Additionally, the polynomial approximation yields a closed-form expression for the Poisson log-likelihood that is quadratic in the GLM parameters, enabling fast and efficient computation. We also propose generalized Laguerre polynomials scaled by an exponential as a new set of basis functions for GLM inference. While these polynomials retain the desirable temporal smoothing properties of the traditionally used raised cosine basis 20 21), they offer orthogonality and closed-form integral solutions, enabling efficient filter identification.
We validate our models on both simulated and real spiking data. In simulations, we find that both MC and PA approaches scale favorably in compute time with recording length and population size, and show improved filter recovery compared to both the discrete polynomial approximate method and traditional GLMs. We then apply our method to real spiking data, where we analyze spike-time recordings from multiple rodent hippocampal regions [22] in a dataset whose size is computationally prohibitive for traditional batched GLMs. We show that recovered coupling filters align with empirical cross-correl corssooorss (CCGs) with sub-millssccond temporal precision, suggesting the model is able to accurately identify monosynaptic coupling between neurons. In addition, we able to use our model to isolate specific coupling filters that identify putative excitatory connections in the rodent hippocampus. We show that these isolated filters coincide with anatomical connectivity structure that is well-established in studies of hippocampal anatomy 23 224,, suggesting GLMs operating at this resolution provides new opportunities in the identification of neural circuitry from spike-train recordings.
2 Background
2.1 Diserete-time Poisson GLMs
Generalized linear models provide useful tool for predicting spiking activity of a single neuron y = (wi..... given recent population spiking activity or external stimuli X ( 2T), and set of model parameters The spike counts y/ conditionally Poisson distributed. yr Poisson yo w. I.). and the model log- likeliooddo is written as:
T
logp(y X, u) Σ ye log( ((x w)) (xx w)
where (ax w) the predicted firing rate time bin / and + : R - TR is a monotonically increasing, convex, and nonnegative function (e.g.. exponential or softplus). The central goal of the Poisson GLM, in the identification of ur, to find smooth time-varying statistical dependencies between either external stimuli individual neuronal spike trains and post-synaptic firing rates in neural population (Fig T)). These filters are typically estimated using a linear combination of a small number
"Wiile one can in principle represent X in sparse matrix format to alleviate computational burden, there is currently limited support for sparse matrix routines in libraries that compatible with modern GPUs.

---

of smooth basis functions and nonlinearity to assure non -og-nivte firing rates. The filters within neural populations reflect temporally delayed correlated firing, so called "functional connectivity," and are often thought of as proxy to anatomical synaptic connections, reflecting how populations of neurons influence each other through either excitatory or inhibitory dynamics.
Throughout this work, we will focus primarily on estimating functional connectivity filters using the GLM, and we will use the exponential nonlinearity, (-)) exp(:). )) as this is common choice in neuroscience and simplifies the log-likilhood objective. However, all of the methods here trivially work with an augmented X include stimuli, and with alternative nonlinearities, such as softplus, which is another common choice in the field (see Supplement S.4 for more details).
The traditional approach described above requires discretization of the time series, with bin size commonly chosen within the range from hundreds of milliseconds to one millisecond, depending on the system and stimulus (features) [
4 25 However, if the goal is to identify functional monosynaptic connections between neurons, which is a common motivation in modern GLMs, even I ms resolution is not sufficient. Electrophyssologici recordings in experimental neuroscience have shown that synaptic dynamics are often highly transient, with the rise and fall in firing occurring within 1-5 ms following presynaptic spike 226 71. This means that even bin sizes as small fail to accurately identify peak amplitude and timing (Fig B and C), which may be important for cell-specific synapse properties or distinguishing correlation firing patterns from synaptic activity.
For discrete-time GLMs, sampling at finer than 1 ms resolution demands prohibitively large memory allocations. The dimensionality of the feature space N depends on the number of neurons N in the recording and the number of basis functions J used to describe each neuron's 's activity history. For a given dataset, this results in design matrix X T For long recordings from a large number of neurons, computing and storing this design matrix with a sufficiently small bin size becomes non-trivial. As shown in Fig E simulating a dataset of 200 neurons at I ms or .l ms resolutions for 10-100 minutes would require an X matrix of -0122 bits, necessitating batched gradient calculations. In contrast, storing only spike times drastically reduces memory usage, making GLM computations far more tractable for modern high-resolution (submillisecond) datasets.
While batching the design matrix X for discrete-time Poisson GLM optimization is a sensible approach, it poses significant problems when practically fitting the model. In particular, due to the sparse firing patterns of neural activity, the variance in gradients across batches can be very large. Even when implementing a state-of-the-art stochastic variance-reduced gradient (SVRG) optimization which guarantees unbiased gradient estimates and minimal memory overhead (27]. we find that in practice the variance of our updates is too large to achieve good fits as compared to discrete GLMs using small enough datasets not require batching (Fig 3 Consequently, batched approaches are not only quite slow-requiring, for example, 5 hours on a dataset of 250 neurons with recording length 1000 seconds binned at 0.1 ms resolution but they can lead inaccurate model fits.
2.2 The Polynomial-Approxime GLM
Previous work has shown that approximating the nonlinearity in the Poisson likelihood with polynomial can be effective tool for scaling GLMs [8] []] These approaches use an orthonormal set of Chebyshev polynomials which provide a good approximation to GLM non-linnasistes over a wide range of values, and are effective even for just second order polynomial approximations Considering the exponential nonlinearity, the approximation can be written as exp(x) A
18
d2x² + axæ + ag, where A is the time bin size and (22 41, an are the optimal Chebyshev coefficients that minimize the mean squared error between the nonlinearity and quadratic approximation across the specified range [zo, _1a.::]].... Using this approximation, the GLM log-likelihood can be written as:
T
logr(y X, w) 2
Σ w x (yr- - a1 1) - a2w T x 0,
where terms that do not depend on w are dropped, and 1 is a vector of ones. Because the log- likelihood quadratic in the parameters, one can directly compute maximum posteriori (MAP) estimate using the sufficient statisties (E St. Σi= y, and E x.2;.) ). For more information on this approach, see [4]
3

---

Polsson Generalized Linear Model
fillers
III
MAE
1001
523 min
nonlinearity
other neurons
firing change
totel date size
III
2.0
Inear Mwithru
LILL
rormalized
$
observed
neuron
history window
M(A
CONERAALS
bin
Figure I: A Schematic of GLM for neuronal filter identification; B Simulation of realistic timescale post-synaptic conductance change estimated by GLM binned at 10, 5 and ms bins; C Mean absolute error (MAE) on filter accuracy from at various bin sizes; Normalized error of discrete- time sufficient statistics from continuously generated Poisson rates estimated using various bin sizes; Memory storage of spike times and X for 200 neurons at various recording lengths and bin sizes.
We find that the second-order polynomial approximation helpful in significantly reducing the computational time of the GLM, but batched sufficient statistics still carry large computational load and can be time-consuming on datasets with fine temporal resolution. Moreover, the binning of in the estimation of the linear and quadratic sufficient statistics
the design matrix introduces that accumulates with increasing number of spikes in the recording (Fig (Fi (see Supplement S.6.2 for more details).
3 The Poisson process GLM model
To improve the scalability and accuracy of these traditional GLM approaches, instead consider continuous-time Poisson Process GLM log-likelihood given by:
S( X(t) dt
logp(y X,
Σ log ((Un)
Here, time-varying Poisson rate X(t) is evaluated at time points designated by observed spike times yk of the post-synaptic neuron = (y1,..., 9Kk,, and the second term integrates the rate over the duration of the entire recording [0, 27. The firing rate at time & is then given by:
X(t: X, w)
Σ
w (bt to²)
where (x 1....... ... Es) represents the full of spikes and each spike æ. (n,, indicates that neuron n,, E l,..., N fired at time to; X't,, H) denotes the set of spikes occurring in the history window 11, t]: w,, R subset of weights associated with neuron n,; and & [0, H] + IR denotes nonlinear mapping onto J temporal basis functions. In this work, we select history window length H of 4-6 ms encompass expected neuronal dynamical effects. While can be easily augmented to include external stimuli, here we restrict our analysis to spike history, primarily focusing on the role of neural interactions and intrinsic dynamics synaptically relevant timescales.
Given that the intensity function X(t) is defined analytically, the first term in the Poisson process log-likelihood be computed exactly. However, the nonlinearity & makes the cumulative intensity function (CIF) fe X(7) dr intractable, and thus the second term of the log-likelihood requires approxi- mation. Here, we propose two methods approximate this integral: 1) Monte Carlo sampling-basee approach (MC) with unbiased estimator the CIF; and 2) a polynomial approximation (PA) that yields expression quadratic in the GLM parameters.

---

3.1 Monte-Carlo sampling for the CIF
To compute the second term in the objective function, 1 X(t)dt, we approximate the integral with Monte Carlo estimate. Instead of simple uniform sampling, we employ stratified sampling: the time support [0, 7] is divided into M equal subintervals, and sample points = (71...., 7M) drawn uniformly from each subinterval. Then,
T M Σ (Tm)) Rs
S' x(1)dt
provides an unbiased estimator of the integral that exhibits lower variance compared to uniform Monte Carlo sampling (see Chapter 8 in (28)) Thus, loss function for fixed sample of r is:
M
K
f(w, T)
T
Σ ((Tm) m=1
Σ R=1
log X(y).)
M
Where the second term be computed exactly. We can employ standard gradient-based optimization procedures on this objective selecting a different t at every iteration.
3.2 The Dolynomial-Approxinaa continuous GLM
Alternatively, can polynomial approximation method inspired by Zoltowski and Pillow 4 and Huggins et al. [8] to derive a tractable, scalable form for the log-likelihood's CIE By fitting second-order polynomial with coefficients 42 @1, dg minimize the mean squared error (MSE) against the true nonlinearity over a specified range, we reformulate the objective into a sum of integrals over linear terms (individual basis functions) and quadratic terms (basis function pairs). Depending on the choice of basis functions, these integrals may admit analytic solutions, enabling efficient evaluation of the log-liielihood. The polynomial-approximaei CIF is written as:
T J A(t)dt=
(' +
Σ 1,EX
Σ
w
Øbt w
to) dt
T Σ
2
02
Σ
T (b - t.a)
di an
Σ JU t.
Σ
wh Φlt - 1.)) dt Too
L.E
aaw Mw + a.m mnm w + Tao
Here, X, denotes the set of spikes from neuron and the linear term includes defined vector m NJ that contains N concatenated &o vectors scaled by respective total number of spikes per neuron, S, (m S sp, S24. SN where &p is vector of precomputed integrals for each of the basis function over
to. That is, 42 = -So (s(T)dr.
The quadratic term is symmetric block matrix M E RNJXNJ with N N blocks of size J x J. Each block M corresponds to a neuron pair (n, n') and accumulates the contributions from all spike pairs (t,, tx) with t, Xr and to' E Xn.. The entry at position (j.j') i') of the block is given by:
[M
Σ LEN.
Bunn, d5(T)op(T Sessty) dd,
(8)
W,E
where &,,,,, 16, 1 Le is the spike time difference. This integral nonzero only when &r,,t,, S H. i.e,, when the spike pair is within the interaction window. Therefore, if these basis function products can be expressed analytically and integrated in closed form, we only need to compute all pairwise spike time differences within the window [t. H, to] and sum the J x integral evaluations.
Given the quadratic expression of the CIF, the first term of the log-likelihood can be computed exactly when using the exponential inverse link function. The contributions from presynaptic spikes are pre-
computed neuron-specific vectors vb,,- = Σk
k= Er, Ex, Φ(y; tz). yiclding the compact
5

---

A
B
Raised Cosine
Laguerre
13
friing rate change [Hd)
RC
1.4
GL
GL
43
MSE
0.
0.
10
4.0 69
time (ms)
time (ms)
time (ms)
time (ms)
basis functions
Figure 2: Visualization of the first 5 RC and GL basis functions; B Best- rest-performing fits both bases onto filters generated from 100 RC bases; C Error on filter reconstruction for varying number of bases for both models.
form
> - Wn th be approximated as:
Th = w k where k RN concatenates all th Now, the full log-likelihood can
K log ((ne) - S/¹ X()dd
logply X, w)
w (k a m) a2 Mw
which admits closed- solution for model parameters For additional details, the full derivation of the quadratic polynomial approximation to the Poisson process eogllikelhhod and its extension to non-canonical link functions (e.g., softplus), please refer Supplement S.4 and S.6.
We define approximation range for the nonlinearity based on estimates of the postsynaptic neuron's 's firing rate. In simulations, where ground-truid binned firing rates available, the approxi- mation range is set between the 2.5th and 97.5th percentiles of these rates, mapped back through the inverse link function (i.c., log(-) when 4 exp). For real data, where firing rate distributions directly accessible, we center the range at the inverse of the mean firing rate and determine its bounds by maximizing cross-validated log-likel likelihood, following the approach of 44. In our analyses of neural recordings, we use an approximation interval spanning 3-7 Hz around the mean rate. Notably, wider intervals accommodate more variability in the estimated filter amplitudes but increase approximation error As result, polynomial approximation methods produce higher error when estimating the true underlying filters (simulated data) CCGs (real data, see Figs. 4B, D, and BB,C).
not
3.3 Generalized Laguerre polynomials as basis functions
We propose using scaled generalized Laguerre (GL) polynomials as basis functions for GLM temporal filters. Unlike raised-cosine (RC) bases, these functions are orthogonal under the weight th and thus can provide more efficient representation of filter variability with fewer basis functions (29). These polynomials have the added feature of following an approximate fammaa-unction envelope, in line with fine time-scale rises and slow decays that correspond to biophysical synaptic and neuronal dynamies (Fig AA). (A) The parameter a > -1 controls the long time-scale delay of the filter, a yielding standard Laguerre polynomials. We additionally add a coefficient
the input variable I that scales the rise-time of the bases. We set 1.5 and a throughout the manuscript based on initial model exploration, but find that varying these values does not dramatically change model performance (Fig. S2D).
These orthogonal polynomials better capture filters in fewer basis functions than the standard RC basis. We demonstrate this on simulated all-to-one coupled GLM whose filters are generated from 100 raised cosine bases. We simulate an 8-neuron population over 1000-second recording, with the postsynaptic neuron's baseline firing rate set to 3 Hz. On these data, we fit the continuous MC GLM using either the standard RC or GL sets of 2-7 bases. We find coupling filters are better matched using GL in fewer bases functions, with the best performing model being 3 GL bases. Figure 2 shows filter matches using 3 GL and 4 RC bases, and shows the mean error + standard deviation across all simulated filters. For more details on the properties the generalized Laguerre basis and comparison to RC, refer to Supplement S.5.

---

A
B
09 eprch
D8 epoch 2000
negative ogg-lkellhood
d scrntt butched
5007
MC 10 samples
101
MC rours)
normalized MSE
2000
750
2000
3000
3060
2003
3000
700
1.49
3506 2000
2000
3000
2000
1990
3000
2000
3330
200
1000
20000
SE0EC6
gradient step
gradient step
sample size
Figure 3: A First 3000 evaluations of negative log-likelihood objective on simulated data for batched and full discrete GLM with batch size 10 time bins; B Same as A for the continuous GLM with MC optimizer for sample sizes comparable to the evaluations in A; C Normalized MSE of the stochastic gradients relative the full gradient at the beginning and end of the optimization procedure for discrete and MC models, across different batch and sample sizes.
These bases also have the advantage of admitting straightforward closed-form solutions for both single and pairwise product basis function integrals. Given generalized Laguerre polynomial of degree n and parameter a, noted by _ (et), integrals I, of these bases have the form:
In - 10" 11 L. (ct) d/2 dt Ec, /"
n C,
th-a/22 dt
(10)
where On is a polynomial constant that depends on n and a. This admits exact integration via the lower incomplete gamma function 7(a, 2>) -h fa-le-fdl, with similar closed-form solutions available for pairwise basis function evaluations (see Supplement S.5 for derivations). While our polynomial approximation framework does not strictly require analytical solutions-an as numerical integration remains computationally efficient-we found that using these closed-form expressions yielded optimal performance both accuracy and speed for our implementation. For the remainder of this work, we run all simulations with 100 RC bases and fit all models with 3 to 5 GL bases.
4 Experiments
4.1 Stochastic gradient variance in discrete and continuous GLM
We first show that naive approach to implementing traditional GLMs on modern datasets-batching the design matrix X-fails to converge to the optimum due high variance of gradient estimates across batches. The discrete batched (DB) approach performs parameter updates on small subsets of data, resulting in highly inaccurate gradients. When comparing DB the full approach (on datasets small enough for the full design matrix X to fit in memory), we find that the GLM ligllikelhoodd converges poorly under gradient descent in the batched case, failing to reach the global optimum achieved by the unbatched version (Fig. A). This gradient variability is function of batch size, but even for batch sizes that push memory Timits, gradient error remains prohibitively high on large datasets (Fig 3C). We therefore look to other approaches for scaling GLMs to large datasets.
Our Monte Carlo (MC) approach also introduces stochasticity in gradient estimates as different samples approximate the CIF integral. However, this variability is substantially lower than that of the discrete batched approach, resulting much more stable inference with better log-likelihood values (Fig. B). This improved stability arises from two key differences: first, the spike term (first term in the likeliie)))ood is always computed exactly over all observed spikes rather than a subset; second, although MC sample size affects the accuracy of the CIF integral estimate (the second term), stratified sampling ensures uniform coverage of the entire recording duration. resulting improvement in gradient accuracy by computing the expected squared error between the true and stochastic gradients, normalized by the squared norm of the initial gradient and averaged over all parameter components:
Fig. 3 we quantify the
7

---

A
B
C
D
10
10a
Pr.
HC louse
10a
Bix (uun)
time (5)
10
10
MSE
10a
time
10
MSE
108
10a
10
10
10
10" # time bins
10
10
of time bins
of neurons
# of neurons
Figure Time to completion for discrete batched GLM (DB), discrete PA (PA-d), and our continuous models recording length increases in size; B, filter accuracy for the models in A: C. same as A as the number of neurons in the population increases; D, same as B for the model fits in C.
U = E V
V
(11)
where Vp is the true gradient at step and V is the corresponding stochastic gradient. Throughout inference, this error remains orders of magnitude higher in the DB model compared to the continuous sampling-base MC approach. Note that the error increases toward the end of training for both methods, as accurately estimating increasingly small gradient steps becomes more difficult as models approach convergence.
4.2 Continuous GLM scalability
We compare model performance and runtime across five approaches: DB GLM with an SVRG optimizer [27] (DB); the polynomial approximation method of Zoltowski and Pillow 4 (PA-d); our continuous- time polynomial approximation (PA-c): our sampling-buse -aaed Monte Carlo method (MC); and hybrid approach that initializes MC inference with PA-c estimates (a "warm start"), reducing optimization steps and accelerating convergence. First, we evaluate performance on simulated data from an all-to-one coupled GLM (N = 8), varying recording duration from 10 to 10 seconds, which spans the range of modern neuroscience recordings, with the bin size to 0.1 ms for discrete models. (Fig. AAB)). Next, assess scalability by simulating random, sparsely (10%) connected GLM with increasing population size (N 10 to N = 350) with a fixed recording length 100 sec (Fig. 4C). We evaluate model performance by computing the mean squared error (MSE) between the estimated and true filters.
While SVRG guarantees convergence given enough passes through the full data, in practice we find that even when its runtime exceeds that of all other models by orders of magniude, DB still underperforms, which is particularly evident at larger population sizes (Fig. 4D). The PA-d method is computationally efficient for smaller dataset sizes but eventually scales poorly in time and neuron number due to the cost of batch- aatch-computing sufficient statistics. In contrast, continuous-time methods utilize upr-lpaallel scans over the data, making them largely insensitive recording length while increase only moderately with population size (Fig. AA,C). In terms of estimation accuracy, the polynomial approximation methods (PA-d and PA-c) are less accurate, as expected, due their approximations in the log-likell ligliikl The MC and hybrid models achieve the best filter recovery, with the hybrid approach offering the best tradeoff between speed and accuracy (Fig. 4B,C). We note here also that PA-c slightly outperforms PA-d due to inaccuracies present in binned data, though both use identical nonlinearity and approximation ranges. Further discussion the discretization error and example filters from all models are provided in Supplement S.2.
4.3 Evaluation on hippocampal data
The hippocampus highly interconnected brain region essential for memory formation and retrieval. Its canonical trisynaptic circuit comprises the dentate gyrus (DG), CA3, and CA1 subregions, with distinct connectivity: the DG projects sparsely to CA3 via mossy fibers, with reciprocal connections back from CA3, and CA3 drives CAI via the Schaffer collaterals (Fig 5 Additionally, CA3 exhibits

---

Schaffer collaterals
CA1 neurors
D
Ca3
DG
Filter
CA3 nourons
Dentate gynus
Mcssy fibers
pre-snappic
500 1000
full
MSE
10--
CA3
filler
DG
10 P4-d
P/cc
Hybrid
(ms)
(ms)
Lime
post-synaplic
Figure 5: A Schematic of hippocampal anatomy; Aligmment filter estimates on subsets of data with CCGs calculated from full dataset; C Example estimated filters with overlayed CCGs selected from high firing rate neuron pairs; D Putative excitatory connections across hippocampal subregions.
dense recurrent excitatory (EE) connectivity hallmark feature supporting autoassociative memory dynamics [23]. While this anatomical framework well-established [24] 30), inferring monosy-
naptic connectivity and population-level spiking dynamics from multi-region electrophysiological recordings remains significant statistical challenge. Cross-correlograms (CCG) based methods are computationally demanding at large scale and require additional processing extract interpretable synaptic coupling patterns [31] 32 225) In contrast, GLMs offer a compact, efficient alternative that reduces parameter count while capturing temporal structure. This setting thus presents opportunity to evaluate our continuous-tin GLM models, which operate at submillisecond temporal resolution.
We use publicly available data from the Allen Institute consisting of 106 neurons (Nc.1 = 62,
Nc = 28, NG 16) recorded with single probe over approximately 2.7 hours 222 All
models are run with ridge regularization (3 1000), common choice for GLMs 3 9 to encourage sparsity in synaptic connections (see Supplement S. for more hyperparameter details). To assess filter accuracy, we compute the MSE between CCGs calculated the full dataset and filter estimates from hybrid PA-MC (H), PA-c, and PA-d models on various subsets of the data. We exclude the discrete batched model (DB) from this analysis, running it convergence the full dataset would be computationally infeasible. We find that our filters empirically match the pairwise CCGs, with the hybrid model showing the closest alignment even with only 500 seconds (8.3 minutes) of data, a small fraction of the full 2.7-hour recording (Fig B, C). While CCGs serve proxy for putative connections and cannot fully isolate synaptic effects from common input indirect pathways, they provide a useful benchmark for evaluating filter estimates. Furthermore, after pre-selecting filters with peaks between 0.3-2.5 ms-indicative of excitatory connections-ww find a connectivity structure that closely reflects known hippocampal anatomy (Fig. DD. Table The CA3 network exhibits the highest density of recurrent excitatory connections (~ 4%), consistent with anatomical estimates [24]. while also showing bidirectional communication with the dentate gyrus (30) and Schaffer collateral projections CA1 (Table I Notably, cross-region couplings tend to exhibit longer temporal delays (measured as time from filter onset to peak) than intra-regions latencies, consistent with axonal conduction times between structures and suggesting physiological validation of our identified connections. Fit results on the full dataset (N 623 neurons across all probes) and a comparison showing improved performance with Generalized Laguerre versus raised cosine basis functions are provided in Supplement S.1.3.
5 Conclusion
We developed continuous-time GLM implementation capable identifying fine-timescale coupling filters in modern large-scale neural recordings. rendering modern datasets (hundreds of neurons recorded for thousands seconds) trainable in minutes with sub-millisecond precision. Our focus has been detecting potential synaptic connections through coupling filters, complementing existing approaches [1] with a key advantage of being able rapidly screen candidate connections in large datasets.

---

Table 1: Putative excitatory connections and synaptic latencies across hippocampal regions.
Block
Pairs Total
Putative E
Fraction (%)
Mean Delay (ms)
CA3-CCA
784
30
3.83
1.75 0.52 0.85 + 0.31
DG-+DG
256
9
3.52
CA3-DDG
448
12
2.68
1.57 0.83
CAI
3844
51
1.33
1.69 = 0.75
DG~AC33
448
5
1.12
2.09 + 0.33
1736
18
1.04
2.15 0.34
CAI-DG
992
0.4
1.86 0.29
CAI-CAA3
1736
3
0.17
2.17 0.06
Our work complements existing continuous-time modeling efforts which have different modeling goals operate in smaller data regimes. In particular, Hawkes processes [16] represent compu- tationally efficient approach to identifying excitatory neuronal connections, but they cannot model inhibitory connections and thus occupy different model class than the general Poisson process GLM. Other models, such as continuous Point-process latent variable models [3]] share a similar likelihood construction but focus on identifying latent structure rather than fine-scale functional connectivity. To our knowledge, the only prior work that actually fits a continuous-time GLM [17] uses Gauss-Lobatto quadrature to approximate the integral in the log-likelihood. However, this approach requires inserting quadrature nodes between every spike time, making it computationally infeasible for the dataset sizes explored here (see Supplement S.3 for details). These fundamental limitations moodel structure mismatch and computational infeasibility precluded direct comparison to these methods in our benchmarks.
Our approach inherits several limitations from the broader class of Poisson GLMs, including the challenge of dissociating monosynaptic connections from correlated firing 2 and the difficulty of identifying true connectivity without overly penalizing weak dependencies or connections involving low-firing neurons 225 32). The Poisson distribution itself may be suboptimal for describing neural spiking due to its variance assumptions; flexible alternatives such as the negative binomial distribution could better capture spiking characteristics. Additionally, there fundamental trade-off' between our two approximation methods: the PA approach enables faster inference through closed- form solutions but inherently less accurate due its global approximation firing rates, while the MC approach is more accurate but requires multiple iterations to converge. Our hybrid model, which uses PA- based initialization followed by MC finetuning, our attempt to balance this trade-off.
19 34
Key future directions include: more thorough evaluation of sparsity priors for population recordings, of additional non-linearities, per-neuron approximation range optimization for our polynomial- approximate approach, and exploring variance reduction techniques [7] 35) for Monte Carlo sampling of the CIF. Additionally, extending the framework to incorporate latent population dynamics-foo instance, by modeling shared low-dimensional trajectories at slower timescales similar to GPFA could help disentangle fast coupling dynamics from slower coordinated population activity. potentially improving both interpretability and generalization to held-out neurons.
33
6 Acknowledgments
This work was supported by the Simons Foundation. AHW was supported by the NIH BRAIN initiative (IRFIMH113377)),
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 complete neuronal population. Nature, 454(7207)9995-999, 2008
[2] Abhranil Das and Ila R Fiete. Systematic emrrr in connectivity inferred from activity in strongly recurrent networks. Nature Neuroscience, 12(10))1286113 1296, 2020.
10

---

[3] Il Memming Park. Miriam LR Meister, Alexander C Huk, and Jonathan W Pillow, Encoding and decoding in parietal cortex during sensorimotor decision-making. Nature neuroxcience, 17(10):1355-1403 2014.
[4] David Zoltowski and Jonathan W Pillow. Scaling the poisson glm to massive neural datasets through polynomial approximations. Advances in neural information processing systems, 31, 2018.
[5] Jacob L Yates, Il Memming Park, Leor N Katz, Jonathan W Pillow, and Alexander C Huk. Functional dissection of signal and noise in mt and lip during decision decision-making, Nature neuroscience, 20(9):1285 1292, 2017.
[6] Daniel Fine English, Sam McKenzie, Talfan Evans, Kanghwan Kim, Euisik Yoon, and György Buzsáki. Pyramidal cell-interneuron circuit architecture and dynamics in hippocampal networks. Neuron, 96(2): 505-520, 2017.
[7] Ian H. Stevenson. Circumstantial evidence and explanatory models for synapses in large-scale spike recordings. Neurons, Behavior, Data analysis, and Theory, 2023.
[8] BL Sabatini and WG Regehr. Timing of synaptic transmission. Annual review of physiology. 61(1):
521-542. 1999.
[9] Eric Hart and Alexander Huk. Recurrent circuit dynamics underlic persistent activity in the macaque frontoparietal network. elife, 9:e52460, 2020.
[10] Robert Kass, Uri Eden, Emery N Brown, et al. Analysis of neural data, volume 491. Springer, 2014.
[11] Liam Paninski. Maximum likelihond estimation of cascade point-process neural encoding models. Nenwork: Computation in Neural Systems, 15(4))243, 2004.
[12] Don H Johnson. Point process models of single- neuron discharges. Journal of computational neuroscience, 3:275-294 2995,299... 1996.
[13] Wilson Truccolo, Uri T Eden, Matthew R Fellows, John P Donoghue, and Emery N Brown. A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects. Journal of neurophysiology, 93(2): 1074-109,, 2005.
[14] Zhe Chen, David F. Putrino, Demba F.. Ba, Soumya Ghosh, Riccardo Barbieri, and Emery N. Brown. A regularized point process generalized lincar model for assessing the functional connectivity in the cat motor cortex. In Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pages 5006-5009, 2009.
[15] Alex Williams, Anthony Degleris, Yixin Wang, and Scott Linderman. Point process models for sequence detection in high--iiensional neural spike trains. Advances neural information processing systems, 33: 14350-14366 14361, 2020.
[16] Scott W. Linderman and Ryan P. Adams. Scalable bayesian inference for excitatory point process networks, 2015. URL URL htttttts/////iiii
hwagg//bs/Siee 03228
[17] Gonzalo Mena and Liam Paninski. On quadrature methods for refractory point process likelihoods. Neural computation, 26(12):2790- 2797, 2014.
[18] Jonathan Huggins, Ryan P Adams, and Tamara Broderick. Pass-glm: polynomial approximate sufficient statistics for scalable bayesian glm inference. Advances in Neural Information Processing Systems, 30. 2017.
[19] Stephen Keeley, David Zoltowski, Yiyi Yu, Spencer Smith, and Jonathan Pillow. gaussian process factor models for spike count data using polynomial approximations. In International conference on machine learning, pages 5177-5186. PMLR, 2020.
Efficient non-conjugate
[20] Jonathan W Pillow, Liam Paninski, Valerie J Uzzell, Eero P Simoncelli, and EJ Chichilnisky. Prediction and decoding of retinal ganglion cell responses with a probabilistic spiking model. Journal of Neuroscience, 25(47): I 1013- 2005.
[21] John G Proakis and Masoud Salehi. Digital communications. McGraw-hill, 2008
[22] Allen Institute for Brain Science. Visual coding - neuropixels, 2023. URL. hattps://pont bww.b-ainnn m-ma org/explore http///xplorre//rriit//////////t wtwtttteeviaaaaoonnnneeeeeeeeeee merroinnnnuropixxel... Dataset includes spike times, L.P,, and be- havior from mouse visual cortex during stimuli presentation.
[23] Alessandro Treves and Edmund T Rolls. Computational analysis of the role of the hippocampus in memory. Hippocampus, 4(3):374-33 391, 1994.

---

[24] David G Amaral and Menno P Witter. The three-dimensional organization of the hippocampal formation: A review of anatomical data. Neuroscience, 31(3):511-591, 1989. doi: 10.1016/0266-55222899992222
[25] Naixin Ren, Shinya Ito, Hadi Hafizi, John M Beggs, and Ian H Stevenson. Model-hased detection of putative synaptic connections from spike recordings with latency and type constraints. Journal of neurophysiolooy. 124(6):1588-1666 2020.
[26] Stephamie Seeman, Luke Campagnola, Pasha A Davoudian, Alex Hoggarth, Travis A Hage, Alice Bosma-Moody, Christopher A Baker, Jung Hoon Lee, Stefan Mihalas, Corinne Teeter, Andrew Ko, Jeftrey G Ojemann, Ryder P Gwinn, Daniel L Silbergeld. Charles Cobbs, John Phillips, Ed Lein, Gabe Murphy, Christof Koch, Hongkui Zeng, and Tim Jarsky. Sparse recurrent excitatory connectivity in the microcireuit of the adult mouse and human cortex. eLife, 2018.
[27] Robert M Gower, Mark Schmidt, Francis Bach, and Peter Richtárik. Variance-reduced methods for machine learning Proceedings of the IEEE, 108(11): 1983, 2020.
[28] Art B. Owen. Monte Carlo theory, methods and examples. https://artowen. su. www.idsmannnn.... /mc/ 2013.
[29] Roelof Koekoek and Hendrik Gerrit Meijer. A generalization of laguerre polynomials. SIAM journal on mathematical analysis, 24(3):768 -7222266-782 1993.
[30] Helen E Scharfman. The ca3 "backprojection" to the dentate gyrus. Progrexs in bruin research, 163: 627-637. 2007.
[31] Zach Saccomano, Sam Mckenzie, Horacio Rotstein, and Asohan Amarasingham. A causal inference approach of monosynapses from spike trains. arXiv preprint arXiv:2405.02780 2024.
[32] Ryota Kobayashi. Shuhei Kurita, Anno Kurth, Katsunori Kitano, K. Mizuseki, Markus Diesmann, B. J. Richmond, and S. Shinomoto. Reconstructing neuronal circuitry from parallel spike trains. Nature Communications, 2019.
[33] Lea Duncker and Maneesh Sahani. Temporal aligmment and latent gaussian process factor inference in population spike trains. Advances in neural information processing systems, 31, 2018.
[34] Jonathan Pillow and James Scott. Fully bayesian inference for neural models with negative-binomial spiking. Advances in neural information processing systems, 25, 2012.
[35] Aaron Defazio, Xingyu Alice Yang. Harsh Mehta. Konstantin Mishchenko, Ahmed Khaled, and Ashok Cutkosky. The road less scheduled. arXiv preprint arXiv:2405. 15682, 2024.
12

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: We elaborate on all claims made in the abstract and introduction in the main paper, ensuring they accurately reflect our contributions.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
We discuss limitations of our work in the main text and in the supplement.
Guidelines:
The answer NA means that the paper has no 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 a separate "Limitations" section 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 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 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 is low 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 fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
13

---

Justification: Our paper does not include formal theoretical results such as theorems or lemmas, and therefore we do not require formal set of assumptions or proofs. While we provide derivations and model formulations, these are not framed as theoretical claims requiring proof.
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 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 provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results 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)?
Question:
Answer: [Yes]
Justification: We have provided the most important information about all parameters, hyperparameters and other modeling details in the main paper with additional information available the supplement.
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 be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how reproduce that algorithm. (b) If the contribution is primarily new 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 (c.g., to registered users), but it should be possible for other researchers to have some path to reproducing verifying the results.
14

---

reproducible, they 5. Open access to data and code
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?
Question:
Answer: [Yes]
Justification: We release code for simulation and model fitting, and we use publicly available data for our results, cited appropriately.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.co
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
/nips. wwtpcc/pccppcpcceeeooooooop /Codeeeeee ieaibbmsssionPo swiniisssiopiiiyyy Policy for more details.
The authors should provide instructions on data access and preparation, including how 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 to 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: Yes, we provide all required details in the main paper text and the supple- ment.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We have included error bars for all results where they are appropriate and meaningful.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
15

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run 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 the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables figures symmetric error bars that would yield results that are out range (e.g. negative error rates).
error bars are reported in tables or plots, The authors should explain 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: We specify the type of computer resorces in the main paper and provide detailes resorce profile in the appendix.
Guidelines:
The answer NA means that the paper does include experiments.
The paper should indicate type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (c.g., preliminary or failed experiments that didn't 't make 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 setiisssttptt...
Answer: [Yes]
Justification:
This paper adheres to the NeurlPs Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
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 (c.g.. if there special consid- eration due to laws or regulations in their jurisdiction..
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
We anticipate any direct societal impact of this work.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
16

---

If the authors answer NA or No, they should explain why their work has no societal impact or 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), faimess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, 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 intended and functioning correctly, 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 attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Question:
Does the paper describe 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)?
Answer: [NA]
Justification: The paper does not involve the release of data or models with a high risk of misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
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 publicly available software libraries and properly cite the benchmark model we compare against. All assets used are under standard open-source licenses, and we respect their terms of use.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
17

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, paperssi thood cww datasets 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: [Yes]
Justification: We provide open-source implementations of the models developed in this work. The code will be shared anonymously at submission and publicly upon acceptance.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moode part of 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:
The paper does not involve research with human subjects nor crowdsourcing.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
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:
The paper does not involve research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
18

---

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]
LLMs do not play a central or novel role the method developed in this research.
Justification:
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy (https://neurl)) aattop://eearrrpccc////oofeeeeee hctp:///eeeeeeeece////////.... for what should or should not be described.
19


Reviews:
[{"review_id": "Ky3BT2eT7b", "reviewer": "Reviewer_D8KZ", "summary": "The authors address an important issue in the analysis of spiking neural data: how to reasonably apply Poisson GLM without discretizing into binned count data in order to estimate functional connectivity at the temporal scale of synaptic events. Current binning approaches degrade temporal resolution and for modern data sets leads to computational limitations. The authors continuous formulation uses two tractable surrogates for the intractable cumulative intensity integral: stratified Monte-Carlo quadrature (MC) and a quadratic polynomial approximation (PA-c) that gives a closed form MAP update. They further introduce an exponentially scaled Laguerre basis as an orthogonal alternative to the raised cosine basis.They validate their methods with one Neuropixels data set, showing that MC/PA-c run faster than batched discrete GLMS, recover filters more accurately than discrete or PA-d methods, and provide better estimates of hippocampal connectivity. They provide GPU-accelerated code.", "questions": "How novel is the Laguerre basis? Is this not used in earlier work?\n\nWhy not use public large-scale data sets such as Steinmetz et al. visual cortex data and baseline against previous continuous Hawkes process models and quadrature methods from Paninksi, Lindermann, and others?\n\nDoes MC becomes prohibitively slow for high firing rates or when PA-c becomes memory-bound (long recordings, hundreds or thousands of neurons)?", "limitations": "The limitations could be better highlighted in the main paper.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n\nThis problem is of increasing importance as neural recording methods improve in speed and scale. \nThe quadratic PA-c derivation is elegant\nThe Laguerre function drop preserves biophysical realism\nThorough wall-time and memory comparisons\nGPU-optimized code provided\n\n**Weaknesses**\n\nThe paper would benefit from more comprehensive baselining. For example, it currently evaluates only against a binned Poisson-GLM variant (PA-d/SVRG) and one small Neuropixels data set. \nHow do things scale as the number of neurons increases? \n\n**Minor**\n\nNotation introduction could be made clearer: for example, not introducing what x_s is at line 140 well after it’s used in Eq 4. This is not the only case of late introduction of variables.", "was_revised": true, "final_justification": "We are satisfied with the responses provided by the reviewer, provided they include the full analysis (N=642) they suggest. We therefore raise our score to a 5.", "rebuttal": "We thank the reviewer for their detailed and insightful feedback. Below are our responses to the reviewer’s comments and questions.\n\n### **Global Response**\n\n**Comparison to existing literature/baseline evaluation**\n\nAll reviewers note that there is prior work on continuous GLMs and note a lack of comparison to these methods. However, much of the cited work are not actual model fits of a continuous GLM, and thus are not appropriate baselines for direct comparison. Some citations develop theoretical tools relevant to continuous Poisson point processes—such as convexity of the log-likelihood [1], or error bounds and information capacity [14]—and others use continuous-time models in latent variable settings [16]. Some of the other work mentioned, particularly [15] and [Chen2019] referenced by 8hfS are actually discrete-time implementations of GLMs, and so are captured in our discrete benchmark in figures 3, 4, and 5. \n\nTo our knowledge, the only prior work that actually fits a continuous time Poisson GLM is [17] which uses Gauss-Lobatto quadrature to approximate the CIF. Unfortunately, the link provided in their manuscript to public code is no longer available. We attempted to implement their method but found that a naive implementation was both computationally infeasible on datasets of the size we consider and less accurate in approximating the CIF and fitting the model on very small subsets of data. Their approach relies on inserting a varying number of quadrature nodes between every spike time, which, for our hippocampal dataset with ~5 million spikes, would require storing and evaluating far more nodes than spikes—leading to significant memory issues and prohibitively slow inference (MC, by contrast, uses many times fewer samples than spikes). We have since tested their quadrature approach on a dataset smaller than our current evaluations in the manuscript and found that the CIF estimate as well as estimated parameter values were less accurate than both our MC and PA approaches even when we use a large number of nodes. We believe this is due to the high-frequency content of the filters, which cannot be captured well by standard quadrature schemes.\n\nNevertheless, to address this concern directly in the manuscript, we will expand our related work section to clarify that nearly all GLM implementations are in discrete time. Regarding [17], although a full reproduction is not feasible, we will include a brief comparison in the supplementary material based on our partial reimplementation on a small simulated dataset, highlighting both the computational challenges and performance limitations we observed. This will help clarify the practical barriers to applying existing quadrature-based methods. We will also note that developing more scalable quadrature methods—those better suited for GPU acceleration or for estimating high-frequency coupling filters—remains a direction for future work.\n\n**Limitations of our approach**\n\nWe agree with the reviewers that the limitations of our approach were not sufficiently discussed, and we have expanded the conclusion section of the manuscript to address this. Many limitations of our model are inherited from the broader class of Poisson GLMs. This includes the challenge of dissociating monosynaptic connections from correlated firing (e.g. [6]) and the difficulty of identifying true connectivity without overly penalizing weak dependencies and connections with low spiking neurons ([24], [31], see response to reviewer 8hfS). Addressing these limitations is a priority for future work. Another limitation is the choice of the Poisson distribution itself, which may be a suboptimal match for describing spiking activity due to its assumptions about variance. Flexible alternatives such as the negative binomial distribution can potentially better capture the characteristics of neural spiking (e.g. [19] and [Pillow&Scott2012]). \n\nIn response to z5Zq’s and 8hfS’s comments about individual contributions and shortcomings of our MC and PA approaches, we note that there is a clear trade-off between speed and performance. As 8hfS noted, the PA method is inherently less accurate due to its global approximation of firing rates. However, it enables much faster inference (especially with the closed-form solution). The sampling-based (MC) approach, while more accurate overall, requires multiple iterations to converge to a solution. The hybrid model, which uses PA-based initialization followed by MC finetuning, is our attempt to enjoy the best of both approaches.\n\n**Notation**\n\nWe thank the reviewers for drawing our attention to the notation errors. We have made the necessary revisions to ensure that our use of notation and abbreviations is introduced clearly and consistently in the correct order. This includes consistent use of $\\textbf{w}$ and $x_t$ in Eq. 1, 2 and throughout the text; a more thorough introduction of the notation used in Equation 4; and consistent use of abbreviations (e.g., MC, GLM, DB). We have also aligned the notation used in the PA-c description with that of the MC method.\n\n### **Reviewer Specific concerns**\n\n>How novel is the Laguerre basis? Is this not used in earlier work?\n\nLaguerre polynomials are typically found in the physics literature (e.g. [Quesne2009]) but have been used as bases for temporal fMRI signals ([Solo2004]). To the best of our knowledge, our paper is the first one to adopt them as basis functions for constructing neuronal filters. \n\n>baseline against previous continuous Hawkes process models and quadrature methods from Paninksi, Lindermann, and others\n\nSee the general response above. To the best of our knowledge, our work presents the first implementation of a fully continuous-time GLM scalable to the dataset sizes analyzed in this paper, and we now include a comparison to quadrature on a small dataset in the supplement.\n\n We also thank the reviewer for pointing out Hawkes processes—a relevant line of work that we had overlooked. We have expanded the review of related work in the manuscript to include a discussion of the distinctions between Hawkes models and our Poisson point process GLM framework. Importantly, Hawkes processes only allow for excitatory interactions between neurons and cannot account for inhibitory interactions. On the other hand, the nonnegativity assumption built into a Hawkes process facilitates certain tricks (in particular, leveraging the superposition property of Poisson point processes). At the end of the day, both models have their virtues. Our work focuses on finding better optimization algorithms for one of them (the point process GLM), and it doesn’t make sense to benchmark optimization methods across the two since the models have fundamentally different structure. [Linderman&Adams].\n\n>Does MC become prohibitively slow for high firing rates or when PA-c becomes memory-bound (long recordings, hundreds or thousands of neurons)?\n\nOur continuous-time GLM framework works with a compact data representation (marked spike-time sequences) whose storage requirements do not scale with the number of neurons or the recording duration directly. Instead, the bottleneck is the total number of spikes in the recording which becomes limiting at the scale of tens of millions of spikes. \n\nThe reviewer is correct that the MC and PA-c methods impose distinct computational constraints. MC-based optimization requires iterating over all spike times and MC samples at each gradient step. However, the number of MC samples is independent of the number of neurons (unlike quadrature-based methods discussed above) and scales only with the recording length. In our analyses of large-scale datasets, we find that the number of MC samples remains well below the number of spikes. In contrast, PA-c approach precomputes sufficient statistics in a single pass of the data. However, this step can lead to memory overhead due to the need to store large intermediate arrays, particularly the matrix of pairwise interactions M, which grows quadratically with the number of neurons. This limits parallelization and may slow down inference in larger (>1000) neural populations.\n\n>The paper […] currently evaluates only one small Neuropixels data set. How do things scale as the number of neurons increases?\n\nIn the main text, we restricted our analysis of hippocampal data to a single probe to eliminate probe-specific artifacts that led to inferring artificial filters and confound connectivity estimation. However, we also ran preliminary experiments on the full recording (N=642 neurons) using both the MC and PA-c models. On this scale, the inference was slower (under an hour for MC and under 30 minutes for PA-c) but models fit without memory issues. If the reviewer believes these results would strengthen the submission, we can include them in the supplementary material along with a discussion of their limitations due to probe artifacts.\n\n\n**Additional References**\n\n[Chen2009] Chen, Z., Putrino, D. F., Ba, D. E., Ghosh, S., Barbieri, R., & Brown, E. N. (2009). A regularized point process generalized linear model for assessing the functional connectivity in the cat motor cortex. IEEE\n\n[Linderman&Adams] Linderman, S. & Adams, R. (2015). Scalable Bayesian Inference for Excitatory Point Process Networks. arXiv:1507.03228.\n\n[Owen2013] Owen, A. B. (2013). Monte Carlo theory, methods and examples. \n\n[Pillow&Scott2012] Pillow, J., & Scott, J. (2012). Fully Bayesian inference for neural models with negative-binomial spiking. Advances in neural information processing systems, 25.\n\n[Solo2004] V. Solo, C. J. Long, E. N. Brown, E. Aminoff, M. Bar and S. Saha. (2004). fMRI signal modeling using Laguerre polynomials. IEEE\n\n[Quesne2009] Quesne, C. (2009). Solvable Rational Potentials and Exceptional Orthogonal Polynomials in Supersymmetric Quantum Mechanics. SIGMA 5, 084."}, {"review_id": "Nu4SlvckU0", "reviewer": "Reviewer_z5Zq", "summary": "This paper used continuous time Poisson process GLM to infer functional connectivity and proposed two methods to approximate the cumulative intensity function in the log likelihood equation: Monte Carlo stratified sampling, and polynomial approximation. They also used Laguerre polynomials as basis functions for GLM temporal filters. Empirical results show that these methods achieve better convergence with less runtime than discrete batched GLM. The methods can uncover some features of hippocampus connectivity after pre-selecting temporal filters within a certain timescale.", "questions": "1. What is the main difference between this paper and previous continuous time Poisson process GLM? \n\n2. Do methods in this paper outperform those in [1] and [15] when used to infer functional connectivity? \n\n3. If so, what is the individual contribution of Monte Carlo stratified sampling, polynomial approximation, choice of temporal filters and choice of nonlinearities?\n\n4. The authors pre-selected filters within a certain time scale (line 297) when using the model on a hippocampal dataset. In general, how much prior is required when using this method on a real dataset?", "limitations": "See above.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 1, "strengths_and_weaknesses": "The framework of this paper is clear, and the methods proposed in this paper achieve better performance than the discrete time batched GLM. However, continuous time Poisson process GLM is not a new concept and has been studied extensively in the field, as is pointed out in the introduction section with citations [1][14][15][16][17], and this paper has no direct comparison to these previous models or methods.\nAdditionally, there are several places that should be written more carefully:\nEquation (1) and (2) summarize methods for discrete time Poisson GLM, but the notations are not treated carefully: $X_t$ is never introduced, and I think the authors meant to write $x_t$ in (1), but then $W^T X_t^T$ will be confusing in (2). These equations directly correspond to equations (1) and (7) of the cited paper [8], which have very clear notations, so maybe the authors should just follow that.\nFigure 3A and B has no y-axis values, making it hard to assess whether batched discrete time GLM ‘converges poorly’ as stated in line 233. Line 237 references Figure 3C instead of Figure 3B as stated. Line 245 references Figure 3B instead of Figure 3A as stated.\nThe hybrid PA-MC model introduced in line 254 should be explained in more detail, especially given that it is the best performing method.", "was_revised": true, "final_justification": "The authors have addressed my my concerns during the rebuttal and I have increase my score as a result.", "rebuttal": "We thank the reviewer for their comprehensive review and detailed feedback. In response to the reviewer’s main criticism regarding insufficient comparison to prior work, we provide the following general reply that addresses this issue across the reviewer’s comments.\n\n### **Global Response**\n\n**Comparison to existing literature/baseline evaluation**\nAll reviewers note that there is prior work on continuous GLMs and note a lack of comparison to these methods. However, much of the cited work are not actual model fits of a continuous GLM, and thus are not appropriate baselines for direct comparison. Some citations develop theoretical tools relevant to continuous Poisson point processes—such as convexity of the log-likelihood [1], or error bounds and information capacity [14]—and others use continuous-time models in latent variable settings [16]. Some of the other work mentioned, particularly [15] and [Chen2019] referenced by 8hfS are actually discrete-time implementations of GLMs, and so are captured in our discrete benchmark in figures 3, 4, and 5. \n\nTo our knowledge, the only prior work that actually fits a continuous time Poisson GLM is [17] which uses Gauss-Lobatto quadrature to approximate the CIF. Unfortunately, the link provided in their manuscript to public code is no longer available. We attempted to implement their method but found that a naive implementation was both computationally infeasible on datasets of the size we consider and less accurate in approximating the CIF and fitting the model on very small subsets of data. Their approach relies on inserting a varying number of quadrature nodes between every spike time, which, for our hippocampal dataset with ~5 million spikes, would require storing and evaluating far more nodes than spikes—leading to significant memory issues and prohibitively slow inference (MC, by contrast, uses many times fewer samples than spikes). We have since tested their quadrature approach on dataset smaller than our current evaluations in the manuscript and found that the CIF estimate as well as estimated parameter values were less accurate than both our MC and PA approaches even when we use a large number of nodes. We believe this is due to the high-frequency content of the filters, which cannot be captured well by standard quadrature schemes.\n\nNevertheless, to address this concern directly in the manuscript, we will expand our related work section to clarify that nearly all GLM implementations are in discrete time. Regarding [17], although a full reproduction is not feasible, we will include a brief comparison in the supplementary material based on our partial reimplementation on a small simulated dataset, highlighting both the computational challenges and performance limitations we observed. This will help clarify the practical barriers to applying existing quadrature-based methods. We will also note that developing more scalable quadrature methods— those better suited for GPU acceleration or for estimating high-frequency coupling filters—remains a direction for future work.\n\n**Limitations of our approach**\n\nWe agree with the reviewers that the limitations of our approach were not sufficiently discussed, and we have expanded the conclusion section of the manuscript to address this. Many limitations of our model are inherited from the broader class of Poisson GLMs. This includes the challenge of dissociating monosynaptic connections from correlated firing (e.g. [6]) and the difficulty of identifying true connectivity without overly penalizing weak dependencies and connections with low spiking neurons ([24], [31], see response to reviewer 8hfS). Addressing these limitations is a priority for future work. Another limitation is the choice of the Poisson distribution itself, which may be a suboptimal match for describing spiking activity due to its assumptions about variance. Flexible alternatives such as the negative binomial distribution can potentially better capture the characteristics of neural spiking (e.g. [19] and [Pillow&Scott2012]). \n\nIn response to z5Zq’s and 8hfS’s comments about individual contributions and shortcomings of our MC and PA approaches, we note that there is a clear trade-off between speed and performance. As 8hfS noted, the PA method is inherently less accurate due to its global approximation of firing rates. However, it enables much faster inference (especially with the closed-form solution). The sampling-based (MC) approach, while more accurate overall, can potentially suffer from slow inference since the number of required samples M grows proportional to the recording length. This limitation is partially overcome in the hybrid model through PA-based initialization.\n\n**Notation**\n\nWe thank the reviewers for drawing our attention to the notation errors. We have made the necessary revisions to ensure that our use of notation and abbreviations is introduced clearly and consistently in the correct order. This includes consistent use of $\\textbf{w}$ and $x_t$ in Eq. 1, 2 and throughout the text; a more thorough introduction of the notation used in Equation 4; and consistent use of abbreviations (e.g., MC, GLM, DB). We have also aligned the notation used in the PA-c description with that of the MC method.\n\n### **Reviewer specific concerns**\n\n> Figure 3A and B has no y-axis values,  Line 237 and Line 245 references \n\nWe apologize for the confusion with Fig.3. We have corrected the references to panels in B (line 245) and C (line 237) as pointed out by the reviewer. We additionally have adjusted the plots to include the log-likelihood values. However, want to emphasize the the purpose of this plot is to demonstrate that the discrete GLM batching approach yields highly variable log-likelihoods that remain far from the optimal log-likelihood even after many gradient steps (dark and light green lines in panel A), but a qualitatively similar stochastic approximation using the MC does not have this same issue (panel B). We have adjusted the wording in the manuscript to clarify this distinction. We also have clarified panel C adding the equation for squared gradient error (see reviewer Adrb response).\n\n> The hybrid PA-MC model [...] should be explained in more detail\n\nIn this approach, we first fit a PA-c model and then use its parameter estimates to initialize the stochastic gradient descent algorithm that leverages  MC. This initialization is closer to the true optimum and thus helps to reduce the number of optimization steps needed and speeds up inference (a.k.a. a “warm start”). We have revised the manuscript to clarify this.\n\n> what is the individual contribution of [...] choice of temporal filters \n\nWe found that our proposed Laguerre basis outperforms the standard raised cosine basis consistently across datasets in a variety of settings and hyperparameter values. In addition to the comparison on simulated data (Fig. 2), we have expanded results on the hippocampal dataset (Fig.5), where the Generalized Laguerre basis achieves higher cross-validated log-likelihoods across a range of basis set sizes (2 to 7). In both cases for these data, performance improves as more basis functions are added. We also added a panel to Fig.S2 in the supplement illustrating the effect of introducing ridge regularization to fits on simulated data. \n\n> and choice of nonlinearities?\n\nAlthough the exponential inverse link provides a clear benefit in terms of closed-form solution for the PA-c model, we find that the softplus nonlinearity generally results in more stable inference for both PA-c and MC. Additional evaluations using softplus can be found in the supplement.\n\n>  In general, how much prior is required when using this method on a real dataset?\n\nWe set the window for filter selection that reflects realistic synaptic dynamics, a transient synaptic delay followed by a fast interaction of 1-2 ms. While these properties are common across systems, we recommend using prior literature describing timescales of the interactions to guide this choice or running preliminary analyses with longer history windows before refining. \n\n> The paper referenced supplement material several times, but I do not see any appendix. \n\nThe supplementary material is available as a downloadable .zip file included with the submission. The folder contains a PDF supplement, along with GPU-optimized model code and an example script. Please let us know if any part of the material is inaccessible.\n\n**Additional references**\n\n[Chen2009] Chen, Z., Putrino, D. F., Ba, D. E., Ghosh, S., Barbieri, R., & Brown, E. N. (2009). A regularized point process generalized linear model for assessing the functional connectivity in the cat motor cortex. IEEE\n\n[Owen2013] Owen, A. B. (2013). Monte Carlo theory, methods and examples. \n\n[Pillow&Scott2012] Pillow, J., & Scott, J. (2012). Fully Bayesian inference for neural models with negative-binomial spiking. NIPS 2025"}, {"review_id": "HtWZxCNiU3", "reviewer": "Reviewer_Adrb", "summary": "The manuscript proposes two new algorithms to optimize continuous-time Poisson Generalized Linear Models (GLM) on spike train data. This allows to infer dynamics at submillisecond timescales and scale this to long spike trains without the memory overhead that is introduced by binning the data into spike counts. For this the authors approximate the cumulative intensity function (CIF) that is part of the continuous-time GLM either using Monte Carlo (MC) or polynomial approximations (PA) of the integral. Using Laguerre polynomials as a basis function for the latter even admits a closed form solution. The method is demonstrated both on real and synthetic data, showing the proposed algorithms scale better in the number of neurons and time-steps. Using their algorithm they find functional connections in the hippocampus that are in line with what we know experimentally.", "questions": "- Do the authors know why the error in Fig. 2 is increasing with more > 3 basis functions and is this something expected? If yes, this is not obvious to me. Could the authors elaborate on this? Would tuning using a Bayesian information criterion make sense?\n- Something that is also not obvious to me is how the stratified sampling leads to lower variance compared to uniform sampling in L158? If the authors could explain this to me, that would be much appreciated.", "limitations": "My feeling is that the authors insufficiently discuss potential shortcomings of their method.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "## Strengths\n- Using a continuous-time formulation of the GLM can take advantage of the sparse nature of spike trains. Instead of operating on discretized spike count arrays of size (1, num_time_bins), one can operatore on lists of spike times [t_1,t_2,...,t_N] whose size is independent of the discretization. While the idea to make these models tractable using polynomial approximations is not new, the authors propose to use a basis of Laguerre polynomials that permit a closed form solution and need fewer features to capture neural filters. Both of this improves computational efficiency, however accuracy does not seam to be effected (according to Fig. 4 and 5). Using a MC approximation is also a straight forward thing to try and gives the experimenter more options when deploying these methods. Both of these are small, but useful contributions in my view.\n- The fact that the authors are able apply this method on hippocampal data and show method can recover known connectivity structures is also a nice demonstration of why I think this method could generally be useful.\n\n## Weaknesses\n- The paper overall feels a bit rough around the edges in its current state and it can be hard at times to follow the thought-process of the authors. While the method and its derivation are sound, the presentation could in my opinion be improved. I found that a lot of information is left for the reader to piece together (see comments and as follows). For example, I find the there could be a more natural flow between paragraphs of the experiments section that makes it easier to follow. In addition I had trouble understanding the exact content of Fig. 3, despite the associated paragraph and figure caption. What is the point it is making? Why is the error increasing after optimization in panel C? If the authors could explain this to me again / make this more clear, I would be appreciate this.\n- Another issue I have with this method is that related work, while being mentioned, is barely discussed and no shortcomings are given.\n\n## Comments\n- L4 why not write Monte Carlo methods?\n- On my printed out version of the manuscript it is hard to tell apart the similar shades of pink and red within the same plot (I.e. Fig. 1 B and D, Fig. 4,). It would be great if the authors could make them more differentiable.\n- L142 The use of W here is a bit confusing, since it means something different from the W used elsewhere. Renaming this to $\\tau$ or $\\Delta T$ similar could make this less confusing.\n- What is $I$ in Eq. 10? Was this explained?\n- While I figured out what DB means in Fig. 3 from Fig 5.'s accompanying text it would be great to mention this somewhere earlier in the manuscript (or pot. the figure caption).\n- In Fig. 4 you plot convergence time, but while you do mention what this means in App. 1.4., this is not linked\n- The use and introduction of abbreviations (i.e. MC) is inconsistent and confusing. Sometimes abbrevs. are used, sometimes not, sometimes they are introduced multiple times. \n- The supplement is sometimes referenced, but it would be nice if the authors could ref. the specific sections directly, i.e. in L273 or in L188.", "was_revised": true, "final_justification": "The authors have addressed all my concerns during the rebuttal. I slightly err on the side of acceptance and have therefore raised my score from 3 to 4.", "rebuttal": "We thank the reviewer for their thoughtful questions, comments, and suggestions. We have made revisions in the manuscript to incorporate their feedback, and we elaborate on their criticisms and requests for clarification below:\n\n### **Global Response**\n\n**Comparison to existing literature/baseline evaluation**\n\nAll reviewers note that there is prior work on continuous GLMs and note a lack of comparison to these methods. However, much of the cited work are not actual model fits of a continuous GLM, and thus are not appropriate baselines for direct comparison. Some citations develop theoretical tools relevant to continuous Poisson point processes—such as convexity of the log-likelihood [1], or error bounds and information capacity [14]—and others use continuous-time models in latent variable settings [16]. Some of the other work mentioned, particularly [15] and [Chen2019] referenced by 8hfS are actually discrete-time implementations of GLMs, and so are captured in our discrete benchmark in figures 3, 4, and 5. \n\nTo our knowledge, the only prior work that actually fits a continuous time Poisson GLM is [17] which uses Gauss-Lobatto quadrature to approximate the CIF. Unfortunately, the link provided in their manuscript to public code is no longer available. We attempted to implement their method but found that a naive implementation was both computationally infeasible on datasets of the size we consider and less accurate in approximating the CIF and fitting the model on very small subsets of data. Their approach relies on inserting a varying number of quadrature nodes between every spike time, which, for our hippocampal dataset with ~5 million spikes, would require storing and evaluating far more nodes than spikes—leading to significant memory issues and prohibitively slow inference (MC, by contrast, uses many times fewer samples than spikes). We have since tested their quadrature approach on dataset smaller than our current evaluations in the manuscript and found that the CIF estimate as well as estimated parameter values were less accurate than both our MC and PA approaches even when we use a large number of nodes. We believe this is due to the high-frequency content of the filters, which cannot be captured well by standard quadrature schemes.\n\nNevertheless, to address this concern directly in the manuscript, we will expand our related work section to clarify that nearly all GLM implementations are in discrete time. Regarding [17], although a full reproduction is not feasible, we will include a brief comparison in the supplementary material based on our partial reimplementation on a small simulated dataset, highlighting both the computational challenges and performance limitations we observed. This will help clarify the practical barriers to applying existing quadrature-based methods. We will also note that developing more scalable quadrature methods— those better suited for GPU acceleration or for estimating high-frequency coupling filters—remains a direction for future work.\n\n**Limitations of our approach**\n\n\nWe agree with the reviewers that the limitations of our approach were not sufficiently discussed, and we have expanded the conclusion section of the manuscript to address this. Many limitations of our model are inherited from the broader class of Poisson GLMs. This includes the challenge of dissociating monosynaptic connections from correlated firing (e.g. [6]) and the difficulty of identifying true connectivity without overly penalizing weak dependencies and connections with low spiking neurons ([24], [31], see response to reviewer 8hfS). Addressing these limitations is a priority for future work. Another limitation is the choice of the Poisson distribution itself, which may be a suboptimal match for describing spiking activity due to its assumptions about variance. Flexible alternatives such as the negative binomial distribution can potentially better capture the characteristics of neural spiking (e.g. [19] and [Pillow&Scott2012]). \n\nIn response to z5Zq’s and 8hfS’s comments about individual contributions and shortcomings of our MC and PA approaches, we note that there is a clear trade-off between speed and performance. As 8hfS noted, the PA method is inherently less accurate due to its global approximation of firing rates. However, it enables much faster inference (especially with the closed-form solution). The sampling-based (MC) approach, while more accurate overall, can potentially suffer from slow inference since the number of required samples M grows proportional to the recording length. This limitation is partially overcome in the hybrid model through PA-based initialization. \n\n**Notation**\n\nWe thank the reviewers for drawing our attention to the notation errors. We have made the necessary revisions to ensure that our use of notation and abbreviations is introduced clearly and consistently in the correct order. This includes consistent use of $\\textbf{w}$ and $x_t$ in Eq. 1, 2 and throughout the text; a more thorough introduction of the notation used in Equation 4; and consistent use of abbreviations (e.g., MC, GLM, DB). We have also aligned the notation used in the PA-c description with that of the MC method.\n\n### **Reviewer Specific concerns**\n> I find there could be a more natural flow between paragraphs of the experiments section\n\nWe have revised the manuscript to improve the overall flow of the experiments section. Specifically, we added a heading to the subsection with the analysis presented in Fig.3 (Stochastic gradient variance in discrete and continuous GLM)  and added clearer transitions at the beginning of the simulated and real data experiments.\n\n>I had trouble understanding the exact content of Fig. 3\n\nWe apologize for the confusion here. The key point of this section is that the variability of the gradient estimates in the discrete GLM is exceptionally high when we are required to update the model in batches. Thus, gradient-based optimization converges very slowly. In contrast, using stratified MC samples yields much more stable inference with much better log-likelihood values. This is because for our MC approach, the first term in the log-likelihood is computed exactly and the CIF integral is sampled over the entire recording. The DB method however performs updates on a small subset of data, resulting in high gradient inaccuracy, as shown in Fig.3 A and C..\n\nWe believe this confusion may be caused by an error in the caption to Fig.3C. This panel does not show the variance of the gradient estimates, but rather the expected squared error between the true and stochastic gradients or, in other words, the variability of the stochastic gradients around the true values. This error increases toward the end of training, since accurately estimating increasingly small gradient steps becomes more difficult as the models approach convergence. We clarify this distinction in the supplemental material and have updated the manuscript text accordingly. We also added the equation for the squared gradient error displayed in panel C: $\\mathbb{E}\\frac{||\\nabla_p-\\tilde{\\nabla}_p||^2_2}{||\\nabla_1||^2_2}$ where $\\nabla_p$ is the true gradient at parameter update p and $\\tilde{\\nabla}_p$ is the corresponding stochastic gradient.\n\n> Comments **(additional notation points)**\n\nWe are grateful for the detailed feedback on notation. We have incorporated the necessary changes into the manuscript, including updating the notation for $W$ (now $H$) and $\\mathbf{W}$ (now  $\\mathbf{w}$), introducing $I$ before Eq.10, and clarifying the wording in line 4. We also added references to relevant sections of the supplementary material where appropriate and fixed citation [30]. We also adjusted the colors used for purple, red, and pink to improve their distinguishability in print, and we rearranged the panels in Figure S2 to a horizontal layout to eliminate excess white space.\n\n> why the error in Fig. 2 is increasing with more > 3 basis functions \n\nThe increase in error observed for both GL and RC bases on the simulated dataset is due to overfitting and was seen specifically for these simulated filters (the introduction of ridge regularization changes this trend somewhat for these, but does not improve overall fits). \nHowever, this is not true for all datasets. We repeated the analysis from Fig. 2 on a subset of the real hippocampal dataset from Fig. 5 and found that, for both basis choices, the log-likelihood on held-out data improved as more basis functions were added. In this case, we use cross-validated log-likelihood as an evaluation metric, though BIC would also be a valid choice. Moreover, the GL basis outperformed the RC basis across all set sizes. We will include these results as an additional panel in Fig. 5 and discuss these points in the manuscript. \n\n> how the stratified sampling leads to lower variance compared to uniform sampling\n\nThe intuition behind stratified sampling is that it ensures a more uniform coverage of the CIF domain, avoiding the clustering and gaps that can occur with purely uniform samples. For example, if we drew M=10 samples, we might get unlucky and all 10 samples would clump together at a similar time. Intuitively, we do better by dividing the recording into 10 equally sized time intervals and drawing one sample per interval. This forces the samples to be more spread out. For more details, see Chapter 8 in [Owen2013]. \n\n**Additional references**\n\n[Chen2009] Chen, Z., Putrino, D. F., Ba, D. E., Ghosh, S., Barbieri, R., & Brown, E. N. (2009). A regularized point process generalized linear model for assessing the functional connectivity in the cat motor cortex.  IEEE \n\n[Owen2013] Owen, A. B. (2013). Monte Carlo theory, methods and examples. \n\n[Pillow&Scott2012] Pillow, J., & Scott, J. (2012). Fully Bayesian inference for neural models with negative-binomial spiking. NIPS 2025."}, {"review_id": "biSg6FBEGt", "reviewer": "Reviewer_8hfS", "summary": "This paper presents novel scalable inference methods for functional neural connectivity using continuous-time Poisson Generalized Linear Models (GLMs), directly addressing the limitations of traditional binned GLMs. Conventional GLMs struggle to capture submillisecond synaptic dynamics due to coarse binning and become computationally prohibitive for large-scale recordings.\n\nThe authors introduce two primary approaches that operate in continuous time: a Monte Carlo (MC) sampling-based method and a polynomial approximation (PA) method. Both methods demonstrate improved accuracy and computational tractability by avoiding the need for a large design matrix. A significant contribution is the proposal of exponentially scaled Laguerre polynomials as orthogonal temporal basis functions, which enhance filter identification and yield closed-form integral solutions for the PA approach.\n\nValidated on synthetic and real rodent hippocampal data, these continuous-time methods achieve high accuracy and scalability, enabling functional connectivity inference with submillisecond temporal precision. The recovered coupling filters align well with empirical cross-correlograms and established anatomical properties of hippocampal subregions, providing new opportunities for neural circuitry identification.", "questions": "- It would be helpful to include more comparisons or discussion with other continuous-time GLM approaches.\n- Maximum likelihood estimation in autoregressive GLMs is often prone to instability issues, such as self-excitation. Does the proposed method encounter similar challenges?\n- How does the overall firing rate influence the performance of the method?", "limitations": "Yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths \n- The work directly tackles the issues of coarse temporal resolution and scalability inherent in traditional binned Poisson GLMs, which struggle to capture submillisecond synaptic dynamics and become computationally prohibitive for large-scale recordings due to the need for a large design matrix.\n- High accuracy and temporal precision: The continuous-time approaches demonstrate superior accuracy and scalability in filter recovery, enabling accurate identification of monosynaptic coupling.\n- Efficient basis functions: exponentially scaled Laguerre polynomials as orthogonal temporal basis functions, which improve filter identification, provide closed-form integral solutions for the PA approach, and better capture biophysical synaptic dynamics compared to raised-cosine bases.\n\nWeaknesses\n-   The polynomial approximation methods (both PA-d and PA-c) are inherently less accurate compared to Monte Carlo (MC) or hybrid approaches. \n-   While the Monte Carlo (MC) approach improves stability compared to discrete batched GLMs, the paper identifies \"exploring variance reduction techniques for Monte Carlo sampling of the CIF\" as a key future direction. \n-   The current work primarily focuses on identifying functional connectivity using an exponential nonlinearity and restricts its analysis to spike history. \n-   Related work: existing work on continuous-time GLM e.g. https://doi.org/10.1109/IEMBS.2009.5334610 is not discussed in the paper.", "was_revised": true, "final_justification": "I am mostly satisfied with the response, and will maintain my score. \nThe paper successfully showcases high accuracy and temporal precision, efficient basis functions, and tackles the scalability challenge for modern large-scale neural recordings. Reproducibility and detailed experimental setup are well-covered. \nThe remaining concern is the limited comparison or discussion with other existing continuous-time GLM approaches. A more explicit contextualization would be valuable. The authors acknowledged self-excitation instability in autoregressive GLMs applies, though overall stability was observed in simulations.", "rebuttal": "We thank the reviewer for their time and effort reviewing our submission, as well as their helpful comments.\n\n### **Global Response**\n\n**Comparison to existing literature/baseline evaluation**\n\nAll reviewers note that there is prior work on continuous GLMs and note a lack of comparison to these methods. However, much of the cited work are not actual model fits of a continuous GLM, and thus are not appropriate baselines for direct comparison. Some citations develop theoretical tools relevant to continuous Poisson point processes—such as convexity of the log-likelihood [1], or error bounds and information capacity [14]—and others use continuous-time models in latent variable settings [16]. Some of the other work mentioned, particularly [15] and [Chen2019] referenced by 8hfS are actually discrete-time implementations of GLMs, and so are captured in our discrete benchmark in figures 3, 4, and 5. \n\nTo our knowledge, the only prior work that actually fits a continuous time Poisson GLM is [17] which uses Gauss-Lobatto quadrature to approximate the CIF. Unfortunately, the link provided in their manuscript to public code is no longer available. We attempted to implement their method but found that a naive implementation was both computationally infeasible on datasets of the size we consider and less accurate in approximating the CIF and fitting the model on very small subsets of data. Their approach relies on inserting a varying number of quadrature nodes between every spike time, which, for our hippocampal dataset with ~5 million spikes, would require storing and evaluating far more nodes than spikes—leading to significant memory issues and prohibitively slow inference (MC, by contrast, uses many times fewer samples than spikes). We have since tested their quadrature approach on dataset smaller than our current evaluations in the manuscript and found that the CIF estimate as well as estimated parameter values were less accurate than both our MC and PA approaches even when we use a large number of nodes. We believe this is due to the high-frequency content of the filters, which cannot be captured well by standard quadrature schemes.\n\nNevertheless, to address this concern directly in the manuscript, we will expand our related work section to clarify that nearly all GLM implementations are in discrete time. Regarding [17], although a full reproduction is not feasible, we will include a brief comparison in the supplementary material based on our partial reimplementation on a small simulated dataset, highlighting both the computational challenges and performance limitations we observed. This will help clarify the practical barriers to applying existing quadrature-based methods. We will also note that developing more scalable quadrature methods— those better suited for GPU acceleration or for estimating high-frequency coupling filters—remains a direction for future work.\n\n### **Reviewer Specific Concerns**\n\n> The polynomial approximation methods (both PA-d and PA-c) are inherently less accurate compared to Monte Carlo (MC) or hybrid approaches.\n\nThe reviewer is correct that the PA method is inherently less accurate due to its global approximation of firing rates. However, it enables much faster inference (especially with the closed-form solution), making it a good initialization procedure in the hybrid PA-MC model. In contrast, the sampling-based (MC) approach, while more accurate overall, can potentially suffer from slow inference since the number of required samples M grows proportional to the recording length. \n\n> the paper identifies \"exploring variance reduction techniques for Monte Carlo sampling of the CIF\" as a key future direction\n\nWe highlight variance reduction for MC as an important direction for future work to address the speed limitation described above. In fact, we are currently working on a control variate technique (e.g., Chapter 8 in [Owen2013]), which significantly reduces the variance of the gradient estimates and decreases the required number of samples, resulting in faster optimization.\n\n> The current work primarily focuses on [...]  using an exponential nonlinearity \n\nWhile we focus on the exponential nonlinearity in the main text, we also use softplus nonlinearity in the supplement. In general, both the MC and PA approaches are compatible with any nonlinearity that does not violate the convexity of the log-likelihood [1]. \n\n> …and restricts its analysis to spike history.\n\nThe input covariates X can be easily extended to include common external variables such as animal position or visual stimulus presentation. This would involve an additional term in the log-rates that include new parameters plus some observed experimental covariate like position or stimulus type. Here, we restrict our analyses to spike history since we are focused on these fast transient dynamics specially as a proxy to anatomical connectivity. \n\n> Maximum likelihood estimation in autoregressive GLMs is often prone to [...] self-excitation\n\nThis phenomenon does indeed apply here. However, for the data simulation process, we ensure that the (autoregressive) self-filters are strictly negative to prevent this very self-excitation. We did not observe any self-excitation issues in the MLE estimates of post-spike filters, in simulated or real datasets.\n\n> How does the overall firing rate influence the performance of the method?\n\nThe effect of overall firing rate in our model is the same as that observed in discrete-time GLMs. High-spiking neurons are generally easier to fit because they provide more data, and even small coupling weights are sufficient to noticeably influence the predicted firing rate. In contrast, low firing rate neurons require a large negative bias term to match their background firing rate. This suppresses the impact of inputs, meaning that only large coupling weights can produce meaningful post-synaptic effects. Overall, low-spiking neurons often exhibit lower signal-to-noise ratio and present a challenge for accurate inference of connectivity. This is an important point, and we have added a sentence to the manuscript to make this limitation clear to the reader.\n\n**Additional references**\n\n[Chen2009] Chen, Z., Putrino, D. F., Ba, D. E., Ghosh, S., Barbieri, R., & Brown, E. N. (2009). A regularized point process generalized linear model for assessing the functional connectivity in the cat motor cortex. Annual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE \n\n[Owen2013] Owen, A. B. (2013). Monte Carlo theory, methods and examples."}]

Author rebuttal:

Meta-review: This paper presents new methods for fitting continuous-time Poisson GLMs to spike train data, avoiding the limitations of binning. The authors propose two approaches—a Monte Carlo estimator and a polynomial approximation with Laguerre polynomials—that together make it possible to recover functional connectivity at submillisecond resolution while scaling to large recordings. The approach is validated on both synthetic data and hippocampal Neuropixels recordings, with open-source GPU implementations provided.

The main strengths are the clear computational benefits of the continuous-time formulation, the elegant closed-form update enabled by the Laguerre basis, and the demonstration on real data. These are practical contributions that address a timely need as datasets continue to grow in size and precision. Reviewers highlighted that the work is technically solid, the methods perform better than discrete baselines, and the release of optimized code increases the potential impact.

The initial submission had weaknesses in presentation and related work, with limited discussion of comparisons to earlier continuous-time models and some notational issues. The rebuttal addressed these points: the authors clarified prior work, added supplementary comparisons, expanded the discussion of limitations, and included results on the full dataset (N=642 neurons). After discussion, most reviewers raised their scores, with consensus that this is a valuable contribution. I recommend acceptance.
Decision: Accept (poster)