Scalable inference of functional neural connectivity at submillisecond timescales

We thank the reviewer for their time and effort reviewing our submission, as well as their helpful comments.

Global Response

Comparison to existing literature/baseline evaluation

All 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.

To 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.

Nevertheless, 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.

Reviewer Specific Concerns

The polynomial approximation methods (both PA-d and PA-c) are inherently less accurate compared to Monte Carlo (MC) or hybrid approaches.

The 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.

the paper identifies "exploring variance reduction techniques for Monte Carlo sampling of the CIF" as a key future direction

We 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.

The current work primarily focuses on [...] using an exponential nonlinearity

While 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].

…and restricts its analysis to spike history.

The 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.

Maximum likelihood estimation in autoregressive GLMs is often prone to [...] self-excitation

This 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.

How does the overall firing rate influence the performance of the method?

The 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.

Additional references

[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

[Owen2013] Owen, A. B. (2013). Monte Carlo theory, methods and examples.

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:

Global Response

Comparison to existing literature/baseline evaluation

All 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.

To 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.

Nevertheless, 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.

Limitations of our approach

We 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]).

In 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.

Notation

We 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 $**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.

Reviewer Specific concerns

I find there could be a more natural flow between paragraphs of the experiments section

We 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.

I had trouble understanding the exact content of Fig. 3

We 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..

We 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.

Comments (additional notation points)

We 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.

why the error in Fig. 2 is increasing with more > 3 basis functions

The 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). However, 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.

how the stratified sampling leads to lower variance compared to uniform sampling

The 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].

Additional references

[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

[Owen2013] Owen, A. B. (2013). Monte Carlo theory, methods and examples.

[Pillow&Scott2012] Pillow, J., & Scott, J. (2012). Fully Bayesian inference for neural models with negative-binomial spiking. NIPS 2025.

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.

Global Response

Comparison to existing literature/baseline evaluation All 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.

To 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.

Nevertheless, 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.

Limitations of our approach

We 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]).

In 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.

Notation

We 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 $**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.

Reviewer specific concerns

Figure 3A and B has no y-axis values, Line 237 and Line 245 references

We 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).

The hybrid PA-MC model [...] should be explained in more detail

In 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.

what is the individual contribution of [...] choice of temporal filters

We 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.

and choice of nonlinearities?

Although 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.

In general, how much prior is required when using this method on a real dataset?

We 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.

The paper referenced supplement material several times, but I do not see any appendix.

The 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.

Additional references

[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

[Owen2013] Owen, A. B. (2013). Monte Carlo theory, methods and examples.

[Pillow&Scott2012] Pillow, J., & Scott, J. (2012). Fully Bayesian inference for neural models with negative-binomial spiking. NIPS 2025

We thank the reviewer for their detailed and insightful feedback. Below are our responses to the reviewer’s comments and questions.

Global Response

Comparison to existing literature/baseline evaluation

All 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.

To 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.

Nevertheless, 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.

Limitations of our approach

We 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]).

In 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.

Notation

We 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 $**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.

Reviewer Specific concerns

How novel is the Laguerre basis? Is this not used in earlier work?

Laguerre 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.

baseline against previous continuous Hawkes process models and quadrature methods from Paninksi, Lindermann, and others

See 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.

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].

Does MC become prohibitively slow for high firing rates or when PA-c becomes memory-bound (long recordings, hundreds or thousands of neurons)?

Our 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.

The 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.

The paper […] currently evaluates only one small Neuropixels data set. How do things scale as the number of neurons increases?

In 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.

Additional References

[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

[Linderman&Adams] Linderman, S. & Adams, R. (2015). Scalable Bayesian Inference for Excitatory Point Process Networks. arXiv:1507.03228.

[Owen2013] Owen, A. B. (2013). Monte Carlo theory, methods and examples.

[Pillow&Scott2012] Pillow, J., & Scott, J. (2012). Fully Bayesian inference for neural models with negative-binomial spiking. Advances in neural information processing systems, 25.

[Solo2004] V. Solo, C. J. Long, E. N. Brown, E. Aminoff, M. Bar and S. Saha. (2004). fMRI signal modeling using Laguerre polynomials. IEEE

[Quesne2009] Quesne, C. (2009). Solvable Rational Potentials and Exceptional Orthogonal Polynomials in Supersymmetric Quantum Mechanics. SIGMA 5, 084.