Fitting summary statistics of neural data with a differentiable spiking network simulator
Year: 2022

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.

heading: 1 Introduction
text: Modelling neural recordings has been a fundamental tool to advance our understanding of the brain. It is now possible to fit recurrent spiking neural networks (RSNNs) to recorded spiking activity [1, 2, 3, 4, 5]. The resulting network models are used to study neural properties [6, 7, 8, 9] or to reconstruct the anatomical circuitry of biological neural networks [4, 5, 10].
Traditionally a biological RSNN is modelled using a specific Generalized Linear Model [11] (GLM) often referred to as the Spike Response Model (SRM) [12, 11, 13]. The parameters of this RSNN 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 to the true network parameters when they exist. However, when fitting neural activity in practice the MLE solutions are often reported to generate unrealistic activity [14, 15, 16] showing that this method is not 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 a 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.
To address this, we optimize sample-and-measure loss functions in 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
∗Equal contributions. Senior authors.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
popular in neuroscience like the peristimulus time histogram (PSTH) and the noise-correlation (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 [10], 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-propagation in stochastic RSNN inspired by straightthrough gradient estimators [17, 18] and numerical tricks from deterministic RSNNs [19]. Previous generative models of spiking activity relying on back-propagation used artificial neural networks [20, 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 [11, 22] or reinforce-style gradients [23, 24, 25, 26] which have a higher variance [27]. There exist other methods to fit neural data using back-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 a posterior distribution [28, 29, 30, 31] or as a GAN discriminator [32, 33, 21]. The absence of a discriminator in the sample-and-measure 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.

heading: 2 A recurrent spiking neural network (RSNN) model
text: We will compare different fitting techniques using datasets of spiking neural activity. We denote a tensor of KD recorded spike trains as zD ∈ {0, 1}KD×T×nV where nV 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 nV+H neurons connected arbitrarily and split into a visible and a hidden population of sizes nV and nH. Similarly to [1, 4, 37, 14] we use a GLM where each unit is modelled with a SRM neuron [12] with ut,j being the distance to the threshold of neuron j and its spike zt,j is sampled at time step t from a Bernoulli distribution B with mean σ(ut,j) where σ is the sigmoid function. The dynamics of the stochastic recurrent spiking neural network (RSNN) are described by:
zt,j ∼ B (σ(ut,j)) with ut,j = vt,j − vthr
vthr (1)
vt,j = nV+H∑ i=1 dmax∑ d=1 W dj,izt−d,i + bj + Ct,j , (2)
whereW defines the spike-history and coupling filters spanning dmax time-bins, b defines the biases, vthr = 0.4 is a constant, and C is a spatio-temporal stimulus filter processing a few movie frames and implemented here as a convolutional neural network (CNN) (this improves the fit accuracy as seen in [14, 38] and in Figure S1 from the appendix). Equations (1) and (2) define the probability Pφ(z) of simulating the spike trains z with this model and φ represents the concatenation of all the network parameters (W , b and the CNN parameters). Traditionally the parameters which best explain the data are given by the MLE: argmax φPφ(zD) [1, 7, 37, 5, 4]. When all neurons are visible, the likelihood factorizes as ∏ t Pφ(zDt |zD1:t−1), therefore the log-likelihood can be written as the negative cross-entropy (CE) between zDt and σ(u D t ) = Pφ(zDt = 1|zD1:t−1) where uDt is computed as in equation (2) with z = zD. So when all neurons are recorded and zD 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 = − logPφ(zD) = CE(zD, σ(uD)) . (3) 2Code repository: https://github.com/EPFL-LCN/pub-bellec-wang-2021-sample-and-measure

heading: 3 The sample-and-measure loss functions
text: In this section we describe the sample-and-measure loss functions which include simulated data in a differentiable fashion in the optimization objective, a direct benefit is to enable the consideration of hidden neurons. We define the sample-and-measure loss functions as those which require sampling spike trains z ∈ {0, 1}K×T×nV+H from the model Pφ 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-measure loss functions is:
LSM = d ( T (zD), EPφ [ T (z) ] ) , (4)
where d is a dissimilarity function, like the mean-squared error or the cross entropy. To compute the expectations EPφ we use Monte-Carlo estimates from a batch of simulated trials z. For example to match the PSTH between the simulated and recorded data, we consider the statistics T (z)t,i = 1 K ∑ k z k t,i and evaluate the expectation with the unbiased estimate σ̄t,i = 1 K ∑ k σ(u k i,t). Denoting the PSTH of the data as z̄Dt,i = 1 K ∑ k z k,D t,i and choosing d to be the cross-entropy, we define the sample-and-measure loss function for the PSTH:
LPSTH = CE(z̄D, σ̄) . (5)
When all neurons are visible, we minimize the loss functionLMLE+SM = µMLELMLE+µSMLSM where µMLE , µSM > 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-likelihood (see appendix D for a derivation inspired by [40, 23, 24, 25]):
LELBO = CE(zD, σ(uV)) , (6)
with σ(uV) being the firing probability of the visible neurons, where the visible spikes zV are clamped to the recorded data zD and the hidden spikes zH are sampled according to the model
dynamics. Hence the implementation of LELBO and LSM are very similar with the difference that the samples used in LSM 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 back-propagation 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, 18, 41, 42, 43]. One way to approximate these gradients is to relax the discrete dynamics into continuous ones [41] or to use relaxed control variates [42], but we expect that the relaxed approximations become more distant from the true spiking dynamics as 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, 18] in the backward pass by defining a pseudo-derivative ∂zt,i ∂ut,i
over the binary random variables zt,i ∼ B(σ(ut,i)). We use here the same pseudo-derivative ∂zt,i ∂ut,i
= γmax(0, 1 − |ut,i|) as in deterministic RNNs [19] because the dampening factor (here γ = 0.3) can avoid the explosive accumulation of approximation errors through the recurrent dynamics [44]. Although the resulting gradients are biased, they work well in practice.
A geometrical description of a sample-and-measure 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 zD are sampled from an unknown distribution PD and we formalize that our goal is to bring the distribution Pφ as close as possible to PD. In this view, we re-write LSM = d(EPD [T (z)] ,EPφ [T (z)]) and we re-interpret LMLE as the Kullback-Leibler divergence (DKL) from Pφ to PD. This is equivalent because the divergence DKL(PD,Pφ) is equal to −EPD [ logPφ(zD) ] up to a constant.
In Figure 1 we represent the losses LMLE and LPSTH in the space of distributions and we can represent LMLE = DKL(PD,Pφ) as the distance between PD and Pφ. To represent the sampleand-measure loss function LSM (or specifically LPSTH in Figure 1), we project the two distributions onto the space of summary statistics T represented in light green. Hence, these projections represent the expected statistics EPφ [T (z)] and EPD [ T (zD) ] and LSM can be represented as the distance between the two projected statistics.
Although minimizing LMLE 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 1D-E). Panels B and C in Figure 1 use the geometrical representation of panel A to summarize the two main scenarios where minimizing LMLE+SM is better than LMLE alone. In panel B, we describe a first scenario in which the model is misspecified meaning that it is not possible to find φ∗ so that PD = Pφ∗ . 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 1E). In this case, adding a sample-and-measure loss function can penalize unrealistic solutions to solve this problem (Figure 1B and F). Another possible scenario is sketched in panel C. It describes the case where the model is well specified but LMLE is flat around PD for instance because too few trials are recorded or some neurons are not recorded at 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-measure loss function To describe formal properties of the sample-and-measure loss function LSM , we say that two distributions are indistinguishable according to the statistics T if the expectation E [T (z)] 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 is true for the mean-squared error and the cross-entropy). Then for any statistics T and associated dissimilarity function d we have:
Property 1. If the RSNN model is expressive enough so that there exists parameters φ◦ for which Pφ and PD are indistinguishable according to the statistics T , then φ◦ is a global minimum of LSM . Reciprocally, if this minimum is reached then Pφ and PD are indistinguishable according to T .
This property is a direct consequence of our assumption on the function d. If T measures the PSTH it means that the optimized simulator produces the same PSTH as measured in the data. This can be true even if the model is misspecified which is why we represented in Figure 1B that the RSNN minimizing LMLE+PSTH lands in the dark green region where the PSTH of the data is matched accurately. As one may also want to target other statistics like the noise correlation (NC), it is tempting to consider different statistics T1 and T2 with corresponding dissimilarity functions d1 and d2 and to minimize the sum of the two losses LSM1+SM2 = µ1LSM1 + µ2LSM2 where µ1, µ2 > 0 are constant scalars. Indeed if d1 and d2 follow the same assumption as previously, we have (see Figure 2 for an illustration):
Property 2. If the RSNN model is expressive enough so that there exists φ◦ for which PD and Pφ are indistinguishable according to both statistics T1 and T2, then φ◦ is a global minimum for LSM1+SM2 . Reciprocally, if this minimum is reached, PD and Pφ are indistinguishable according to T1 and T2.
This is again a direct consequence of the assumptions on d1 and d2. Additionally Figure 1C conveys the idea that LSM and LMLE are complementary and LSM can be interpreted as a 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 to recover the true model by minimizing LMLE , 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 = Pφ∗ and in the limit of infinite data, then the global minimum of LMLE+SM exists, it is unique and equal to φ∗.
To prove this, we first note that all the conditions are met for the consistency of MLE so φ∗ is the unique global minimum of LMLE . Also the assumption PD = Pφ∗ is stronger than the assumption required in Properties 1 and 2 (previously φ◦ only needed to match summary statistics: with parameters φ∗ the model is indeed matching any statistics) so it is also a global minimum of LSM . As a consequence it provides a global minimum for the summed loss LMLE+SM . This solution is also unique because it has to minimize LMLE which has a unique global minimum. It may seem that those properties are weaker than the classical properties of GLMs: in particular LSM+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+SM 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 LMLE . We think that it happens because the conditions for the consistency of the MLE and Property 3 (identifiability and infinite data) are not fully met in practice. On the contrary, the Properties 1 and 2 hold with a misspecified model or a limited amount of recorded data which may explain the success of the sample-and-measure loss functions in practice. See Figure 1 for a geometrical interpretation.

heading: 4 Numerical simulations without hidden neurons
text: 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 V1-dataset. It contains the spike trains of 69 simultaneously 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 early-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 3A-C the learning curves obtained when minimizing LMLE , 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 3B, the PSTH correlation saturates at a sub-optimal level and drops unexpectedly when LMLE decreases. In contrast, with the sample-andmeasure loss function, the PSTH correlation improves monotonously and steadily (see Figure 3A).
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) on 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 as a measure of functional connectivity, and each coefficient is defined for the neuron pair i and j as the correlation of their activity. Concretely it is proportional to the statistics T (z)i,j = 1KT ∑ k,t(z k t,i − z̄t,i)(zkt,j − z̄t,j) where z̄t,i is the PSTH (see appendix C for details). Therefore the natural sample-and-measure loss function for NC is the mean-squared error between the coefficients T (zD)i,j and the Monte-carlo estimates 1KT ∑ k,t(σ(u k t,i)− σ̄t,i)(σ(ukt,j)− σ̄t,j). We denote the resulting loss function as LNCMSE . We also tested an alternative loss LNC which uses the cross entropy instead of mean-squared error and compares: T (zD)i,j = 1KT ∑ k,t z k,D t,i z k,D t,j
with the Monte-Carlo estimate 1KT ∑ k,t σ(u k t,i)σ(u k t,j).
We compare quantitatively the effects of the loss functions LMLE , LPSTH , LNC and LNCMSE and their combinations on the V1-dataset. The results are summarized in Table 1 and NC matrices are shown in Figure 3E. 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 LNC or LNCMSE . In fact a network minimizing LMLE alone yields lower performance than minimizing LPSTH+NC for both metrics. When combining all losses into LMLE+PSTH+NC or LMLE+PSTH+NCMSE the log-likelihood on the test set is not compromised and it fits better the PSTH and the NC: the coefficient of determination R2 of the NC matrix improves by a large margin in comparison with the MLE solution. Analyzing the failure of MLE we observe in Figure 3 that the NC coefficients are overestimated. We wondered if the fit was mainly impaired by trials with unrealistically high activity as in Figure 1E. But that does not seem to be the case, because the fit remains low with MLE (R2 = −0.78) even when we discard trials where the firing probability of a neuron is higher than 0.85 for 10 consecutive time steps.
We report in Figure S2 the PSTH correlation and noise-correlation in a different format to enable a qualitative comparison with the results obtained with a spike-GAN on the same dataset (see
Figure 3B,C from [21]). The fit is qualitatively similar even if we do not need a 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 a 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 [14]. 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 a quantitative comparison is summarized in Table S6. Under a moving bar stimulus condition we achieve a higher noise correlation (3% increase) and a 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 a CNN but we still achieve a better noise-correlation (5% improvement) with a 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 as done in the following section.
In summary, this section shows that using a differentiable simulator and simple sample-and-measure 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 D and Figure S3.

heading: 5 Model identification
text: Beyond simulating realistic activity statistics, we want the RSNN parameters to reflect a 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 V1-dataset and sample from this model a synthetic dataset of similar size as the V1-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 φ∗ of the target networks.
Well specified model without hidden neurons As a first comparison we consider the simplest case where the target network is fully observed: the target network consist of 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 ( ∑ dW d i,j) for the target network and two students networks. We do not show the PSTH because all methods already fit it well on the V1-dataset (see Table 1). In this setting, combining LNC and
LPSTH with LMLE brings almost no advantage: the MLE already provides a good reconstruction of the NC and connectivity matrices.
Model misspecification when ignoring hidden neurons From these results we hypothesize that this fully identifiable setting does not capture the failure of MLE observed with real data because the recorded neurons are 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 V1-dataset and the remaining 431 are only regularized to produce a 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 is well described by the scenario sketched in Figure 1B.
The results are shown in Figure 4B. We found that MLE is much worse than the sample-and-measure method in this misspecified setting and the results resemble better what has been observed with real data. With MLE the noise-correlation coefficient are over estimated and the overall fit is rather poor (negative R2), but it significantly improves after adding the sample-and-measure loss functions (R2 = 0.95). This suggest that ignoring the impact of hidden neurons can explain the failure of MLE experienced in the real V1-dataset. We find little relationship between the student and teacher connectivity matrices (only the connectivity between visible neurons are compared, see Figure 3). This suggests that the standard strategy, where the hidden neurons are ignored, is unlikely to be informative about true cortical connectivity.
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 4C). Since the model is well specified but data about the hidden neurons is missing, this experiment is well summarized by the scenario of Figure 1C. We use LELBO for the visible units and we add a sample-and-measure loss function LSM−h 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 (R2 = 0.49) and connectivity matrix (R2 = 0.55) compared to the previous misspecified models which did not include the hidden neurons. When we add sample-and-measure loss functions to fit the PSTH and NC of the visible neurons, the noise-correlation matrix and connectivity matrix are fitted even better (R2 = 0.97 and R2 = 0.63). Quantitatively, the R2 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 a 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.

heading: 6 Discussion
text: We have introduced the sample-and-measure method for fitting an RSNN to spike train recordings. This method leverages deep learning software and back-propagation for stochastic RSNNs to minimize sample-and-measure loss functions. A decisive feature of this method is to model simply 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 our approach paves the way towards better models with neuroscientifically informed biases to reproduce accurately the recorded activity and functional connectivity. Although we have focused here on GLMs, PSTH and NC, the method is applicable to other spiking neuron models and statistics.
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 or hidden. This problem was tackled partially in [46, 47]. It is natural to implement this with our method and we believe that it is a decisive step to produce models with a dense coverage of the recorded areas.
To investigate if our method is 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 on a 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_184615) 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.

Review summaries:
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.

This paper introduces a novel approach to fitting spiking network models to neural data. The proposed method uses differentiable recurrent spiking networks, which are fit to the data by minimising a loss function with the help of recent techniques such as the straight-through gradient estimate. The key novelty is to generalise the loss function to be minimized, to include not only the likelihood of the data, but also a number of regularizing terms that depend on desired spike train statistics, such as the PSTH or the noise correlations. These domain-specific terms in the loss function help with breaking degeneracies in the model, and with improving fit quality in the face of model misspecification. After a theoretical analysis of the properties of the proposed loss function(s), the paper shows examples of applications of the technique to some simulated data. A comparison to a state-of-the-art generative model for population spike trains, and an example of usage of the proposed method on a large database of experimental recordings, are presented in the supplementary material.