A Statistical Model for In Vivo Neuronal Dynamics

Single neuron models have a long tradition in computational neuroscience. Detailed biophysical models such as the Hodgkin-Huxley model as well as simplified neuron models such as the class of integrate-and-fire models relate the input current to the membrane potential of the neuron. Those types of models have been extensively fitted to in vitro data where the input current is controlled. Those models are however of little use when it comes to characterize intracellular in vivo recordings since the input to the neuron is not known. Here we propose a novel single neuron model that characterizes the statistical properties of in vivo recordings. More specifically, we propose a stochastic process where the subthreshold membrane potential follows a Gaussian process and the spike emission intensity depends nonlinearly on the membrane potential as well as the spiking history. We first show that the model has a rich dynamical repertoire since it can capture arbitrary subthreshold autocovariance functions, firing-rate adaptations as well as arbitrary shapes of the action potential. We then show that this model can be efficiently fitted to data without overfitting. We finally show that this model can be used to characterize and therefore precisely compare various intracellular in vivo recordings from different animals and experimental conditions.

, i, j = 1, ..., n, where I denotes the imaginary unit. In practice, discrete Fourier transforms are not actually computed by matrix multiplication, but by means of a Fast Fourier Transform (FFT) algorithm.

Circulant matrices
In order to reduce the computational complexity of the likelihood estimation, we approximate the autocovariance matrix K (which is a Toeplitz matrix K) with a circulant matrix C. By definition a circulant matrix can be expressed as C ij = c (i−j mod n)+1 (S2) we write C = C n (c). All circulant matrices of dimension n can be diagonalized by the unitary discrete Fourier transform matrix U = 1 √ n F n : where † is the conjugate transpose. This implies thatĉ is the vector of eigenvalues of C, and it is a vector with real entries. This makes calculation of inverse and determinant of C extremely cheap, as is multiplication of C −1 by a vector x ∈ R n , which simplifies to where the vector in brackets is the component-wise quotient of the vectors U n x and F n c.

Circulant approximation
The task is now to find a circulant matrix which is as close as possible to the covariance matrix K. This can be formalized as the following minimization problem: where N denotes a multivariate Gaussian with specified mean vector and covariance matrix. This problem has the unique solution Proof: The Kullback-Leibler divergence between two Gaussians is given by We have and hence where the constant does not depend on c. We obtain the derivative and therefore, at the stationary point we havê The sum of roots of unity over j only gives a non-zero value if the integer q = l − m + 1 − i is a multiple of n. Since q has a maximum of q = n − 1 when l = n, m = i = 1 and a minimum of q = 2 − 2n when l = 1, m = i = n, only q = 0 and q = −n are eligible. Hence, The second equality is obtained by reparametrizing r = l − m.

Sampling using FFTs
In order generate a sample u of length n from a multivariate Gaussian with mean vector m and circulant covariance matrix C = C n (c), one generates a zero-mean white noise vector x with unit variance and then uses FFTs to compute u, i.e.
as the following calculation shows, the covariance comes out correctly

Optimization method
The optimization scheme used is a quasi-Newton method, where the vector of parameters Θ is updated according to where f is the function to be minimized (e.g. − log p(u som , s)), ∇f is the gradient, and the matrix B is chosen to be Where H f is the Hessian of f , γ (k) denotes a learning rate, and G is a metric tensor on the parameter space which is used to rescale the parameters to lie in similar ranges. The learning rate γ (k) < 0 is increased when the previous step was successful (typically, by 10 percent), and reduced when the optimizer either runs into boundaries of the admissible parameter region or increases the value of the function (we used a reduction by a factor of 2). Below, we derive the formulae for the gradient and Hessian required for the optimization. The derivations hold for the case where the GP covariance function k is parametrized arbitrarily by θ i and the spike-shape kernel and adaptation kernel are given by linear combinations of basis functions α (k) , k = 1, ..., n a and η (k) , k = 1, ..., n w respectively.

Gradient
The likelihood function has the form where q i = ∆te βui+Ai+log r0 (S18) is the probability of a spike in bin i and depends on all parameters except the ones that parametrize the covariance function k. The membrane potential is given implicitly by It is worthwhile to write the derivatives of log p in the following form: Now, let us evaluate all the terms one by one. The Fourier transform of the circulant covariance c only depends on the GP kernel parameters θ, i.e.
Let us turn to the q terms next. Their differential is where by (S19) and using the fact that α is a linear combination of basis kernels nα k=1 a k α (k) Moreover, A i is also a linear combination, so Therefore (S22) can be written as Lastly, by (S23) we also have Using all the previous results can be regrouped to yield

Hessian
For the Hessian, we mainly need the following where The components of the Hessian matrix are computed as follows ∂ 2 log p(u som , s) ∂ 2 log p(u som , s) ∂u 2 The components as given by equations (33-51) are then combined to form the Hessian matrix H f , which is used in Eq. (S16).