Inferring and validating mechanistic models of neural microcircuits based on spike-train data

The interpretation of neuronal spike train recordings often relies on abstract statistical models that allow for principled parameter estimation and model selection but provide only limited insights into underlying microcircuits. In contrast, mechanistic models are useful to interpret microcircuit dynamics, but are rarely quantitatively matched to experimental data due to methodological challenges. Here we present analytical methods to efficiently fit spiking circuit models to single-trial spike trains. Using derived likelihood functions, we statistically infer the mean and variance of hidden inputs, neuronal adaptation properties and connectivity for coupled integrate-and-fire neurons. Comprehensive evaluations on synthetic data, validations using ground truth in-vitro and in-vivo recordings, and comparisons with existing techniques demonstrate that parameter estimation is very accurate and efficient, even for highly subsampled networks. Our methods bridge statistical, data-driven and theoretical, model-based neurosciences at the level of spiking circuits, for the purpose of a quantitative, mechanistic interpretation of recorded neuronal population activity.


Exponential integrate-and-fire neuron model
The exponential I&F model is specified by Eqs. (5) and (6) in Methods section "I&F neuron models" with function f defined by where ∆ T is the threshold slope factor and V T denotes the effective threshold voltage. The exponential term in Eq. (1) effectively approximates the rapidly increasing Na + current at spike initiation [1] and yields an improved fit to intracellular measurements of current-voltage relationships (typically obtained in-vitro) compared to the leaky I&F model [2,3]. The value of V s does not play an important role as long as it is sufficiently larger than the effective threshold V T . In fact, when V increases beyond V T , it diverges to infinity in finite time due to the exponential term, which defines a spike. In practice, however, the spike is said to occur when V reaches V s . The exponential I&F model used here is further equipped with an absolute refractory period of length T ref , during which V (t) is clamped at V r . A change of V T or ∆ T in the exponential I&F model can be completely compensated in terms of spiking dynamics by appropriate changes of µ, σ, V r and V s . This can be seen using the change of variables V := (V − V T )/∆ T . Consequently, the relevant parameters for inference from spike train data are µ, σ, τ m , V r and T ref .

Numerical solution for ISI density p ISI
Below we present two numerical solution schemes: the first one uses a finite volume discretization and is generally applicable for arbitrary variations of the input parameters, the second one uses the Fourier transform and is suitable for small amplitude variations of the parameters. Depending on the model scenario the latter scheme may be computationally more efficient than the former. Python code for both numerical schemes is available at https://github.com/neuromethods/inference-for-integrate-and-fire-models.

Solution based on finite volume discretization
This method employs a recently developed finite volume numerical scheme with implicit time discretization and Scharfetter-Gummel flux approximation, adapted from [4] for the first passage time problem. It provides an accurate solution of Eqs. (9)-(13) in Methods section "Method 1: conditioned spike time likelihood" for arbitrary variations of the mean input. In brief, we first discretize the domain Within each cell the numerical approximation of p V (V, s) is assumed to be constant and corresponds to the average value denoted by p V (V m , s). Integrating Eq. (9) over the volume of cell m, and applying the divergence theorem, yields where ∆V = V 2 − V 1 . To solve Eq. (2) forward in time (represented by the ISI variable s) the fluxes at the borders of each cell are approximated using the first order Scharfetter-Gummel flux [5], where v m+ 1 2 (s) = f (V m+ 1 2 ) + µ ISI (s) and D = 1 2 σ 2 denote the drift and diffusion coefficients, respectively. For the time discretization we rewrite Eq. (2) (with Eq. (3)) in vectorized form and approximate the involved time derivative as first order backward difference to ensure numerical stability. This yields in each time step of length ∆s a linear system for the vector p n+1 of probability density values at s n+1 , given the values p n at the previous time step s n , with vector elements p n m = p V (V m , s n ), where I is the identity matrix and G n ∈ R N V ×N V is a tridiagonal matrix that contains the discretization of the membrane voltage (cf. Eqs. (2), (3)), including the absorbing and reflecting boundary conditions (Eqs. (12) and (13) of Methods section "Method 1: conditioned spike time likelihood"). The ISI density, Eq. (14), in this representation is obtained by evaluating the flux, Eq. (3) above, at the spike voltage V s , taking into account the absorbing boundary condition, Eq. (12), and introducing an auxiliary ghost cell [6] with center V N V +1 , which yields For additional details we refer to [4] (incl. supplemental material therein). Note that this method also allows for the variation of other parameters, in addition to the mean input, within the ISI. Naturally, the finite volume scheme can be used to compute the first order approximation (15) for small amplitude variations of the mean input in a straightforward way: using J = 0 for p 0 ISI , and a small value J > 0 to calculate p ISI , to obtain p 1 ISI ≈ (p ISI − p 0 ISI )/J. An alternative method to compute p 0 ISI and p 1 ISI , which may be computationally more efficient in certain cases, is presented in the next section.

Solution for weak input perturbations based on Fourier transform
The method described below is inspired by [7] and provides a solution of Eqs. (9)-(13) in Methods section "Method 1: conditioned spike time likelihood" for constant mean input, or small amplitude variations of it, through a perturbative approach. It extends the technique proposed in [7] to the case of time-varying mean input within the ISI, and can be adjusted for small amplitude variations of other parameters in a straightforward way.
We consider µ ISI (s) = µ 0 ISI + Jµ 1 ISI (s) with small |J| according to Eq. (16). The solution of the first passage time system to first order in J can be expressed as We insert this in the Fokker-Planck system, Eqs. (9)-(13), and neglect terms of order > 1 in J. The resulting system is then Fourier-transformed over time using separation of variables,

with Fourier transform and the inverse transform defined bŷ
where ω denotes angular frequency. This yields the following two (coupled) systems, one for the steady-state and one forp 1 The ISI density (in the frequency domain) is then given bŷ We solve the linear ordinary differential equations (7)-(13) above for each fixed frequency ω by splitting the solutions,p and analogously forp 1 These solution components can be conveniently computed by backward integration from V s to a sufficiently small lower bound V lb < V r . We then obtainp 0 ISI andp 1 ISI by satisfying the reflecting boundary conditions (in (9) and (13)): . (16) The ISI density in the time domain is finally calculated by the inverse transform, p ISI (s|µ ISI [0, s], θ) = F −1 [p ISI (ω)] using Eqs. (14) and (16). In practice, we first solve the steady-state system (7)-(10), evaluate the functionĥ which appears in Eq. (12) and then solve the system for the change in ISI density due to the perturbation, (11)-(13). Note that knowledge of µ 1 ISI (s) for s ≥ 0 is required to calculateĥ. If the input perturbations are given by delta pulses (as for the networks in Results section "Inference of synaptic coupling") the calculation ofĥ is greatly simplified; e.g., for a pulse at s = s 0 , µ 1

Related I&F-based method for benchmarks
We considered the approach from [8] for comparison. This method infers the parameters of a generalized leaky I&F model with spike history dependence in response to a fluctuating background input and a timevarying stimulus by maximizing the spike train likelihood. For this purpose an approximation based on the Fokker-Planck equation is used [9]. The method is applicable to the model scenarios in Results sections "Inference of background inputs", "Inference of input perturbations" and "Inference of neuronal adaptation". We ensured equal conditions for the benchmarks in terms of model dynamics, parameters to be estimated, amount of observed data and computational resources. An efficient implementation in Matlab with integrated C/C++ code for accelerated program execution was used for the method from [8]; the code is available at http://pillowlab.princeton.edu/code GIF.html. All tests were performed on the same personal computer. Specifically, the model from [8] adapted for our benchmarks involves Eqs. (5) and (7) in Methods section "I&F neuron models", where the time-varying mean input µ(t) consists of three additive components: a constant baseline µ, perturbations given by the convolution of the kernel k(t) with the sequence of perturbation pulses x(t) = l δ(t −t l ), i.e., k * x(t) := ∞ 0 k(τ )x(t − τ )dτ , and spike history-dependent input given by the convolution of kernel h(t) with the spike train s(t) = k δ(t − t k ), i.e., h * s(t). The kernels k(t) and h(t) were parametrically described using a particular set of basis functions (raised cosines with log scaling of the time axis; for details see [8]). To ensure equal numbers of estimated parameters for the different methods we used two basis functions for each kernel, whose weights were the parameters to be inferred. The basis functions were selected to yield, when equally weighted, a kernel shape that is similar to the shape of the respective equivalent kernel in our description. The selected basis functions and resulting kernels are shown in Supplementary Figs. 3e and 8d.
For benchmarks on inference of background inputs (cf. Results section "Inference of background inputs") we omitted the perturbation and spike history components, k(t) = h(t) = 0; for inference of input perturbations (cf. Results section "Inference of input perturbations") we set h(t) = 0 and measured detection sensitivity as for methods 1a and 2; for inference of neuronal adaptation (cf. Results section "Inference of neuronal adaptation") we set k(t) = 0.

Ground truth connection labels from in-vivo data
(Mono)synaptic connections were assumed to produce excess synchrony above baseline co-modulation in a short interval following PYR spikes. To generate estimated ground truth connection labels we adopted the approach in [10]. A CCG (0.4 ms binning) was computed for each PYR-INT pair using the evoked PYR spikes. For positive labels the peak of the CCG in the interval [0.8, 2.8] ms needed to exceed that from the slowly co-modulated baseline, and it needed to be significantly larger than the largest peak in the anticausal direction (short negative lags). The lower frequency baseline was computed by convolving the observed CCG with a partially hollow Gaussian kernel, with a standard deviation of 10 ms and a hollow fraction of 60%. We estimated the probability of obtaining a spike count sc of n as observed (or higher) in the m th time lag within [0.8, 2.8] ms of the CCG, given the expected, low frequency baseline rate λ s (m) using the Poisson distribution with a continuity correction, Similarly, we estimated the probability of obtaining the observed spike count n (or higher) in the m th time lag within [0.8, 2.8] ms of the CCG as expected from the maximum rate λ a across negative lags within [−2, 0] ms using the Poisson distribution with a continuity correction, Connections were labeled as synapses if P fast < 0.001 and P causal < 0.0026 for all (binned) lags in [0.8, 2.8] ms, according to the rigorous criterion defined in [10]. m e a n i n p u t i n p u t s t d .

Supplementary Figures
m e a n i n p u t i n p u t s t d .
m e a n i n p u t i n p u t s t d .

PYR INT PYR INT
Supplementary Figure 2. Correlation between estimated input parameters and empirical input statistics from in-vitro ground truth data. Pearson correlation coefficient between estimated and empirical input statistics for all PYRs and INTs, and for two stimulus durations: 5 s-long simuli (cf. Fig. 2c) and 15 s-long stimuli (corresponding to 3 repetitions of 5 s stimuli with equal statistics). and τ (center) as a function of true J using methods 1a (red) and 2 (blue), and for the two perturbation parameters using the method from [8] (right). Estimates are based on spike trains of length 100 s (for details see Supplementary Methods section 3). True τ = 10 ms. b: central 50 % of relative errors between estimated and true parameter values for µ, σ and τm from our approach (left) and the one from [8] (right) using the same data as in a. Note that the parameters of input perturbations as well as µ, σ and τm were simultaneously inferred. c: detection sensitivity as a function of J for the methods used in a (solid lines) and a CCG-based method (dashed line, cf. Fig. 3d). d: central 50 % of computation times (left) and numbers of observed spikes and events (right) as a function of J, corresponding to the estimations in a-c. Note that slight asymmetries in a-d for excitatory compared to inhibitory perturbations may be due to differences in the number of observed spikes. e, left: true time series of mean input perturbation triggered at time 0 for J = 0.2 mV/ms (green) and reconstructions from the methods used in a (gray; best 75% according to the combined relative errors in a); right: selected basis functions for the method from [8].   Fig. 5g as a function of recording duration for the three methods (I&F, GLM, CCG); average mean absolute error MAE between true and estimated coupling strengths and delays, respectively, for the I&F method (bottom). c-f : detailed inference results from the three methods (as shown in Fig. 5d, 5e right and 5f) corresponding to the networks in Fig. 5h- Figure 6. Estimation results for strong synaptic coupling using simulated data. a,b: inference results from the three methods (as shown in Fig. 5d, 5e right, 5f and 5h-k left) for subsampled networks of Ntot = 1000 (500 excitatory, 500 inhibitory) neurons with connection probability of 0.1 (a) and 0.2 (b); otherwise as in Fig. 5h  . True τw = 100 ms. b: central 50 % of relative errors between estimated and true parameter values for µ and σ from method 1 (left) and the method from [8] (right) using the same data as in a. Note that the adaptation parameters as well as µ and σ were simultaneously inferred. c: central 50 % of computation times as a function of ∆w, corresponding to the estimations in a and b. d, left: true time series of adaptation current triggered by a spike at time 0 for ∆w = 0.5 mV/ms (green) and reconstructions from the methods used in a (gray; best 75% according to the combined relative errors in a); right: selected basis functions for the method from [8].
m e a n i n p u t i n p u t s t d .
PYR INT m e a n i n p u t i n p u t s t d . Supplementary Figure 9. Correlation between estimated input parameters and empirical input statistics for adaptive and nonadaptive models. Pearson correlation coefficient between estimated and empirical input statistics for the I&F model with adaptation (magenta symbols) and without adaptation (black symbols) for all PYRs and INTs.

Supplementary Tables
Supplementary Table 1 < 0.2 < 10 < 2 < 17 < 117 < 7 All computations were performed on a hexa-core personal computer. † Applies to each of the two examples. ‡ Computation time scaled sub-linearly with recording time (hence with number of observed spikes) and supra-linearly with number of observed neurons.