Inferring the collective dynamics of neuronal populations from single-trial spike trains using mechanistic models

Multi-neuronal spike-train data recorded in vivo often exhibit rich dynamics as well as considerable variability across cells and repetitions of identical experimental conditions (trials). Efforts to characterize and predict the population dynamics and the contributions of individual neurons require model-based tools. Abstract statistical models allow for principled parameter estimation and model selection, but possess only limited interpretive power because they typically do not incorporate prior biophysical constraints. Here we present a statistically principled approach based on a population of doubly-stochastic integrate-and-fire neurons, taking into account basic biophysics. This model class comprises an idealized description for the dynamics of the neuronal membrane voltage in response to fast independent and slower shared input fluctuations. To efficiently estimate the model parameters and compare different model variants we compute the likelihood of observed single-trail spike trains by leveraging analytical methods for spiking neuron models combined with inference techniques for hidden Markov models. This allows us to reconstruct the shared input variations, classify their dynamics, obtain precise spike rate estimates, and quantify how individual neurons couple to the low-dimensional overall population dynamics, all from a single trial. Extensive evaluations based on simulated data show that our method correctly identifies the dynamics of the shared input process and accurately estimates the model parameters. Validations on ground truth recordings of neurons in vitro demonstrate that our approach successfully reconstructs the dynamics of hidden inputs and yields improved fits compared to a typical phenomenological model. Finally, we apply the method to a neuronal population recorded in vivo, for which we assess the contributions of individual neurons to the overall spiking dynamics. Altogether, our work provides statistical inference tools for a class of reasonably constrained, mechanistic models and demonstrates the benefits of this approach to analyze measured spike train data.


Introduction
Cortical computations are represented in the collective spiking activity of multiple neurons. The growing interest to uncover how these neuronal populations process and transform complex incoming information into decisions and motor actions has brought about cell-resolving activity measurements at an increasing scale and pace. Although these activity patterns can, in principle, span a high-dimensional space, often a large fraction of neural variability is captured by low-dimensional manifolds [1][2][3][4][5].
Statistical inference based on generative models is a powerful approach to interpret the measured spike train data and characterize the hidden, low-dimensional, collective dynamics [6,7]. For this purpose typically abstract, time-dependent latent variable models are fitted to the data [6,[8][9][10][11][12][13][14][15]. These approaches are well suited for quantifying the structure in the data, and benefit from statistically principled parameter estimation and model selection methods. However, their interpretive power and capacity for identifying neural mechanisms are limited as the underlying models typically do not incorporate prior biophysical constraints.
Mechanistic models of neural populations, on the other hand, involve variables and parameters that can be biophysically interpreted, and have proven useful to dissect neural dynamics. Well-known exponents are the detailed Hodgkin-Huxley-type neuron models [16,17], which can involve complex morphology and patterns of ionic currents; however, they are not well suited to analyze spike-train data, because multiple combinations of parameters give rise to the same firing patterns [18,19]. An alternative, prominent class of models with reduced complexity are spiking neuron models of the integrate-and-fire (I&F) type, which implement in a simplified manner essential biophysical constraints. These models have been advanced in recent years to accurately predict neuronal activity [20][21][22] and classify multiple neuron types [17,23,24]. They have become state-of-the-art models for describing neural activity in in-vivo-like conditions and have been applied in a multitude of studies on neural network dynamics. Yet, while I&F models are biophysically more faithful than abstract statistical models, fitting such models to multi-neuronal spike trains with account of variability and latent, low-dimensional population dynamics, is a substantial challenge.
Here we consider a population of doubly-stochastic I&F neurons to model highly non-stationary, collective spiking dynamics as typically observed in vivo. The hidden neuronal inputs, which impinge on the hidden membrane potentials, contain fast independent fluctuations that give rise to spiking variability, and slower shared variations that dominate the low-dimensional, latent population dynamics. Shared variability across observed neurons, which typically reflect only a small fraction of the local population, is captured in the model by a common, latent process rather than by putative direct coupling. We present a statistically principled approach based on derived likelihood functions to fit this type of model to single-trial spike trains and quantitatively compare different model variants, including simpler phenomenological models. Specifically, we efficiently compute the likelihood of a given spike train by exploiting analytical methods for stochastic I&F neuron models [25] combined with inference techniques for hidden Markov models [26]. This allows us to infer the model parameters by likelihood maximization, classify the latent dynamics via likelihood-based model selection, and estimate hidden time series from the time-varying probability density over the latent state.
We evaluate our approach extensively on synthetic data in terms of reconstruction of the true latent time series, classification of their dynamics, and estimation of the model parameters. We then validate our method using in-vitro ground truth recordings of cortical neurons stimulated by current signals with additive noise [27]. We successfully reconstruct the true dynamics of the signal and demonstrate improved fitting performance compared to a classical inhomogeneous Poisson point process model. Finally, we apply our approach to extracellular recordings from macaque primary visual cortex in vivo [28] to characterize the low-dimensional population dynamics and the contributions of individual neurons.

Results
Our results are structured as follows. The model is explained in section 1. In the following two sections we outline our inference methods and evaluate them exhaustively on synthetic data. For reasons of clarity and comprehensibility we first focus on single neurons in section 2 and consider neuronal populations in section 3. In section 4 we validate our approach based on in-vitro recordings. In section 5 we apply our methods to population spike train data from extracellular recordings in vivo.

Statistical modeling with doubly-stochastic I&F neurons
We consider having observed cell-resolved spike trains (ordered sets of spike times) from a population of N neurons. Such data are typically obtained from extracellular multi-electrode recordings after pre-processing that includes spike sorting [29]. The activity of each neuron (with index i) of the population shall be described by the classical leaky I&F model [30] with membrane time constant τ m , where the compound synaptic input is given by a Gaussian white noise process with mean and standard deviation σ i . This process effectively models synaptic bombardment impinging on the neuronal membrane voltage (see, e.g., [31][32][33][34]). The input includes shared, slow variations x(t) which may be caused by a common component in the external drive or by network interactions, and independent, rapid fluctuations with strength σ i .μ i is a cell-specific offset. C i quantifies the extent to which neuron i is affected by x(t).
Notably, in the context of inferring x(t) from observed activity data, C i measures how much the activity of neuron i contributes to the shared variations. In the following, we simply refer to it as coupling strength (for the coupling between the individual activity of a neuron and the shared dynamics).
Since the dynamics of the common input are not known we describe x(t) by a stochastic process. In particular, we consider two qualitatively different Markov models, an Ornstein-Uhlenbeck process (OUP) and a Markov jump process (MJP), which gives rise to two model variants. For the OUP x(t) varies continuously with time constant τ , whereas for the MJP x(t) is piece-wise constant intermittently jumping to different values with rate τ −1 . The stationary probability density of both processes is standard normal and their autocorrelation functions are identical. For further details on the generative model see Methods section 1. In Results sections 4 and 5 below we also consider a simpler, classical model for comparison, which describes the spike train of neuron i by a Poisson point process with rate exp(µ i (t)), where µ i (t) is given by Eq (1).
It is important to point out that not all model parameters need to be estimated. The membrane voltage can be scaled such that the remaining parameters of interest for inference are those for the input together with the membrane time constant (see Methods section 1). Moreover, a change of τ m can be well compensated for in terms of spiking dynamics by appropriate changes ofμ i and σ i [25]. Therefore, we focus on the following parameters for inference C i ,μ i , σ i for i ∈ {1, . . . , N } and τ , with particular interest on the coupling strengths (C i ) and the time constant (τ ).

Outline of inference approach
It is instructive to first consider a single neuron. For improved readability we omit the neuron index i and group the parameters as ϑ := {C,μ, σ, τ m } and τ . We collect the measured spike-train data in the increasing sequence of spike times t 0:K := (t 0 , . . . , t K ) and define the k-th interspike interval (ISI) by s k := t k − t k−1 . K is the number of observed ISIs. For notational ease below we use that s k , the duration of the k-th ISI, implicitly contains information about the start time and end time of the interval. As spike emission in the leaky I&F model is a renewal process the likelihood of observing a given spike train from the model is completely determined by the likelihoods of observing the constituent ISIs. To tackle the inference problem we approximate the time series of the process x(t) for each ISI by one value, x k for ISI s k , which is justified by the assumption that x(t) varies slowly, i.e., τ is large compared to the average ISI. Defining x 0 as the value of the process at t = 0 and x k as the value at t = t k for k ≥ 1 we obtain a hidden Markov model with sequence of latent states x 0:K := (x 0 , . . . , x K ). The observations are contained in the sequence s 1:K := (s 1 , . . . , s K ). The joint likelihood of observing a given spike train and realization of x(t), approximated by sequence x 0:K , is given by with effective mean input µ k = Cx k +μ. p(s k |µ k , ϑ) is the ISI probability density of the leaky I&F neuron exposed to Gaussian white noise input (with mean µ k and standard deviation σ), evaluated at ISI s k . This density can be accurately computed by solving a Fokker-Planck partial differential equation, which can be achieved numerically in efficient ways [25]. p(x k |x k−1 , τ ) is the transition probability density for x(t), from state x k−1 to state x k , which depends on the time constant τ and on the time duration between states x k−1 and x k . For uncluttered notation this duration is not explicitly indicated. The transition densities for both model variants (OUP and MJP) are known (for details see Methods section 2). p(x 0 ) is the prior probability density of the process at t = 0 < t 0 (prior to t 0 ) for which we assume its stationary distribution, a standard normal. Note that in Eq (2) we have used the renewal property of leaky I&F neurons. We obtain the likelihood of the observations by marginalizing Eq (2) with respect to x 0:K , p(s 1:K |ϑ, τ ) = p(s 1:K , x 0:K |ϑ, τ ) dx 0:K , which is efficiently performed by an iterative procedure known as forward filtering (see Methods section 3).
Parameter estimatesθ,τ are determined by maximizing the likelihood,θ,τ := argmax ϑ,τ p(s 1:K |ϑ, τ ). For this purpose it is an important advantage that we can evaluate p(s 1:K |ϑ, τ ) with high accuracy and low computational effort. For maximization we use an established simplex-based method [35] (for details see Methods section 3).
After having inferred the parameter values we compare the fitting performance of both model variants (OUP and MJP) to identify the best model and thereby classify the dynamics of the slow input variations. Specifically, we use the log-likelihood ratio (LLR) for comparison, equivalently expressed as the difference LLR(OUP, MJP) := log p OUP (s 1:K |θ,τ ) − log p MJP (s 1:K |θ,τ ), where subscripts OUP and MJP indicate the respective process. Positive values of the LLR indicate that the OUP model variant is favored, while negative values indicate that the MJP model variant better describes the data. We further compare the doubly-stochastic I&F model to a Poisson point process whose principal parameter, the rate, is given by λ(t) = exp(µ(t)) with µ(t) = Cx(t) +μ (cf. Eq (1)) and x(t) is described by an OUP or MJP as specified above. The exponential function is a typical choice for the mapping between input and output rate in Poisson models of this type, see e.g. [10,14,36]. We infer the parameters of this model similarly as for the I&F model, the only two differences are that for the Poisson model the ISI density is explicitly expressed by p(s k |µ k , ϑ) = λ −1 k exp(−s k λ k ) with λ k = exp(µ k ), µ k = Cx k +μ, and ϑ contains two parameters, C andμ, instead of four. For model comparison we use the Akaike information criterion (AIC), which takes into account both goodness of fit and complexity of a model (for details see Methods section 4).
Finally, to reconstruct the time series of the process x(t) and estimate the instantaneous spike rate of the neuron we infer p(x k |s 1:K ,θ,τ ), the probability density of the state x k during the ISI s k given the observations s 1:K and the inferred parameter values, for each k. We compute this density using an iterative technique known as backward smoothing, which in addition to the above-mentioned forward filtering method constitutes the established forward-backward algorithm for hidden Markov models [26] (see Methods section 5). Our estimate for the time series of x(t) is then obtained by the sequence of expected values ( x 0 , . . . , x K ) calculated using these densities. Additionally, we estimate the spike rate time series by the sequence (r 0 , . . . , r K ) of expected inverse mean ISIs, where the mean ISI for the k-th interval is calculated using p(s|µ k , ϑ) and the expectation is with respect to p(x k |s 1:K , ϑ, τ ) (for details see Methods section 5). Note that r k represents the instantaneous spike rate at the time that corresponds to x k .

Evaluation on synthetic data
For evaluation we consider synthetic spike-train data from numerical simulations of the underlying generative model. We assess the performance of the proposed method to reconstruct the time series of x(t), classify its dynamics and estimate the parameters (Fig 1). Example spike trains and true time series of x(t) together with the inferred sequences of conditional density over the latent state are shown in Fig 1A,B. The evolution of the density inferred from the spike train well resembles the true time series. Using the incorrect model variant for inference leads to an approximation that appears slightly less accurate.
A systematic quantification of classification accuracy for a wide range of parameter values and amount of observed data is presented in Fig 1C-E. Classification according to the LLR is correct throughout the tests. The decision becomes more obvious as the coupling strength C, the time constant τ or the length of the spike train increases. Intuitively, it is challenging to classify the dynamics if the neuron is only weakly affected by . Surprisingly, a few hundred observed spikes are already sufficient to identify the latent process. Note the small bias towards the MJP, which may be caused by the piece-wise constant approximation of the process in our inference method.
Accuracy of parameter inference is visualized in Fig 1F-H. The true values of the parametersμ, C, σ and τ are well recovered, and estimation errors decrease with increasing spike train length. Note that the membrane time constant τ m was not estimated here, but set to the true value. Setting τ m instead to a wrong value within a biologically plausible range affects the maximal likelihood only very little (see Fig S1), because a change of τ m can be well compensated for in terms of spiking dynamics by suitable changes of the input parameters (cf. [25]). We, therefore, keep τ m fixed for the rest of this study. Inference results for single neurons using synthetic data. A: Example spike train from the doubly-stochastic I&F model with OUP x(t) (red traces below) together with sequence of estimated density over the latent variable (color map) when using the OUP (center) or MJP (bottom) for inference. B: Example as in A with true x(t) given by a MJP (blue traces). C: Log-likelihood ratio (LLR, cf. Eq (4)) for different values ofμ and C with true x(t) given by an OUP (top) or by a MJP (bottom). We considered allμ, C-pairs, for which the expected spike rate remains in the interval [1,110] Hz for at least 98 % of the simulation time. This was analytically determined from the steady-state spike rate and the quantiles of the stationary distribution of µ(t). Results represent averages over 10 realizations. In D-H red color indicates that the true x(t) is given by an OUP, blue color indicates that it is a MJP. For each parametrization results from 10 realizations are shown, colored areas denote mean ± standard deviation. D: LLR as a function of time constant τ . E: LLR as a function of number of observed spikes. F: Difference between maximal log-likelihood for different values of σ and the one obtained for the true value. G: Estimated versus true parameter values ofμ (left), C (center) and τ (right). The results forμ and C correspond to the fits in C, the results for τ correspond to the fits in D. H: Relative and absolute errors between estimated and true parameter values as a function of number of observed spikes, respectively. The results correspond to the fits in E. If not stated otherwise, parameter values for generating the data wereμ = −4.7 mV ms , C = 0.6 mV ms , σ = 4 mV √ ms and τ = 500 ms; 5 × 10 3 observed spikes were used for fitting.
In sum, our method allows for accurate parameter inference and characterization of the latent dynamics based on synthetic data of single neurons. Even relatively short spike trains lead to satisfactory results.

Outline of inference approach
We now turn to a larger population of neurons as described in Results section 1. We collect the observed data in terms of (cell-resolved) ISIs of all N neurons in the single sequence s 1:K := (s 1 , . . . , s K ) which is ordered according to the last spike time of each ISI. Similarly as for single neurons we effectively approximate the time series of x(t) for each neuron and each ISI by one value, which is justified by the slowness of the process. Specifically, we define the sequence of latent states x 0: is the value of the process at the k-th spike time across the population for k ≥ 1. Let i k indicate the neuron that corresponds to the k-th ISI s k in the sequence s 1:K , ending at spike time t k . The mean input for neuron i k across ISI s k is thus approximated by µ k,i k = C i k x k +μ i k . That is, each neuron "samples" the latent state sequence x 0:K (which can contain multiple values within an ISI) in a distinct way according to its spike times. This approximation allows us to extend the approach for a single neuron (cf. Results section 2) to a neuronal population in a straightforward way. We express the joint likelihood of the observed data and the latent state sequence for this hidden Markov model as where Note that the three probability density functions on the right hand side of Eq (5) are identical to those of Eq (2) explained in Results section 2. To obtain the likelihood of observing the data from the model we marginalize Eq (5) with respect to x 0:K , and then continue analogously to the single neuron case for inference of model parameters, classification of the dynamics of x(t), reconstruction of its time series and estimation of the time-varying population spike rate (for details see Methods sections [3][4][5]. Notably, although we can accurately and rapidly evaluate the (marginalized) likelihood using numerical methods (cf. Methods sections 2 and 3), maximization with respect to 3N + 1 parameters becomes a serious challenge for a large population. Therefore, we use an efficient, approximate technique for optimization (see Methods section 3).

Evaluation on synthetic data
We now evaluate our method using simulated data from the doubly-stochastic I&F population model (Fig 2). Example population spike trains, spike rate histograms, hidden time series of x(t), and the results from our method are shown in Fig  Accuracy of parameter inference from multiple tests is shown in Fig 2F-I. The true values ofμ i , C i , and τ are well approximated, and errors between true and estimated values increase only mildly as recording time is reduced from 10 min to 2.5 min.
Computation times for parameter inference are visualized in Fig 2J. From the tests considered in this section computation time appears to increase linearly with the number of neurons in the population due to the efficient optimization scheme that allows for parallelization over observed neurons.
We next examine whether our inference method serves to capture qualitatively distinct shared dynamics. For this purpose we consider synthetic data from simulations of the population model where x(t) is described by a sine wave instead of the OUP or MJP (Fig 3). The results from these examples demonstrate that  the oscillatory dynamics of x(t) can be well recovered from a population of N = 20 neurons using our method based on either the OUP or MJP. Reconstruction accuracy deteriorates with increasing oscillation frequency; furthermore, the task becomes increasingly challenging for small populations and small neuronal spike rates. In sum, our method allows for accurate parameter inference, reconstruction and classification of the shared input dynamics, and spike rate estimation based on synthetic population data. It further serves to reconstruct oscillatory dynamics of the shared input, which demonstrates the dynamical flexibility of the doubly-stochastic I&F model.

Validation based on in-vitro ground truth recordings
We validate our method using ground truth data from patch-clamp recordings of cortical pyramidal neurons [27]. The cells received an injected current composed of a weak sinusoidal signal with strong additive noise such that neuronal spike rates oscillated roughly between 2 and 6 spikes/s (for details see Methods section 6). We fitted the doubly-stochastic I&F model using only the spike times and compared the inferred oscillatory dynamics of the mean input with the true oscillation of the current signal (Fig 4). The available data from three example recordings with different oscillation period and the corresponding results from our method are shown in Fig 4A-C. The magnitude of the rapid additive fluctuations (noise) dominates the oscillation amplitude of the injected current, so that spiking activity is only weakly modulated in an oscillatory manner. We measured the strength of this modulation by the resultant vector length (RVL) for the distribution of phases of the signal at spike times ( Fig 4A-C right). The RVL quantifies the degree of concentration of that circular distribution; a large value of the RVL close to 1 indicates a strongly peaked, narrow distribution. Our inferred sequences of conditional density over the mean input clearly exhibit oscillatory dynamics. Importantly, the respective amplitude-normalized oscillations we extracted from these inferred sequences match well with the (amplitude-normalized) true signals (for details see Methods section 6).
We next quantified the synchronization between the extracted and true oscillatory signals using the synchrony index (SI), a complex number whose radius (magnitude) indicates the strength of locking and the angle represents the average phase shift between the oscillations (for details see Methods section 6). The SI for all cells is shown in Fig 4D. The angle values are concentrated around 0 and most magnitudes are large, close to 1, which means that the oscillations are aligned locking is strong for most cells. The locking strength, which measures the quality of signal reconstruction, clearly depends on the extent of oscillatory structure in the spiking activity, as quantified by the RVL (Fig 4E). Those few cells for which locking is weak exhibit small RVL values, for RVL values 0.1 signal reconstruction is excellent. Weak locking is further reflected by a small value of inferred coupling strength C (Fig 4F).
Finally, we compared the fitting performance of the I&F model with that of an inhomogeneous Poisson process using the AIC on these data (Fig 4G, cf. Results section 2). The preferred model is indicated by the lower AIC value. For all recordings the I&F model outperforms the Poisson model according to this measure.
In this section we used the OUP for the latent process x(t) (in the I&F and Poisson models); using instead the MJP did not significantly change the results. Note that the models do not include any prior information about oscillatory dynamics.
In sum, our method allows to accurately recover the ground truth dynamics of hidden, weak input signals from spike trains of cortical neurons and outperforms a classical model-based approach on these data.

Application to in-vivo multi-electrode recordings
We now apply our method to extracellular recordings from primary visual cortex of a macaque monkey under anaesthesia [28]. The data consist of spontaneous spiking activity of a population of N = 20 single units over a duration of 10 minutes, which was split into two sets of 5 minutes each (one for fitting and one for testing; for details on preprocessing see Methods section 7). We fitted both I&F model variants, with either OUP or MJP for the shared input dynamics, to these spike trains and compared their fitting performance on the test set. The MJP is favored according to the LLR evaluated on the test data (LLR test = −443.0). Notably, a benchmark test against a Poisson point process model demonstrates a clear preference for the I&F population (LLR test = −1514.9, using the MJP in both models). A 30 s segment of the spike-train data from the test set and the results from our method are shown in Fig 5A-C. The spike trains are ordered according to the inferred value of coupling strength C i (Fig 5A). Units which strongly couple to the shared input component exhibit coordinated spiking activity whereas spiking of units with weak coupling appears more disordered. Using the inferred coupling strength C i , which quantifies the neuronal contribution to the low-dimensional population dynamics, we can thus separate "chorister" from "soloist" neurons [37] in a statistically principled way.
Population rate estimates versus empirical rate histograms are depicted in Fig 5B,C. Surprisingly, the OUP model variant yields a slightly better match in this respect despite its inferior fitting performance. Note, however, that approximating the population rate histogram is not an explicit objective of the fitting procedure.
Histograms over the sequence of expected values of x across the entire recording (which we refer to as latent state distributions) are shown in Fig 5B,C right. Note that the distributions differ from a standard normal, which is the stationary distribution of x as expected from the generative model. This indicates that the inferred low-dimensional dynamics deviate from both pure model variants for x, and it shows that our model flexibly captures different dynamics (as clearly demonstrated above, cf. Figs 3 and 4). Using only the values of the sequence of expected x at the observed spike times of a neuron we obtain N = 20 distinct histograms (one for each unit). This allows us to measure the extent to which individual neurons are tuned to the estimated shared dynamics in a way that does not directly depend on the inferred coupling strengths. Comparison between that measure of tuning and the coupling strengths enables another (indirect) validation. Indeed, units with a large value of C exhibit strong tuning, shown by distributions that are skewed towards large values of x (and clearly differ from the latent state distribution), while units with weak coupling also exhibit weak tuning, shown by distributions that are indistinguishable from the latent state distribution (Fig 5D,E).
In sum, our method allows to characterize the low-dimensional population dynamics and the contribution/coupling of individual neurons to those dynamics in a principled way.

Discussion
We presented a statistically principled approach to fit an I&F population model with doubly-stochastic input to spike-train data. The input contains neuron-specific white noise fluctuations and a shared Markovian component that dominates the low-dimensional overall population dynamics. We extensively evaluated our methodology on synthetic data, validated it using ground truth in-vitro and in-vivo recordings, and compared the I&F population to a classical Poisson point process model in terms of fitting performance. Altogether, the results demonstrated the benefits of our approach to identify and classify the latent low-dimensional population dynamics, and quantify how individual neurons couple to those dynamics.

Related Work
At the core of our methodology we exploit efficient numerical schemes to compute the ISI probability density of an I&F neuron exposed to Gaussian white noise input with high precision. The ability to rapidly compute this density allows to evaluate the likelihood of observed spike trains, which has been previously used for inference purposes using single stochastic I&F neurons [25,38] and networks [25]. The nonstationarity considered here (doubly-stochastic population model) constitutes a relevant and highly nontrivial extension. A related approach on the level of single neurons has been introduced in [33], which uses a diffusion process for the dynamics of the mean input and involves an approximation of the ISI probability density to avoid its numerical computation. An advantage of our approach in this regard is that it allows for the application of various, reasonable processes for the mean input, as demonstrated using the OUP and MJP. Note that the variance of these processes remains fixed over time, whereas that of a diffusion process increases without bounds. Importantly, we further considered neuronal populations.
Methodologically related work previously focused on the latent dynamics that underlie the spiking activity of single neurons, using Poisson point process models with continuous or jumping dynamics of the underlying rate, similar to those considered here [39,40]. In our comparisons the doubly-stochastic I&F model clearly outperformed such Poisson point process models with latent rate dynamics. Whether and in what way this could affect classification results [40] is currently not clear.
Abstract models with larger complexity, in particular, generalized linear point process models that account for the dynamics of latent variables can be designed and optimized to fit well observed spike train data [9,[41][42][43][44]. However, the dynamical flexibility of these models is typically associated with large numbers of parameters that need to be optimized; hence, this approach is prone to overfitting unless strong contraints and/or parameter regularization are enforced [7,36,[45][46][47]. An advantage of our approach in this regard is that the parameter space for optimization is comparably low-dimensional, which diminishes the risk of overfitting, without sacrificing essential aspects of neuronal spiking dynamics. Furthermore, the variables and parameters of I&F models can be interpreted in a straightforward way.

Limitations & potential extensions
In our inference method the mean input is approximated by one value for each neuron and each ISI. This approximation leads to an "event-based" (spike-based) temporal binning of the inferred sequences, as was previously used for single neurons [33]. It implies that the latent dynamics cannot be captured precisely across extended intervals without spikes.
We aimed to capture the latent low-dimensional collective dynamics using a stochastic process for the variations of shared inputs. This process describes the dynamics of the external drive and picks up effects due to synaptic interactions on the population-average level. Detailed effects of synaptic coupling between the observed neurons are not explicitly included. The assumption that inputs from unobserved neurons dominate the collective population dynamics is well justified as long as the observed neurons make up only a small fraction of the local population (which is a typical scenario) [10]. Nevertheless, the model we considered could be extended to a doubly-stochastic I&F circuit with an explicit description of synaptic interactions between observed neurons. In particular, a combination of our approach with recent inference methods for I&F circuits [25] -for example, using the inferred (sequences of) parameters of our population model for subsequent estimation of connectivity -may be beneficial for the estimation of synaptic couplings.
We used the classical leaky I&F model to describe the dynamics of the neuronal membrane voltage (with focus on spiking). Our methodology, however, can be easily extended to other I&F model variants for which spike emission is a renewal process, such as the exponential I&F model [48]. Furthermore, an absolute refractory period can be included in a straightforward way.
We considered two qualitatively different Markov processes (OUP and MJP) for the shared input fluctuations. Other processes may be used instead: for example, a differentiable process whose temporal correlation-function is described by a Matérn kernel or periodic kernel (see, e.g., [49]).
Here we considered a one-dimensional process for the dynamics of shared input variations, which dominate the population dynamics. Depending on the size of the observed population and application it may be relevant to account for multi-dimensional latent dynamics; our framework can be extended accordingly in a similar way as for generalized linear models [9].

Applicability
Using in-vitro ground truth data we have demonstrated that our approach accurately extracts the latent dynamics of input signals that are masked by strong noise, outperforming a more classical model-based method. These results support the validity of our method for application to in-vivo data that typically do not entail ground truth about the properties to be inferred. For such a scenario, using spike train data from extracellular recordings, we have shown that the inferred coupling strengths allow to distinguish between cells which clearly participate in coordinated population activity and those which spike rather independently. Such a characterization of "chorister" and "soloist" neurons has recently been shown to correlate with the underlying synaptic connectivity [37]. In future applications our method may serve, for example, to relate classification results (between model variants) and inferred parameter values across different experimental conditions.

Materials and methods 1 Generative model
We consider a population of N leaky I&F neurons driven by fluctuating inputs. The dynamics of the membrane voltage V i of neuron i (i ∈ {1, . . . , N }) are governed by the stochastic differential equation where τ m denotes the membrane time constant. The input variations are described by the time-dependent mean µ i (t) and Gaussian white noise process ξ i (t) scaled by σ i , where ξ i (t)ξ j (t + ∆) = δ i,j δ(∆) for i, j ∈ {1, . . . , N } with expectation · . The neuron fires a spike when the membrane voltage reaches the spike threshold value V s , subsequently V i is reset to the value V r . The dynamics of the mean input µ i (t) are determined by a common process x(t) through where C i denotes the strength of coupling of neuron i to x(t) andμ i is the offset for neuron i. We separately consider two processes with qualitatively different dynamics for x(t): an Ornstein-Uhlenbeck process (OUP), described by with time constant τ and Gaussian white noise process ξ(t), ξ(t)ξ(t + ∆) = δ(∆); alternatively, x(t) is described by a Markov jump process (MJP) which is piece-wise constant across intervals with exponentially distributed duration of mean τ and takes values drawn from a standard normal distribution N (0, 1), This process is also known as marked (homogeneous) Poisson point process. For both processes the stationary distribution is standard normal, lim t→∞ x(t) ∼ N (0, 1), and the autocorrelation function is given by Hence, our doubly-stochastic model incorporates fast independent input fluctuations by σ i ξ i (t) in Eq (6) with Gaussian white noise process ξ i (t) and slower shared input variations by µ i (t) via Eq (7) with common OUP or MJP x(t).
It is not meaningful to estimate all parameters of this model. A change of V s or V r can be completely compensated in terms of spiking dynamics by appropriate changes of C i ,μ i and σ i , which can be seen using the change of variablesṼ i := (V i − V r )/(V s − V r ). Consequently, we exclude V s and V r from the parameters to be inferred and instead set them to reasonable values. Furthermore, we fix τ m for most of our results, since a change of that parameter can also be well compensated for by appropriate changes of the input parameters (see Results section 2). The parameter values were τ m = 10 ms, V s = −40 mV and V r = −65 mV.

Spike-based discretization
We consider having observed the spike trains of N neurons and collect these data in a sequence of ISIs s 1:K := (s 1 , . . . , s K ) that are ordered increasingly according to the last spike time of each ISI. Under the assumption that the process x(t) evolves slowly compared to the duration of the ISIs we approximate the mean input for each neuron across each ISI effectively by one value: µ i k (t) ≈ µ k,i k = C i k x k +μ i k for all times t within the k-th ISI s k in the sequence s 1:K , where x k is the value of x(t) at the end spike time of s k and i k indicates the neuron that corresponds to that ISI. Note that in order to simplify notation we use that s k , explicitly the duration of the k-th ISI, informs also (implicitly) about the start and end times of the ISI.
Hence, for the inference problem x(t) is replaced by the sequence x 0:K := (x 0 , . . . , x K ) with x 0 := x(0) and x k := x(t k ) for k ≥ 1, t k being the k-th spike time across the population, which greatly facilitates the task.

Markov property
Due to the fact that both processes, OUP and MJP, are Markovian we can factorize the probability density of the sequence given the time constant parameter as where p(x 0 ) is the probability density for the initial process value (prior). We assume p(x 0 ) ∼ N (0, 1) which corresponds to the stationary distribution of the process. The transition probability density for the OUP is given by where ∆ k = t k − t k−1 . For notational convenience we do not explicitly indicate ∆ k in p(x k |x k−1 , τ ) and instead use that x k also contains information about the corresponding time point. For the MJP, on the other hand, we have where the probability of jumping is P jump = 1−exp(−∆ k /τ ), δ x k ,x k−1 is the Kronecker delta and p N ∼ N (0, 1). Note that P jump denotes the probability that at least one jump occurs during an interval of duration ∆ k .

Conditional likelihood
Since spike emission for an I&F neuron is a renewal process, the likelihood of the observed data s 1:K conditioned on the sequence x 0:K (which approximates the process x(t)) and neuron-specific parameter values can be factorized as with µ k,i k = C i k x k +μ i k . ϑ = {ϑ 1 , . . . , ϑ N } contains the parameters for each neuron, where ϑ i = {C i ,μ i , σ i , τ m }. Each factor on the right hand side is the ISI probability density of an I&F neuron exposed to Gaussian white noise input with constant parameter values, evaluated at one point. We can compute the ISI density p(s k |µ k,i k , ϑ i k ) efficiently with high precision for a range of values s ≥ 0 at once. This is achieved by solving a Fokker-Planck partial differential equation that describes the so-called first passage time problem for the stochastic I&F model [25,38,50]. Here we apply the finite volume solution method described in [25]. In practice, we pre-compute p(s k |µ k,i k , ϑ i k ) on a sufficiently fine grid of values for s k and µ k,i k , which is then used as a look-up table in the optimization procedure described below. For values of the mean input µ k,i k encountered during the inference procedure that are not on the grid we use linear interpolation. For additional details we refer to our Python implementation available at Github: https://github.com/neuromethods/inference for doubly stochastic IF models

Joint and marginal likelihoods
By combining Eqs (11) and (14) we obtain the joint likelihood of the observed data and the process sequence as with µ k,i k = C i k x k +μ i k . For the marginal likelihood of observing the data from the model, integration over x 0:K is required: Note the difference between this marginal likelihood and the likelihood in Eq (14) which is conditioned on knowledge of the process sequence.

Parameter estimation
Estimates for a particular set of parameters are obtained by maximizing the marginal likelihood with respect to those parameters (maximum likelihood estimation). In the following we describe how we evaluate this likelihood from Eq (16) and how we maximize it with respect to the parameters of interest. For additional details we refer to our code available at Github.

Evaluation of the marginal likelihood
Integration over x 0:K in Eq (16) can be efficiently performed by an iterative procedure. Knowing p(x k−1 |s 1:k−1 , ϑ, τ ) allows for the prediction step using the Chapman-Kolmogorov (forward) equation The iteration is initialized for k = 1 and p(x 0 |s 1:0 , ϑ, τ ) = p(x 0 ). In the following step we incorporate the next observation s k by the filtering p(s k , x k |s 1:k−1 , ϑ, τ ) = p(s k |µ k,i k , ϑ i k )p(x k |s 1:k−1 , ϑ, τ ) with µ k,i k = C i k x k +μ i k . By marginalizing out x k we obtain p(s k |s 1:k−1 , ϑ, τ ) = p(s k , x k |s 1:k−1 , ϑ, τ )dx k (19) and calculate in the last step p(x k |s 1:k , ϑ, τ ) = p(s k , x k |s 1:k−1 , ϑ, τ ) p(s k |s 1:k−1 , ϑ, τ ) , which completes the iteration k − 1 → k. After K iterations the marginal likelihood is given by Practically, the steps of the iteration described above are performed using a reasonable discretization for x; here we used bins of width 0.05 in the range [−3.5, 3.5]. Consequently, the step in Eq (17), for example, involves a square matrix for the transition probability density and a sum for the integral.

Likelihood maximization for a single neuron
We maximize the logarithm of the likelihood from Eq (21) with respect to the parameters C,μ and τ using a simplex optimization method [35]. Note that our approach is not restricted to the simplex method, other (e.g. gradient-based) techniques may also be applied. The parameter σ is estimated via an outer loop, where we iterate over a range of values for σ and maximize the likelihood for each of those. The membrane time constant τ m is not optimized as justified in Results section 2.

Likelihood maximization for a population
The parameter space for optimization has 3N + 1 dimensions, where N is the number of neurons. To reduce computational efforts and obtain parameter estimates in reasonable time we applied the following approximate optimization scheme. First, we estimate the noise amplitude σ i for i ∈ {1, . . . , N } by fitting a doubly-stochastic I&F model with independent process x i (t) to each spike train separately, as described in the previous paragraph. Then, we re-estimate the offsetμ i by fitting a model neuron with C i = 0 and σ i from the previous step for each neuron. Note that these two steps can be performed in parallel across neurons. Next, we use the full population model and maximize the likelihood with respect to all coupling strengths C 1:N := (C 1 , . . . , C N ) given τ and vice versa in an alternating way, keeping σ 1:N andμ 1:N fixed. The optimization with respect to C 1:N is done in parallel across neurons: for optimizing C i the couplings for other neurons are fixed to the values found in the previous optimization step. Finally, after convergence of the previous iteration procedure for C 1:N and τ , our estimates for the parameters C 1:N andμ 1:N are improved by maximizing the likelihood with respect to C i andμ i using the previously determined values as starting points, keeping all other parameters fixed, for i ∈ {1, . . . , N }. This last step is performed again in parallel. Although this scheme appears rather approximate it yields accurate parameter estimates (cf. Results section 3).

Classification of latent dynamics and model comparison
To classify the dynamics of the slow (shared) input variations we compare both fitted model variants, one using the OUP and one using the MJP, by applying the log-likelihood ratio (LLR). This quantity can be expressed as the difference LLR(OUP, MJP) := log p OUP (s 1:K |θ,τ ) − log p MJP (s 1:K |θ,τ ), where the subscripts OUP and MJP indicate the considered process, andθ,τ denote the parameter estimates determined prior to model comparison. Note that the LLR does not take into account the model complexity, which is not needed because both compared models have the same number of estimated parameters.
To compare models with different complexities in terms of number of parameters to be estimated, specifically the I&F and the Poisson models considered here, we apply the Akaike information criterion (AIC) [51]. This measure is given by 2Nθ ,τ − 2 log p(s 1:K |θ,τ ), where Nθ ,τ is the number of estimated parameters. The preferred model is indicated by the smallest AIC value.
The probability density p(x k |s 1:k ,θ,τ ) on the right hand side of Eq (23) is already known from the forward (filtering) iteration (cf. Eq (20)) at this point. We calculate the probability density p(s k+1:K |x k ,θ,τ ) iteratively by backward smoothing, withμ k,i k =Ĉ i k x k +μ i k . We initialize p(s K+1:K |x K ,θ,τ ) = 1 (uniform distribution), since s K is the last observed ISI. The denominator in Eq (23) is obtained by numerical marginalization. In this way we obtain the sequence (p(x 1 |s 1:K ,θ,τ ), . . . , p(x K |s 1:K ,θ,τ )) that we visualize in Figs 1-5. Using the elements from this sequence we estimate x 0:K by the sequence of expected values ( x 0 , . . . , x K ), where x 0 is calculated with respect to p(x 0 ).
The sequence of marginal densities further allows us to estimate the (instantaneous) population spike rate time series by the sequence (r 0 , . . . , r K ), ∞ 0 s p(s|μ k,i ,θ i )ds −1 (25) withμ k,i =Ĉ i x k +μ i , where the expectation · is calculated with respect to p(x k |s 1:K , ϑ, τ ). Each summand in Eq (25) is the expected inverse mean ISI a neuron at the time that corresponds to x k .

In-vitro data
For validation purposes we considered recorded activity of cortical pyramidal neurons from mouse brain slices [27]. Each recorded neuron was stimulated by a fluctuating current with where I 0 denotes the offset, ∆I s the amplitude of the sinusoidal component with period T , and ∆I n is the magnitude of the noise component N (t), which is an Ornstein-Uhlenbeck process with zero mean, unit variance and correlation time 3 ms. The parameters I 0 , ∆I s and ∆I n were tuned such that the neuronal spike rate oscillated between 2 and 6 spikes/s. We considered recordings that resulted in an average spike rate of at least 4 spikes/s. Oscillation period was T ∈ {4, 8, 16} s. Each recording lasted 68 s: for the first 4 s the input current was constant, for the remaining 64 s it was described by Eq (26). For each of these 64 s segments we fitted a doubly-stochastic I&F model using the OUP variant. We aimed to capture the oscillatory dynamics by the OUP x(t) and approximated the fast fluctuations of the true input (note the small correlation time) by the Gaussian white noise process in the model. We determined the optimal noise amplitude σ in the model by maximizing the likelihood for a range of values of σ in the interval [0.1, 3] mV/ √ ms. To quantify oscillations in the inferred input and compare them with the true signal we first computed the time series of the expected mean input µ(t) = C x(t) +μ using the density p(x k |s 1:K ,θ,τ ) (cf. Methods section 5) and linear interpolation for values of the time series between spike times. Due to strong onset transients in the spiking activity we disregarded the first 8 s of the time series and linearly detrended µ(t) (using the function signal.detrend from the Python package scipy) for these comparisons.
The instantaneous phase of this time series, φ est (t), was then obtained using the Hilbert transform. Analogously, the instantaneous phase of the true signal, φ true (t), was calculated using the Hilbert transform of the sinusoidal input component. To visualize the oscillations in the true and inferred time series we used sin(φ true (t)) and sin(φ est (t)), respectively (Fig 4A-C). We next computed the synchrony index (SI) [52] as where t 1 , . . . , t L are the time points of experimental measurements. The radius (or magnitude) of the SI measures the strength of locking and the angle indicates the average phase shift between the two oscillations.
To measure the strength of locking of neuronal spiking to the (true) oscillatory signal we calculated the resultant vector length (RVL) of the circular mean of phase at spike times t 0 , . . . , t K , 7 In-vivo data We applied our method to extracellular recordings from primary visual cortex of an anaesthetized monkey [28,53]. Extracellular potentials were recorded using an implanted multi-electrode array. After spike sorting, the data consisted of spike times of several single and multi-units. Following [53], we discarded units with signal-to-noise ratio < 2.75 (average wave form amplitude divided by the standard deviation of the wave form noise) as multi-units. We further excluded units with average spike rate > 100 Hz and omitted unphysiologically short ISIs < 3 ms. This resulted in 26 single units that we used for fitting. To account for potential spike sorting errors in this dataset we replaced the factors of the conditional likelihood in Eq (14) by P error u(s k ) + (1 − P error )p(s k |µ k,i k , ϑ k ), where P error is the probability of a spike sorting error and u(s) is the density of a uniform distribution over a reasonably bounded interval. We set P error = 0.05 in consistency with experimental observations [54]. We first fitted a doubly-stochastic I&F neuron for each of these units and determined the optimal value of σ i using the OUP and MJP separately. For population analyses we then considered the 20 single units with the slowest inferred processes (i.e., those with the largest time constant τ , using averages across the results from the OUP and MJP model variants). We continued the fitting procedure for that population as described in Methods section 3.