Fitting summary statistics of neural data with a differentiable spiking network simulator
Authors: Guillaume Bellec
Shuqi Wang
Alireza Modirshanechi
Johanni Brea
Wulfram Gerstner
Venue: NeurIPS 2021 Poster (2021)
Keywords: fitting
spikes
MLE
likelihood
GLM
spiking
recurrent
RNN
RSNN
neural data
surrogate gradients
straight through
Bayesian
prior
statistics
binary networks

Abstract: Fitting network models to neural activity is an important tool in neuroscience. A popular approach is to model a brain area with a probabilistic recurrent spiking network whose parameters maximize the likelihood of the recorded activity. Although this is widely used, we show that the resulting model does not produce realistic neural activity. To correct for this, we suggest to augment the log-likelihood with terms that measure the dissimilarity between simulated and recorded activity. This dissimilarity is defined via summary statistics commonly used in neuroscience and the optimization is efficient because it relies on back-propagation through the stochastically simulated spike trains. We analyze this method theoretically and show empirically that it generates more realistic activity statistics. We find that it improves upon other fitting algorithms for spiking network models like GLMs (Generalized Linear Models) which do not usually rely on back-propagation. This new fitting algorithm also enables the consideration of hidden neurons which is otherwise notoriously hard, and we show that it can be crucial when trying to infer the network connectivity from spike recordings.

Fitting summary statistics of neural data with a differentiable spiking network simulator
Guillaume Bellec", Shuqi Wang", Alireza Modirshanechi, Johanni Brea", Wulfram Gerstner
Laboratory of Computational Neuroscience
École polytechnique fédérale de Lausanne (EPFL)
first. lastname@epfl.c..
Abstract
Fitting network models to neural activity is an important tool in neuroscience. A popular approach is to model brain with probabilistic recurrent spiking net- work whose parameters maximize the likelihood of the recorded activity. Although this is widely used, we show that the resulting model does not produce realistic neural activity. To correct for this, we suggest to augment the log-likelihood with terms that measure the dissimilarity between simulated and recorded activity. This dissimilarity is defined via summary statistics commonly used in neuroscience and the optimization is efficient because relies on pack-propagation through the stochastically simulated spike trains. We analyze this method theoretically and show empirically that generates more realistic activity statistics. We find that it improves upon other fitting algorithms for spiking network models like GLMs (Generalized Linear Models) which do not usually rely on back proppaatpttnnn This new fitting algorithm also enables the consideration of hidden neurons which is otherwise notoriously hard, and we show that can be crucial when trying to infer the network connectivity from spike recordings.
Introduction
Modelling neural recordings has been fundamental tool to advance our understanding of the brain. It is now possible fit recurrent spiking neural networks (RSNNs) recorded spiking activity 2 4 5). The resulting network models are used to study neural properties 6
7 8 or
reconstruct the anatomical circuitry of biological neural networks [4 [0]
Traditionally a biological RSNN is modelled using a specific Generalized Linear Model 11 (GLM) often referred as the Spike Response Model (SRM) 122 model are fitted to data with the maximum likelihood estimator (MLE). The MLE is consistent, meaning that if the amount of recorded data becomes infinite it converges parameters when they exist. However, when fitting neural activity in practice the MLE solutions are often reported generate unrealistic activity [14 15 [6) showing that this method is perfect despite it's popularity in neuroscience. We also find that these unrealistic solution emerge more easily when hidden neurons outside of the recorded units have substantial impact on the recorded neurons, This is particularly problematic because the likelihood is not tractable with hidden neurons and it raises the need for new methods to tackle the problem.
11 13. The parameters of this RSNN
the true network
To address this, we optimize sample-and aoalll--nd-mmssur loss functions addition to the likelihood: these loss functions require sampling spiking data from the model and measuring the dissimilarity between recorded and simulated data. To measure this dissimilarity we suggest to compare summary statistics
"qqual contributions. Senior authors.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

popular in neuroscience like the peristimulus time histogram (PSTH) and the correlcornl (NC). Without hidden neurons, this method constrains the network to generate realistic neural activity but without biasing the MLE solution in the theoretical limit of infinite data. In practice, it leads to network models generating more realistic activity than the MLE. With hidden neurons, the sample-and-measure loss functions can be approximated efficiently whereas the likelihood function is intractable. Although recovering the exact network connectivity from the recorded spikes remains difficult [O], we show on artificial data that modelling hidden neurons in this way is crucial to recover the connectivity parameters. In comparison, methods like MLE which ignore the hidden activity wrongly estimate the connectivity matrix.
In practice the method is simple to optimize with automatic differentiation but there were theoretical and technical barriers which have prevented earlier attempts. The first necessary component is to design an efficient implementation of back pnooppropggttion in stochastic RSNN inspired by straight-
through gradient estimators 17 18 and numerical tricks from deterministic RSNNs [19]. Previous
generative models of spiking activity relying on back- propppropngttion used artificial neural networks 220 21 which are not interpretable model in terms of connectivity and neural dynamics. Previous attempts to include hidden neurons in RSNN models did not rely on back-prop but relied on expectation maximization ITI [22] or reinforce-style gradients [3] 24 25 26 which have a higher variance PT There exist other methods to fit neural data using back- pock-propagation and deep learning frameworks but they do not back- propagate through the RSNN simulator itself, rather they require to engineer and train a separate deep network to estimate posterior distribution [28] [9] 30 [1] or as a GAN discriminator 32 B3 21. The absence of discriminator in the aaaplll-addmmssuse loss function connects it with
other simple generative techniques used outside of the context of neural data 34 35 36).
Our implementation of the algorithm is published online openly
2 A recurrent spiking neural network (RSNN) model
We will compare different fitting techniques using datasets of spiking neural activity. We denote a tensor of KP recorded spike trains as x² E f0, 10,1).2 1] xTxny where ny is the total number of visible neurons recorded simultaneously and T is the number of time steps. To model the biological network which produced that activity, we consider a simple model that can capture the recurrent interactions between neurons and the intrinsic dynamics of each neuron. This recurrent network contains ny H neurons connected arbitrarily and split into a visible and a hidden population of sizes my and ny- Similarly [T 37 [4] we use GLM where each unit is modelled with a SRM neuron [12] with Ur.j being the distance to the threshold of neuron j and its spike Zt,j is sampled at time step from a Bernoulli distribution B with mean (ux,s) where is the sigmoid function. The dynamics of the stochastic recurrent spiking neural network (RSNN) are described by: Ut.j Uthr
25,5 ~ B (o
with
(1)
Uthr
Wd =, Cori++00 +,, Cr.s
VtJi
where W defines the spike-history and coupling filters spanning dmax time-bins, b defines the biases, Uthr = 0.4 is a constant, and C a spatio-temporal stimulus filter processing a few movie frames and implemented here as convolutional neural network (CNN) (this improves the fit accuracy as seen in [14 38 and Figure(ST] from the appendix). Equations ① and (2) define the probability (²(x) of simulating the spike trains& with this model and é represents the concatenation of all the network parameters (W, and the CNN parameters). Traditionally the parameters which best explain the data given by the MLE: argmax a Po(xD) 7 5 4. When all neurons are visible, the likelihood factorizes as II, P= therefore the log likelihood can be written as the negative cross-entropy (CF) between z," and (uP) = Po(=? = 1|=P,-1) where uP is computed as in equation (2) with = =D. So when all neurons are recorded and =0 are provided in the dataset the computation of the MLE never needs to simulate spikes from the model and it is sufficient to minimize the loss function:
LMLE = - log Po CE (u ())
(3)
"Code repository: htttp:://gitee wwt/ bellec- wang- 2021-sample and-andd-aass measure

---

A
B
C
T
LSM
PD
x DD
(())
P
MLL
Pe
summary T
RENN model
RENN modals unraasisi
Region where Veluee
Loss ctors represented as distances
RS4N solution minimazing EALER
Distrbbtiio is generating the same stat stics
RSW solution mimmating EALE4775T
D
Racter pict recorded data
E
Rester piot simulated with the MLE salution
F
Raster PSTH solution
Time
Time 15
Time 1s1
Figure I: A) The distribution Ps represents the RSNN model and PD represents the true biological network, the goal is to bring Pe close to Pp. The loss function LLE is represented as the distance between Pe and Pd because is equal up to constant k ( (Pr, Pa). We draw the space of summary statistics T to represent the loss function LISTH as the distance between the statistics T(z) simulated from the RSNN model Pe and measured from the data Pp. B) Even the model is misspecified and Ps cannot perfectly match the true distribution PD. when 4 PSTH is minimized the statistics T are indistinguishable between simulation and data (the green dot lands in the dark green area). When minimizing MLE alone, the solution (black dot) might generate unrealistic activity. C) When the RSNN model expressive enough to represent the true distribution but the data is insufficient, there is some uncertainty about the true network parameters (black ellipse): minimizing LM favours solutions with realistic PSTH. D) Spikes recorded simultaneously from the V VI-dataset [39]. E) Simulated spiking response from an RSNN minimizing LMLE. F) Same as E but with 4MLL+PSTH-
3 The sample-and-measure loss functions
In this section we describe the sample-and -ed-uuur loss functions which include simulated data in a differentiable fashion in the optimization objective, direct benefit is to enable the consideration of hidden neurons. We define the aosell--nd-mmann loss functions as those which require sampling spike trains z E 1] X from the model Pa and measuring the dissimilarity between the recorded and simulated data. This dissimilarity is defined using some statistics T(z) and the generic form of the sample-and aoaele-and-masurre loss functions is:
Lsm a( T(zP), Ep, [T()))
(4)
where dl a dissimilarity function, like the mean mnea--squreed error or the cross entropy. To compute the expectations Ep, we use Monte-Carlo estimates from batch of simulated trials z. For example to match the PSTH between the simulated and recorded data, we consider the statisties T(z),,; =
1 K Σx 2 and evaluate the expectation with the unbiased estimate o,, = K Σ. alut,). Denoting the PSTH of the data as xP. K Σx Zt, sample-and- sumple-and-messure loss function for the PSTH:
k.D and choosing to be the cross-entropy, we define the
LPSTH =CE(=P,0)
When all neurons are visible, we minimize the loss function LALEESM = where PALLE ASM >> 0 are constant scalars. When there are hidden neurons. the log-likelihood is intractable. Instead we minimize the negative of a lower bound of the log- likelihoooo (see appendix D for derivation inspired by 40
24 251):
LELEO CE(2P ((("))
(6)
with (u²) being the firing probability of the visible neurons, where the visible spikes 2 are clamped to the recorded data z" and the hidden spikes z are sampled according to the model
3

---

dynamics, Hence the implementation of LELBO and LSM are very similar with the difference that the samples used in L are not clamped (but all our results about LsM are also valid when they are clamped if we use the extended definition given in Appendix D). To compute the gradients with respect to these loss functions we use propaapopno which requires the propagation of gradients through the stochastic samples z. If they were continuous random variables, one could use the reparametrization trick [27], but extending this to discrete distributions is harder [17 [8 41 42 43
One way to approximate these gradients is to relax the discrete dynamics into continuous ones [4]
or to use relaxed control variates 42], but we expect that the relaxed approximations become more distant from the true spiking dynamics the network architecture gets very deep or if the network is recurrent. Instead, we choose to simulate the exact spiking activity in the forward pass and use straight-through gradient estimates [17 [8] in the backward pass by defining pseudo -eruvvvverv dzui over the binary random variables žt, ~ B(o(u,,))) We use here the same pseudo-derivative Bugg
8 = y max(0,1 - |u,,e) as in deterministic RNNs 19 because the dampening factor (here dhe.s
= 0.3) can avoid the explosive accumulation of approximation eπrors through the recurrent dynamics 44 Although the resulting gradients are biased, they work well in practice.
A geometrical description of a sample-and-measuee melsure loss function In the remaining paragraphs of this section we provide a geometrical representation and a mathematical analysis of the loss function LsM. For this purpose, we consider that the recorded spike trains z² are sampled from an unknown distribution Pp and we formalize that our goal is to bring the distribution Pe as close as possible to Pp. In this view, we re-write LSM = d(Ep,, [T(2)] Ep,. [=))) and we re-interpret LMLE the Kullback-Leibler divergence (Dx,.) from Pe Pp. This is equivalent because the divergence DKL(PD. Po) is equal to -Ero [log (=(= (=P) up to a constant.
In Figure D we represent the losses LILE and LPSTH in the space of distributions and we can represent LIL DKL((p Pa) the distance between PD and Pg. To represent the sample- and-measure loss function LSM (or specifically LPSTH in Figure T), we project the two distributions onto the space of summary statisties T represented in light green. Hence, these projections represent the expected statistics Ep, [T()) and Epp (T(xP) and LsM can be represented as the distance between the two projected statistics.
Although minimizing LLE should recover the true distributions (i.e. the biological network) if the recorded data is sufficient and the model is well specified, these ideal conditions do not seem to apply in practice because the MLE solution often generates unrealistic activity (see Figure [D-E). Fig.re(D-EE. Panels B and C in Figure I Jsee the geometrical representation of panel A to summarize the two main scenarios where minimizing LMIE+SM is better than L MML alone. In panel B, we describe a first scenario in which the model is misspecified meaning that it is not possible to find d² that Pp = Pa-. In this misspecified setting, there is no guarantee that the MLE solution yields truthful activity statistics and it can explain why the MLE solution generates unrealistic activity (Figure [[].. In this case, adding sample-and-measure loss function can penalize unrealistic solutions to solve this problem (Figure T and F). Another possible scenario is sketched in panel C. It describes the case where the model is well specified but Cas LE is flat around Pd for instance because too few trials are recorded or some neurons are not recorded all. In that case we suggest to minimize LMLE+SM to nudge the solution towards another optimum where LMLE is similarly low but the statistics T match precisely. In this sense, LsM can act similarly as a Bayesian log-prior to prefer solutions producing truthful activity statistics.
Theoretical analysis of the sample-and-measur loss function To describe formal properties of the sample-and-measun mee-uus loss function Lsas, we say that two distributions are indistinguishable according to the statistics T if the expectation E. [T()) is the same for both distributions. We assume that the dissimilarity function d(T, ,T') reaches a minimum if and only if T T (this true for the mean- -qqares error and the cross-entropy). Then for any statistics T and associated dissimilarity function d we have:
Property 1. Ifthe the RSNN model is expressive enough so that there exists parameters &e for which Po and pD are indistinguisabble according to the statistics T, then 6i is a global minimum of Lsas- Reciprocally, if this minimum is reached then Pe and PD are indistinguishable according to T-

---

This property is direct consequence of our assumption on the function If T measures the PSTH means that the optimized simulator produces the same PSTH mea- sured in the data. This can be true even if the model misspecified which is why we represented in Figure TB that the RSNN minimizing MLLE+PPT lands in the dark green region where the PSTH of the data is matched ac- curately. As one may also want to target other statistics like the noise correlation (NC), it is tempting consider different statistics T1 and Tz with corresponding dissim- ilarity functions d1 and d2 and to minimize the sum of the two losses LsM,+ &M, = Pr LsM, + (12 Ms, where J1, M2 0 are constant scalars. Indeed if d1 and d2 follow the same assumption as previously, we have (see Figure for illustration):
PD x
Distributions matching the true PSTH Distributions matching the true NC
RSNN solution minimizing LYSTH RSNN solution minimizing Lvc RSNN solution minimizing LPSTH+NG
Figure Combining LPSTH and LNC-
Property If the RSNN model is expressive enough so that there exists & for which PD and Ps are indistinguishable according both statistics Ti and Tz. then de is global minimun for SS Reciprocally, if this minimum is reached, PD and Ps are indistinguishable according to Ti and T2.
This is again a direct consequence of the assumptions on d1 and d2. Additionally Figure IC conveys the idea that L and LMLE are complementary and LsM be interpreted log-prior. This interpretation is justified by the following Property which is similar to an essential Property of Bayesian log-priors. It shows that when it is guaranteed recover the model by minimizing LLEE. minimizing the regularized likelihood LMLE+SM will also recover the true model.
Property 3. If the RSNN is well specified and identifiable so that Pd Ps- and in the limit of infinite data, then the global minimum of MML++SM exists, is unique and equal to 6*.
To prove this, we first note that all the conditions are met for the consistency of MLE so o" is the unique global minimum of LMLE. Also the assumption PD Pa- is stronger than the assumption required in Properties]] [ 6 only needed to match summary statisties: with parameters o² the model is indeed matching any statisties) it is also global minimum of LsM- As a consequence it provides a global minimum for the summed loss L.M.E This solution is also unique because it has to minimize LALE which has a unique global minimum.
may seem that those properties are weaker than the classical properties GLMs: in particular CSM+MLL MLE is not convex anymore because of LSM and the optimization process is not guaranteed to find the global minimum. This could be a disadvantage for LMLE+S but we never seemed to encounter this issue in practice. In fact, as we argue later when analyzing our simulation results, the optimum of LMLE+SM found empirically always seem to be closer to the biological network than the global minimum of LM We think that it happens because the conditions for the consistency of the MLE and Property (identifiability and infinite data) are not fully met in practice. On the contrary, the Properties and hold with misspecified model a limited amount of recorded data which may explain the success of the sample-and-mease loss functions in practice. See Figure for a geometrical interpretation.
4 Numerical simulations without hidden neurons
For our first quantitative results we consider a single session of in-vivo recordings from the primary visual cortex of a monkey watching repetitions of the same natural movie [39]. We refer to this dataset as the VI-dataset. It contains the spike trains of 69 simultancously recorded neurons for 120 repetitions lasting 30 seconds. We only keep the first 80 repetitions in our training set and 10 repetitions are used for carly-stopping. Performances are tested on the remaining 30. In our first numerical results we do not include hidden neurons.
To illustrate that minimizing LMLE alone does not fit well the statistics of interest, we show in Figure -CC the learning curves obtained when minimizing LMLES LPSTH and LMLE+PSTH. We evaluate the PSTH correlation between simulated and recorded activity every time the training loss function reaches a new minimum. With MLE in Figure 3 B, the PSTH correlation saturates at sub-optimal level and drops unexpectedly when LMLE decreases. In contrast, with the sample-and- measure loss function, the PSTH correlation improves monotonously and steadily (see Figure 3A).
5

---

A
Training Plat
B
Training MLEE
C
Training Plat MLE+
0.4
PSTH Nuason's Comeoocin
Putsonn' Comslatton
Cigrittaingg
6.42
Eage pogcttannng
Enga
2393
4291
2903
4000
9000
2,00
a000
6000
D
E
Noise Correlation matrices and R on the test set
Recorded
MLE Nº 009
R² -33
PSTH+NGust PLE+ R² 0.73
1153
1330
Figure 3: Learning curves and NC performance summary on the VI-dataset. A-C) Negative log- likelihood (i.e. LMLE) and PSTH correlation on the training set for three networks trained with CPSTH. LMLE and LML++PSTH The red curves represent the correlation between the PSTH of the recorded and simulated data. To show how the training loss influences the resulting PSTH correlation we plot a new point of the red curve only when the training loss reaches a new minimum. The vertical blue line represents the best network achieving the lowest validation losses (i.e. CPSTH- LMLE and LMLE+PSTH for the three corresponding plots). D) Noise-correlation (NC) matrix as recorded in the data. X- and y-axis represent the neuron identities. E) NC matrices when the spikes are simulated from the model, it uses the same colorbar as in the recorded NC and the simulated ones, more results are provided in Table [
D). R²2 values are reported to compare
In Figure 3C, one sees that minimizing LMLE+SM produces low values of LMLE and maximizes efficiently the PSTH correlation.
We then fit simultaneously the PSTH and the noise-correlation (NC) the V1-dataset. The NC matrix is complementary to the PSTH and it is used regularly to measure the fit performance [21 37 14]. Its entries can be viewed the neuron pair and j as the correlation of their activity. Concretely it is proportional to the statistics T(=),, KT Σx, (=,,
measure of functional connectivity, and each coefficient is defined for
(=+, Zt,)) where žt, is the PSTH (see appendix C for details). Therefore the natural sample-and-mease
loss function for NC is the mean-squared error between (a
Fe.))- We denote the resulting loss function as LNCMSS We also tested an alternative loss LNG which KDD K uses the cross entropy instead of mean-squared error and compares: T(= )ks KT Ex. 24, Ftj with the Monte- Carlo estimate KT
the coefficients T(z"),, )s, and the Monte-carlo estimates KT (a
(o
We compare quantitatively the effects of the loss functions LMLE, LYSTH, LNC and LNCWE and their combinations on the Vl-dataset. The results are summarized in Table NC matrices are shown in Figure E. The network fitted solely with LPSTH shows the highest PSTH correlation while its noise correlation is almost zero everywhere (see Figure 3E) but this is corrected when adding Lve or LNCMEEE In fact network minimizing LMLE alone yields lower performance than minimizing LPSTH+NG for both metries. When combining the log-likelihood on the test set is not compromised and it fits better the PSTH and the NC: the coefficient of determination R² of the NC matrix improves by a large margin in comparison with the MLE solution. Analyzing the failure of MLE we observe in Figure that the NC coefficients are overestimated. We wondered if the fit was mainly impaired by trials with unrealistically high activity in Figure [.. But that does not seem to be the case, because the fit remains low with MLE (R² -0 0.78) even when we discard trials where the firing probability of a neuron is higher than 0.85 for 10 consecutive time steps.
losses into LMLL+
or
SMLE+PST++NCHPL
We report Figure
PSTH correlation and noise-correlation in a different format to enable qualitative comparison with the results obtained with spike-GAN on the same dataset (sec

---

Table 1: Performance summary on the test set when fitting RSNN models to the VI-dataset. The precise definition of the performance metrics are given in Appendix C The standard deviation across neurons provided for the PSTH correlation. The variability across different network initialization is relatively low in comparison with the difference across algorithms, for instance we computed the standard deviation of LALLE over 3 seeds for MLLE and MI.E+PSTH+NC MSE and found respectively 2- 1 and 3 10 4. For the noise correlation R², the standard deviation was
9 100 for MLE+PSTH+NC
Noise Correlation (R²)
Negative log- -igllikeli
Method
PSTH correlation
LMLE (on test set)
MLE
0.67 0.16
-0.91
0.370
PSTH
0.72 + 0.15
-3.3
0.44
MST
0.69 + 0.15
0.80
0.50
MLE+PSTH+NC
0.69 + 0.15
0.73
0.370
Figure 3B,C from (2I). The fit is qualitatively similar even if we do not need separate discriminator network. Also our RSNN model is better suited to make interpretations about the underlying circuitry because it models explicitly the neural dynamics and the recurrent interactions between the neurons (whereas generic stochastic binary CNN without recurrent connections was used with the spike-GAN).
We also compare our approach with the 2-step method which is a contemporary alternative to MLE for fitting RSNNs [4]. The PSTH and noise correlation obtained with the 2-step method were measured on another dataset of 25 neurons recorded in-vitro in the retina of the Rat [9] We trained our method on the same dataset under the two stimulus conditions and quantitative comparison is summarized in Table S6 Under a moving bar stimulus condition we achieve a higher noise correlation (3% increase) and higher PSTH correlation (19% increase). But this difference might be explained by the use of a linear-simulus filter [14] instead of a CNN. Under a checkerboard stimulus condition, the 2-step method was tested with CNN but we still achieve better noise-correlation (5% improvement) with slightly worse PSTH correlation (2% decrease). Another difference is that it is not clear how the 2-step method can be extended to model the activity of hidden neurons done in the following section.
In summary, this section shows that using differentiable simulator and simple sample-and-ne - loss functions leads to a competitive generative model of neural activity. The approach can also be generalized to fit single-trial statistics as explained in the Appendix [] and Figure S3
5 Model identification
Beyond simulating realistic activity statistics, want the RSNN parameters to reflect truthful anatomical circuitry or realistic neural properties. To test this, we consider a synthetic dataset generated by a target network for which we know all the parameters. We build this target network by fitting it to the VI-dataset and sample from this model a synthetic dataset of similar size as the VI-dataset (80 training trials of approximately 30 seconds). Since our target network can generate as much data as we want, we simulate a larger test set of 480 trials and a larger validation set of 40 trials. We then fit student networks on this synthetic dataset and compare the parameters é of the student networks with the ground truth parameters o² of the target networks.
Well specified model without hidden neurons As first comparison we consider the simplest case where the target network is fully observed: the target network consist 69 visible neurons and each student network is of the same size. This is in fact the ideal setting where the log-likelihood is tractable and the MLE enjoys strong theoretical guarantees. In particular if the CNN weights are not trained and are copied from the target-network, the loss function LMLE is convex with respect to the remaining RSNN parameters Φ and the target network is identifiable [45) The resulting fitting performance is summarized in Figure 4A where we show the NC matrix and the connectivity matrix
(E, W for the target network and two students networks. We do not show the PSTH because
all methods already fit it well on the VI-datasct (see Table I). In this setting, combining Lvc and
7

---

A
Identifiable setting
C
Misspecified setting
E Modelling the hidden activity
Recorded teacher netwo & with 65 nourons
Studant network with 69 neurers
Teacher with recorded and 421 hidden neurons
Studient nabwork relurons
Stucent metwork, with some hidden neurens INLE not aspilicatll
D
B
D
:..9
R
R²
Noise Correi Coreistion
01
-0.0
-00
0.65
11
n² 0.0.5
18 :.55
Connectvity Matrix
r9.55
0.25
0.00
0.00
0.2
025
MLE++P5TH-NCC
ELBO+ PSST--NGux
(ours)
Figure 4: A, B) Results of the model identification experiment in the fully identifiable setting. A summarizes the simulation setup. The results shown in B: in the first row, we show the noise correlation matrix. In the second row, we show the connectivity matrix Ed Wd where the x- axis indicates pre-synaptic neuron and y-axis indicates post-synaptic neuron j. C, D) Same as A, B but in a misspecified setting: the target network has 500 neurons and student networks have 69 neurons. E, F) Same as before with the same target network as hidden neurons. When averaging R2 of the connectivity matrices over 5 different student network initialization, we find R2 0.56 0.0069 for ELBÖ+SM-h and 0.63 + 0.013 for ELBO+SM- h*PSTH+NCMsSR For the noise correlation, we find 0.66 + 0.1 for ELBO+SM-1 and 0.95 + 0.2 for ELBOOSMM-A+STHHNNC.cc
but the student has some
Here number of hidden neurons the student network same as in the target network, but this is not necessary as seen in Table S4 The connectivity matrices are displayed with higher resolution in Figure S5 and S6
LPSTH with LALE brings almost no advantage: the MLE already provides good reconstruction of the NC and connectivity matrices.
Model misspecification when ignoring hidden neurons From these results hypothesize that this fully identifiable setting does not capture the failure of MLE observed with real data because the recorded neurons embedded in a much larger biological network that we cannot record from. To model this, we construct another synthetic dataset based on a larger target network of 500 neurons where the first 69 neurons are fitted to the neurons recorded in the VI-dataset and the remaining 431 are only regularized to produce realistic mean firing rate (see appendix for simulation details). As in the standard setting where one ignores the presence of hidden neurons, we first consider that the student networks model only the first 69 visible neurons. This model is therefore misspecified because the number of neurons are different in the target and student networks, hence this setting well described by the scenario sketched in Figure TB.
The results are shown in Figure 4B. We found that MLE much worse than the sample-and- -easure method in this misspecified setting and the results resemble better what has been observed with real data. With MLE the roise-correlate poor (negative R²), but it significantly improves after adding the sunple-and-measun loss functions (R2 0.95). This suggest that ignoring the impact of hidden neurons can explain the failure of MLE experienced in the real V VI-dataset. We find little relationship between the student and teacher connectivity matrices (only the connectivity between visible neurons are compared, see Figure This suggests that the standard strategy, where the hidden neurons are ignored, is unlikely to be informative about true cortical connectivity.
coefficient are over estimated and the overall fit is rather
Well specified model with hidden neurons To investigate whether including hidden neurons leads to more truthful network models, we take the same target network of 500 neurons and fit now student

---

networks of the same size (500 neurons) but where only the first 69 are considered visible (Figure CC). Since the model is well specified but data about the hidden neurons is missing, this experiment is well summarized by the scenario of Figure IC. We use LELBO for the visible units and we add sample-and-measure loss function LSM-N to constrain the average firing rate of the hidden neurons which are completely unconstrained otherwise (see appendix). As seen in Figure 4C, it yields more accurate NC matrix (R²² = 0.49) and connectivity matrix (A²² = 0.55) compared to the previous misspecified models which did not include the hidden neurons. When we add sample-and-measuie loss functions to fit the PSTH and NC of the visible neurons, the noise-correlation matrix and connectivity matrix are fitted even better (R² 0.97 and R 0.63). Quantitatively, the R for the connectivity matrices are almost as high as in the easy setting of panel A where the target network is fully-visible and identifiable. Although the student network had exactly the same number of hidden neurons in Figure 4 E-F, the success is not dependent on the exact number of hidden neurons as shown in Table S4 Rather, assuming small hidden population size of only 10 neurons was enough to alleviate the failure observed in the absence of hidden neurons in Figure 4 C-D. Table [S4] also shows however that the accuracy of the reconstruction improves substantially if the hidden population is made larger in the student network.
6 Discussion
We have introduced the
method for fitting an RSNN to spike train recordings. This method leverages deep learning software and prop-propn for stochastic RSNNs to min- imize sample-and-measee and efficiently the activity of hidden neurons. We have shown that this is important to reconstruct trustworthy connectivity matrices in cortical areas. We believe that towards better models with neuroscientifically informed biases to reproduce accurately the recorded activity and funetional connectivity. Although we have focused here GLMs, PSTH and NC, the method is applicable to other spiking neuron models and statistics.
loss functions. A decisive feature of this method is to model simply
approach paves the way
Perspective One of the promising aspects of our method is to fit models which are much larger. One way to do this, is to combine neurons from separate sessions in a single larger network by considering them alternatively visible hidden. This problem natural to implement this with our method and we believe that it is a decisive step to produce models with dense coverage of the recorded areas.
tackled partially in [46] (7). It is
To investigate if our method viable in this regime we carried out a prospective scaling experiment on a dataset from the Mouse brain recorded with multiple Neuropixels probes across 58 sessions 48] The goal of this scaling experiment is only to evaluate the amount of computing resources required to fit large networks. We ran three fitting experiments with 2, 10 and 20 sessions respectively. Focusing on neurons from the visual cortices, it yielded models with 527, 2219 and 4995 neurons respectively. Each simulation could be run single NVIDIA V100 GPU and running 100 training epochs took approximately 4, 12 and 36 hours respectively. We conclude that this large simulation paradigm is approachable with methods like ours and we leave the fine-tuning of these experiments and the analysis of the results for future work.
Acknowledgments and Disclosure of Funding
This research was supported by Swiss National Science Foundation (no. 200020_ 1146011 and the Intel Neuromorphic Research Lab. Many thanks to Christos Sourmpis, Gabriel Mahuas, Ulisse Ferrari, Franz Scherr and Wolfgang Maass for helpful discussions. Special thanks to Stéphane Deny, Olivier Marre and Ulisse Ferrari for sharing with us the Retina dataset and to Matthew Smith and Adam Kohn for making their dataset publicly available.
Authors contributions GB and SW conceived the project initially. SW did most of the simulations under the supervision of GB. All authors contributed significantly to the theory and the writing.

---

References
[1] Jonathan W Pillow, Jonathon Shlens, Liam Paninski, Alexander Sher, Alan M Litke, EJ Chichilnisky, and Eero P Simoncelli, Spatio-temporal correlations and visual signalling in a complete neuronal population. Nature, 454((207):995-999) 2008
[2] Caroline A Runyan, Eugenio Piasini, Stefano Panzeri, and Christopher D Harvey. Distinct timescales of population coding across cortex. Nature, 548(7665);92-66, 2017.
[3] Rajeev Rikhye, Aditya Gilra, and Michael M Halassa. Thalamic regulation of switching between cortical representations enables cognitive flexibility. Nature neuroscience, 21(12):1753 1763, 2018.
[4] Felipe Gerhard, Tilman Kispersky, Gabrielle J Gutierrez, Eve Marder, Mark Kramer, and Uri Eden. Successful reconstruction of physiological circuit with known connectivity from spiking activity alone. PLoS Comput Biol, 9(7)ee003 2013.
[5] Ryola Kobayashi, Shuhei Kurita, Anno Kurth, Katsunori Kitano, Kenji Mizuseki, Markus Diesmann, Barry J Richmond, and Shigeru Shinomoto. Reconstructing neuronal circuitry from parallel spike trains. Nature communications, 10(1):1 13, 2019.
[6] Skander Mensi, Richard Naud, Christian Pozzorini, Michael Avermann, Carl CH Petersen, and Wulfram Gerstner. Parameter extraction and classification of three cortical neuron types reveals two distinct adaptation mechamisms, Journal of neurophysiology, 107(6):175--117) 2012.
[7] Christian Pozzorini, Richard Naud, Skander Mcnsi, and Wulfram Gerstner. Temporal whitening by power-law adaptation in neocortical neurons. Nature neuroscience, 16(7):942 948, 2013.
[8] Corinne Teeter, Ramakrishnan Iyer, Vilas Menon, Nathan Gouwens, David Feng. Jim Berg, Aaron Szafer, Nicholas Cain, Hongkui Zeng, Michael Hawrylycz, et al. Generalized leaky integrate-and-fire models classify multiple neuron types. Nature communications, 911::1115, 2018.
[9] Stephane Deny, Ulisse Ferrari, Emilie Mace, Pierre Yger. Romain Caplette, Serge Picaud, Gašper Tkačik, and Olivier Marre. Multiplexed computations in retinal ganglion cells of a single type. Nature communica-- rions, 8(1):1-17, 2017.
[10] Abhranil Das and Ila Fiete. Systematic errors in counectivity inferred from activity in strongly recurrent networks. Nature Neuroscience, 23 23(10): 1206::222222:1 1296, 2020.
[11] Jonathan W Pillow and Peter Latham. Neural characterization in partially observed populations of spiking neurons. Ady Neural Information Processing Systems, 20(3.5), 2008.
[12] Wulfrrm Gerstner. Time structure of the activity in neural network models. Physical review E, 51(1):738, 1995.
[13] Wulfram Gerstner, Werner M Kistler, Richard Naud, and Liam Paninski. Neuronal dynamics: From single neurons to networks and models of cognition. Cambridge University Press, 2014.
[14] Gabriel Mahuas, Giulio Isacchini, Olivier Marre, Ulisse Ferrari, and Thierry Mora. A new inference approach for training shallow and deep generalized linear models of noisy interacting neurons. Advances in neural information processing systems, 2020.
[15] David Hocker and Memming Park. Multistep inference for generalized linear spiking models curbs runaway excitation. In 2017 8th International IKEF//MBS Conference on Neural Engineering (NER), pages 613-616. IEFE, 2017.
[16] Felipe Gerhard, Moritz Deger, and Wilson Truccolo. On the stability and dynamics of stochastic spik- ing neuron models: Nonlinear hawkes process and point process glms. PLoS computational biology. 13())ee005330, 2017.
[17] Yoshua Bengio, Nicholas Léonard, and Aaron Courville. Estimating or propagating gradients through stochastic neurons for conditional computation. arXiv preprint 130.:33338.332 2013.
[18] Tapani Raiko, Mathias Berglund, Guillaume Alain, and Laurent Dinh. Techniques for learning binary stochastic feedforward neural networks. arXiv preprint 1406 2989, 2014.
[19] Guillaume Bellec, Darjan Salaj, Anand Subramoney, Robert Legenstein, and Wolfgang Maass, Long short- term memory and Saminngtt-leeae in networks of spiking neurons. In Advances in Neural Information Processing Systems, pages 787-797, 2018
[20] David Sussillo, Rafal Jazefowicz, L. Abbott, and Chethan Pandarinath. I.fads-latent factor analysis via dynamical systems. arXiv preprint arXiv:1608.00 06315. 2016.
[21] Poornima Ramesh, Mohamad Atayi, and Jakob II Macke. Adversarial training of neural encoding models population spike trains. NeurlPS workshop 2019: Real Neurons and Hidden Units, 2019.
[22] Vernon Lawhern, Wei Wu, Nicholas Hatsopoulos, and Liam Paninski. Population decoding of motor cortical activity using a generalized linear model with hidden states. Journal of neuroscience methods. 189(2):267222) 280, 2010.
10

---

[23]
Johanni Brea, Walter Senn, and Pean-Pascal Pfister. Sequence learning with hidden units in spiking neural networks. Advances in neural information processing systems, 0111222-4300 2011.
[24) Johanni Brea, Walter Senn, and Jean-Pascal Pfister. Matching recall and storage in sequence learning with spiking neural networks. Journal of neuroscience, 33(23):9565 00025595555555 2013.
[25] Danilo Jimenez Rezende and Wulfram Gerstner. Stochastic variational learning in recurrent spiking networks. Frontiers in computational neuroscience, 8:38, 2014.
[26] Diego M Arribas, Yuan Zhao, and II Memming Park. Rescuing neural spike train models from bad mle. NeurIPS 2020, 2020.
[27] Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra. Stochastic hackpropagation and approxi- mate inference in deep generative models. arXiv preprint /4XX::4401,4482 2014.
[28] George Papamakarios and lain Murray. Fast e-free inference simulation models with bayesian conditional density estimation. In Advances in Neural Information Processing Systems, pages 1028-1036, 2016.
[29] Jan-Matthss Lueckmann, Pedro Goncalves, Giacomo Bassetto, Kaan Öcal, Marcel Nonnenmacher, and Jakob H Macke. Flexible statistical inference mechanistic models of neural dynamics. Advances in neural information processing systems, pages 1289 1299, 2017.
[30] Pedro J Gonçalves, Matthis Lueckmann, Michael Deistler, Marcel Nonnenmacher, Kaan Öcal, Giacomo Bassetto, Chaitanya Chintaluri, William Podlaski, Sara A Haddad, Tim P Vogels, al. Training deep neural density estimators to identify mechanistic models of neural dynamics. bioRxiv, page 838383, 2020.
[31] Sean R Bittner, Agostina Palmigiano, Alex T Piet, Chunyu A Duan, Carlos D Brody, Kenneth D Miller, and John P Cunningham. Interrogating theoretical models of neural computation with deep inference. bioRxiv, page 837567, 2019.
[32] Ian Goodfellow, Jean Abadie, Mehdi Mirza, Bing Xu, David Warde- Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. In Advances in neural information processing systems, pages 2672-2680, 2014.
[33] Manuel Molano-Mazon, Arno Onken, Eugenio Piasini, and Stefano panzeri. Synthesizing realistic neural population activity patterns using generative adversarial networks. arXiv preprint arXiv: 1803.0033333 2018.
[34] Yujia Li, Kevin Swersky, and Rich Zemel. Generative moment matching networks. In International Conference on Machine Learning, pages 1718-1727. PMIR, 2015.
[35] Jesse Engel, Lamtharn Hantrakul, Chenjie Gu, and Adam Roberts. Ddsp: Differentiable digital signal processing. arXiv preprint arXiv:2001, 04643, 2020.
[36] Kyle Cranmer, Johann Brehmer, and Gilles Louppe. The frontier of simulation-based inference. Proceed- ings of the National Academy of Sciences, 1 13744)):0005-33002, 2020.
[37] Jakob H Macke, Lars Buesing, John P Cunningham. Byron M Yu, Krishna V Shenoy, and Maneesh Sahani. Empirical models of spiking in neural populations. Advances in neural information processing systems. 14:1350-1358, 2011.
[38] Lane T McIntosh, Niru Maheswaranathan, Aran Nayebi, Surya Ganguli, and Stephen A Baccus. Deep learning models of the retinal response to natural scenes. Advances in neural information processing systems, 29:1369, 2016.
[39] Matthew A Smith and Adam Kohn. Spatial and temporal scales of neuronal correlation in primary visual cortex. Journal of Neuroscience, 28(44): 2008.
[40] Michael I Jordan, Zoubin Ghahramani, Tommi S Jaakkola, and Lawrence K Saul. An introduction variational methods for graphical models. Machine learning, 37(2):83322) 1999.
[41] Chris J Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A continuous relaxation of discrete random variables. arXiv preprint 11Xiv:1611.0077 2016.
[42] George Tucker, Andriy Mnih Chris J Maddison, Dieterich Lawson, and Jascha Sohl- bial-Dicsstein. Rebar: Low- variance, unbiased gradient estimates for discrete latent variable models. arXiv preprint 1003v:/703.07377 2017.
[43] Emre O Neftci, Hesham Mostafa, and Friedemann Zenke. Surrogate gradient learning in spiking neural networks. IEEE Signal Processing Magazine, 36:61-63,, 2019.
[44] Yoshua Bengio, Patrice Simard, and Paolo Frasconi. Learning long-term dependencies with gradient descent is difficult. IEEE transactions on neural networks, 5(2):157 1994.
[45] Liam Paninski. Maximum likelihood estimation of cascade point-process neural encoding models. Network: Computation in Neural Systems, 15(4):243-262, 2004.
[46] Oleksandr Sorochynskyi, Stéphane Deny, Olivier Marre, and Ulisse Ferrari. Predicting synchronous firing of large neural populations from sequential recordings. PLoS computational biology. 17(1):1008501) 2021.
11

---

[47] Srini Turag, Lars Buesing, Adam M Packer, Henry Dalgleish. Noah Petúi,, Michael Hausser. and Jakob H Macke. Inferring neural population dynamics from multiple partial recordings of the same neural circuit. NIPS, 2014.
[48] Joshua H Siegle, Xiaoxuan Jia, Séverine Durand, Sam Gale, Corbett Bennett, Nile Graddis, Greggory Heller, Tamina K Ramirez, Hannah Choi. Jennifer A Luviano, al. Survey of spiking in the mouse visual system reveals functional hierarchy. Nature, pages 1-7, 2021.
[49] Adam Kohn and Matthew A Smith. Utah array extracellular recordings of spontaneous and visually evoked activity from anesthetized macaque primary visual cortex (v1). CRCNS. org, 2016.
[50] William F Kindel, Elijah D Christensen, and Joel Zylberberg. Using deep learning to probe the neural code for images in primary visual cortex. Journal of vision, 19(4)299-29 2019.
12


Reviews:
[{"review_id": "lnoZMOSeMSO", "reviewer": "Reviewer_nRQi", "summary": "This paper introduces a novel approach to fitting spiking network\nmodels to neural data. The proposed method uses differentiable\nrecurrent spiking networks, which are fit to the data by minimising a\nloss function with the help of recent techniques such as the\nstraight-through gradient estimate. The key novelty is to generalise\nthe loss function to be minimized, to include not only the likelihood\nof the data, but also a number of regularizing terms that depend on\ndesired spike train statistics, such as the PSTH or the noise\ncorrelations. These domain-specific terms in the loss function help\nwith breaking degeneracies in the model, and with improving fit\nquality in the face of model misspecification. After a theoretical\nanalysis of the properties of the proposed loss function(s), the paper\nshows examples of applications of the technique to some simulated\ndata. A comparison to a state-of-the-art generative model for\npopulation spike trains, and an example of usage of the proposed\nmethod on a large database of experimental recordings, are presented\nin the supplementary material.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "In my opinion, the paper brings an interesting perspective to a topic\nthat is of sufficiently broad interest within the NeurIPS community, and\nthe idea behind the proposed technique is, as far as I am aware, a\nnovel combination of known methods. The effectiveness of the\nsample-and-measure method in improving fit quality, especially in the\nface of model misspecification, is convincingly presented, and I\nenjoyed the geometrical interpretation of the technique. the authors\ndo a good job of comparing with another relevant method from the\nrecent literature, and of giving the reader an idea of how the\ntechnique could scale up to large datasets, which are becoming more\nand more common. The paper is clearly written, and I was not able to\nspot any technical issue or shortcoming. Overall, I think this is a\nsolid contribution.\n\nMy only critique to the work is that the claims on the usefulness of\nthis method for connectivity reconstruction seem not to be very\nstrongly supported by the evidence presented. More specifically, the\nclaim that \"the sample-and-measure method [...] is important to\nreconstruct trustworthy connectivity matrices in cortical areas\"\n(lines 284--288) rests upon the fact that a correctly-specified model\nwith hidden neurons can recover the connectivity matrix of the visible\npopulation with an R² \"almost as high\" (line 281) as in the case where\nthe full population is visible. I have two objections to this:\n1. the R² for the proposed method with hidden neurons is 0.64, while\n   in the fully-visible case it is 0.66. For the baseline MLE method,\n   these numbers are 0.59 vs 0.65. In an absolute sense the numbers\n   don't look that far off, and in absence of systematic statistical\n   assesments it is hard to tell if the difference between the two\n   methods is really up to the difference between them, or it depends\n   on the specific dataset/network being simulated.\n2. even taking the claim at face value, what the data in figure 4c\n   show is that connectivity can be recovered in presence of hidden\n   neurons (whether better or not than by simple MLE is up for debate,\n   see previous point), *if we work with a perfectly well-specified\n   model*. Indeed, and unsurprisingly, the only results we have for a\n   misspecified model (figure 4b) show that, even with the\n   sample-and-measure method, connectivity estimates should not be\n   trusted. But in practice, with real data, our models will be always\n   misspecified, so it is not clear how useful it is to show that\n   connectivity can be recovered in the well-specified case.\nTo summarise my two points above, I think that the good performance of\nthe model at recovering the connectivity can be explained mostly by\nthe fact that the model is well specified (and indeed, MLE does pretty\nwell too), rather than as a specific property of the proposed method.\n\n\nMinor points and typos\n- Line 220: \"instead of a CNN\" (remove \"the\").\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have adequately addressed the limitations and potential negative societal impact of their work.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "XLYfjxmixFV", "reviewer": "Reviewer_k9D3", "summary": "This paper proposes a new method to fit a recurrent spiking neural network (RNN) to neural data, which consists in simultaneously maximizing the likelihood of the spike trains observed (MLE), while fitting other observables such as averaged activity over trials (PSTH) and noise correlations (NC). On real neural data (monkey V1, rat retina), the method is shown to better capture these other observables compared to models only trained on MLE. On semi-synthetic data (artificial network fit to match the V1 dataset with additional hidden neurons), the method is shown to slightly better reflect the underlying connectivity profile compared to a method trained only on MLE.", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "__Summary of the review:__ The proposed optimization procedure to fit RNNs to neural data is sophisticated and theoretically well-justified. However the quantitative results of the method are not entirely convincing, the potential of the method to bring new insights to neuroscience is not clearly demonstrated, and the work would greatly benefit from being better tied to prior work. For these reasons, I am not fully convinced that this paper is ready for acceptance at NeurIPS in its current state.\n\n__Strengths__\n- The optimization procedure proposed is sophisticated, thoroughly described and well-justified.\n- The theoretical insight that fitting these additional observables can be seen as a Bayesian prior on the space of solutions is clear and interesting.\n\n__Limitations__\n- The quantitative benefits of the method are not entirely convincing. On real data, the method is shown to indeed improve the fit of the observables compared to a pure MLE model, which is expected, but perhaps disappointingly the MLE fit itself in not improved by the method. On synthetic data, the network connectivity (R2 = 0.64) recovered by fitting PSTH and NC is barely improved compared to a model which does not fit these additional observables (R2 = 0.59).\n- The potential of this method to bring new insights to neuroscience is not clearly demonstrated.\n- The paper would benefit greatly from being better connected to prior art: \n\nLFADS is a prior effort of fitting RNN to neural data (it would be nice to understand the difference with this work):\n\n*LFADS - Latent Factor Analysis via Dynamical Systems\nDavid Sussillo, Rafal Jozefowicz, L. F. Abbott, Chethan Pandarinath\nhttps://arxiv.org/abs/1608.06315*\n\nPrevious work has shown that the underlying connectivity of a network is very difficult to estimate from the neural activity:\n\n*Systematic errors in connectivity inferred from activity in strongly coupled recurrent circuits\nAbhranil Das,  Ila R. Fiete\nhttps://www.biorxiv.org/content/10.1101/512053v1*\n\nThere has been a lot of efforts in fitting various observables of population statistics, which are not discussed in the paper. For a review:\n\n*Modeling the correlated activity of neural populations: A review\nhttps://arxiv.org/abs/1806.08167*\n\nFor a specific example: the K-pairwise model is the maximum entropy model reproducing the mean spiking activities,\npairwise correlations, population count histogram P(K):\n\n*Tkacik, G., Marre, O., Amodei, D., Schneidman, E., Bialek,\nW., and Berry II, M. J. (2014). Searching for Collective\nBehavior in a Large Network of Sensory Neurons. PLoS\nComput Biol, 10(1), e1003408.*\n\n*Tkacik, G., Marre, O., Mora, T., Amodei, D., Berry II,\nM. J., and Bialek, W. (2013). The simplest maximum entropy model for collective behavior in a neural network.\nJournal of Statistical Mechanics: Theory and Experiment,\n2013(03), P03011.*\n\nOn fitting sequential recordings of overlapping populations (a goal mentioned in the Perspective section):\n\n*Oleksandr Sorochynskyi, Stéphane Deny, Olivier Marre, and Ulisse Ferrari. Predicting synchronous firing of large neural populations from sequential recordings. arXiv preprint arXiv:1904.04544, 2019.*\n\n__Minor__\n- It would be interesting to see a fit of the PSTHs in the plots.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"ethical_concerns": "No.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "7pOQKsZ8cw", "reviewer": "Reviewer_cqK7", "summary": "This work proposes to regularize the MLE with loss on statistics of model-simulated spike train to improve the estimation of recurrent spiking neuron network. It also introduces hidden neurons into the model in hope that it helps inference.", "questions": "", "limitations": "", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality:  Previous work was discussed in the submission. However, it loses track of some recent work on the same topic.\n- [Arribas et al. Neurips 2020] proposed similar idea to the sample-and-measure loss function. They optimize the maximum mean discrepancy (MMD) between data and simulated spike train or MMD-regularized MLE.\n- Park et al. Neural computation 2012] proposed a general kernel-based framework that defines dissimilarity and divergence between spike trains.  \n\nQuality: Generally the submission shows theoretical and empirical analysis and evaluation to support the claim. However these properties lack rigorous proof and more evaluations are needed for comprehensive comparison.\n- The model in this work looks more like a (autoregressive) GLM to me (It is called GLM in the citations e.g. [1]). Yet I will not dwell on arguing about the definition. \n- MLE of such a model in the submission is a convex problem. However the sample-and-measure loss function is not necessarily convex w.r.t. the parameters in general. Then its existence and uniqueness of the global minimum is questionable. The authors should provide rigorous proof for these properties. Practically,  even if a global minimum exists, it is still unknown if it can be reached by the optimization algorithm used in this work.\n- The V1 data used in the experiment is input-driven. It is also necessary to see the result of spontaneous activities, especially the long-term simulation.\n- In real data, the recorded neurons always have input from other sources. The hidden neurons in the model could identify them, but it can hardly be enough for single neurons. It is more important to compare the model performance with/without hidden neurons for the same data. \n- The empirical results only compared different losses on the proposed model. Comparisons with other work are needed.\n\nClarity: Overall the submission is clearly written and well organized. Please clarify the following\n- How to determine the number of hidden neurons\n- Choice of $\\mu_{MLE}$ and $\\mu_{SM}$\n- Table 1: error of NC and NLL\n\nSignificance: One issue of MLE on such models with autoregressive components is that the model is actually different for training and simulation. During training the log-likelihood is conditional on spike history that is considered fixed, but the simulation is conditional on spike history that is generated randomly (random variables). This often results in unrealistic simulation like Fig.1E. The proposed method could help to deal with this issue.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The limitations are not adequately addressed. Comparisons with other work are needed.\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "10"}, "extra_text": {}}, {"review_id": "6JAgR39Myc", "reviewer": "Reviewer_faLe", "summary": "This study proposes an approach to fitting spiking networks to neural population data, by combining a typical maximum likelihood loss (in case of a fully observed network, but an ELBO loss for a network with unobserved neurons) with an additional loss term measuring the distance between summary statistics of recorded data and simulated data. This approach is shown to capture the summary statistics of interest better than a typical maximum likelihood approach (without this additional term) in a few scenarios.\n", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "### Originality\n\nThe novelty of the work resides in the design and analysis of a new loss function for fitting spiking neural networks. Previous work is appropriately cited and it is clear how this work differs from previous contributions.\n\n\n### Quality\n\nThe paper is technically sound, but some of the claims are not sufficiently backed by empirical evidence:\n\n-the empirical results of Figures 3 and 4 seem to be for one run of the model. This is somewhat insufficient to assess the real difference in performance between losses, and increasing this to at least 10 repeats (e.g. different seeds) would be critical;\n\n-the hidden neuron scenario of Figure 4C - where the neural activity is unobserved but the exact number of neurons is known - is quite unrealistic and I would recommend the authors to test a scenario where the exact number of hidden neurons is also unknown.\n\n\n### Clarity\n\nThe manuscript is for the most part clearly written. However, it sometimes contains typos and imprecisions, such as:\n\n-instead of \"has been a fundamental tool to progress our understanding\", \"has been a fundamental tool to progress in our understanding\";\n\n-instead of \"like the peri-stimulus histogram (PSTH)\", \"like the peri-stimulus time histogram (PSTH)\";\n\n-line 67, z is missing index k in the factorisation of the likelihood;\n\n-line 83, \\mu_{MLE} and \\mu_{SM} are not defined (similarly, see line 152);\n\n-\"the visible spikes zV are clamped to the recorded data zD and the hidden spikes zH are sampled causally using our model\". The term \"causally\" is slightly confusing. I would suggest the authors to either remove the term or rephrase the sentence;\n\n-I would suggest the authors to improve the clarity of the sentence \"Although minimizing LMLE recovers the true distribution PD under ideal conditions it does not seem to apply with neural data because the MLE solution often generates unrealistic activity\";\n\n-the sentence in lines 159-162 (\"This interpretation...ideal for MLE\") is not very rigorous and hard to parse, and would benefit from a rephrasing.\n\n\n\n### Significance\n\nThe results suggest that the developed approach could be a promising step towards scalable methods for fitting spiking networks to neural population data. However, the study generally falls short of empirical evaluations of the method performance (note that the manuscript is one page less than the maximum, and would have the space for additional material).\n\n### After rebuttal\n\nThank you for the answers to my review. Given the other reviews and the authors' rebuttal, and in particular the willingness of the authors to clarify the main contributions of the paper, I increased my rating to better reflect the quality and significance of the work.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors do not sufficiently discuss the limitations of the study. In particular, the authors discuss the computational cost of the method as a function of the dimensionality of the fitted network, but do not include the performance of the method as a function of the dimensionality, which is arguably more important to report and discuss than the computational cost.\n\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: Dear authors, 

Congratulations on your paper being accepted at Neurips. Your submission was discussed extensively, (and, at times, controversially) amongst the reviewers. We appreciated that your approach has the potentially to substantially increase the usability of a common analysis tool in neuroscience, by ensuring that GLMs fitted to multi-neuron spike recordings lead more faithful models of neural firing. The current ’state of the art’ is still to use MLE, and there is mounting empirical evidence that this often yields models that do not provide accurate generative samples (e.g. exploding firing rates, or inaccurate representations of PSTHs and correlations). The approach to extend the loss function and to use differential simulators for optimisation seems to provide a clear path for overcoming this general  problem — therefore, the final decision was to accept this paper at Neurips. 

At the same time, the discussion and review process revealed multiple points in which the paper could be strengthed— please see the reviews and comments by the reviewers for specific points. However, a general theme relates to the framing of the paper— it is clear that approach leads to better generative samples, and this is demonstrated empirically. It is not immediately clear why this approach would also lead to better reconstructions of 'functional connectivity’ (aside from the many issues plaguing this concept), and the empirical evidence is similarly not completely convincing. 

We would ask the authors to take the feedback by the reviewers seriously, and adjust the paper accordingly. 

Best, 

Your AC
Decision: Accept (Poster)