Neurons transform time-varying inputs into action potentials emitted stochastically at a time dependent rate. The mapping from current input to output firing rate is often represented with the help of phenomenological models such as the linear-nonlinear (LN) cascade, in which the output firing rate is estimated by applying to the input successively a linear temporal filter and a static non-linear transformation. These simplified models leave out the biophysical details of action potential generation. It is not a priori clear to which extent the input-output mapping of biophysically more realistic, spiking neuron models can be reduced to a simple linear-nonlinear cascade. Here we investigate this question for the leaky integrate-and-fire (LIF), exponential integrate-and-fire (EIF) and conductance-based Wang-Buzsáki models in presence of background synaptic activity. We exploit available analytic results for these models to determine the corresponding linear filter and static non-linearity in a parameter-free form. We show that the obtained functions are identical to the linear filter and static non-linearity determined using standard reverse correlation analysis. We then quantitatively compare the output of the corresponding linear-nonlinear cascade with numerical simulations of spiking neurons, systematically varying the parameters of input signal and background noise. We find that the LN cascade provides accurate estimates of the firing rates of spiking neurons in most of parameter space. For the EIF and Wang-Buzsáki models, we show that the LN cascade can be reduced to a firing rate model, the timescale of which we determine analytically. Finally we introduce an adaptive timescale rate model in which the timescale of the linear filter depends on the instantaneous firing rate. This model leads to highly accurate estimates of instantaneous firing rates.
Neurons encode stimuli by emitting trains of actions potentials in response to sensory inputs. To uncover the corresponding neural code, the mapping between sensory inputs and output action potentials needs to be described with the help of a quantitative model [1]. In the recent years, generalized linear models have become a popular class of models for that purpose [2]–[4]. The most basic version of these models is the linear-nonlinear (LN) cascade, in which the instantaneous firing rate of the neuron is estimated by applying to the sensory signal successively a linear temporal filter and a static non-linear function. Phenomenological models of that kind are attractive because they are simple and efficient, and yet allow for enough freedom to fit experimental data. A drawback of this approach is however a lack of a direct relationship between the parameters of a LN cascade and the underlying biophysics, and it has been debated to what extent such descriptions capture the temporal dynamics of spike trains of real neurons [5], [6]. In more detailed models of the neural input-output mapping, membrane potential dynamics play the role of the intermediate between input currents and output action potentials [7]. While more biophysically faithful than linear-nonlinear models, these spiking neuron models are also significantly more complex and a significant amount of effort has been invested in reducing the dynamics of populations of spiking neurons to an effective mapping between their input and the output firing rate [8]–[17]. In the firing rate models (see [18] chapter 7. 2 and [7] chapter 6), the input-output mapping of individual units is essentially a linear-nonlinear cascade, the linear filter being usually a simple exponential. Hence the two problems of relating a LN cascade to biophysical parameters and representing dynamics of spiking neurons by a firing rate model are very closely related. In this communication, we examine to what extent a linear-nonlinear cascade can quantitatively reproduce the firing rate dynamics of spiking neuron models. To this end, we exploit known analytic results for integrate-and-fire models to obtain parameter-free expressions for the linear filter and static non-linearity. We then compare quantitatively the estimates of instantaneous firing rates obtained from various LN models with results from simulations of spiking neurons. For both the leaky integrate-and-fire (LIF) and exponential integrate-and-fire (EIF) models, in most of parameters space we find a good match between the estimate and the simulation results. In the case of the EIF, we show that a single exponential provides a good approximation for the linear filter, so that the LN cascade reduces to a firing rate model, the time constant of which we compute analytically. We then introduce an adaptive timescale rate model in which the decay time of the linear filter depends on the instantaneous firing rate, and show that this model provides a significant improvement with respect to both the basic rate model and the LN cascade. Finally, we examine a conductance-based spiking neuron and find that in this case also the adaptive timescale rate model provides an excellent description of firing rate dynamics.
We wish to approximate the trial-averaged firing rate of the neuron using a linear-nonlinear cascade: (1) where is the signal input, is a temporal linear filter, is a static non-linearity, and is the convolution between and. Moreover, , where is a Gaussian process of zero mean, unit variance and correlation time. The LN approximation of firing rate dynamics becomes exact in two extreme cases: (i) the linear limit of vanishing signal amplitude; (ii) the adiabatic limit of very long signal correlation time, . We first determine the linear filter and the static non-linearity in these two limits. To extend the obtained expressions to the full parameter space, we simply set the linear filter and the static non-linearity to be independent of the input parameters and, in which case the linear and adiabatic limits fully determine and. In this section we describe the general approach, which we later apply to specific neural models. We start by examining the leaky integrate-and-fire (LIF) model [28], in which action potentials are generated when the membrane potential crosses a fixed threshold value, the dynamics of the membrane potential being governed only by a leak current. Despite its apparent simplicity, this model is capable of reproducing quantitatively the transfer function of neocortical pyramidal cells in presence of in vivo like noisy inputs [29]. The LIF model has been studied extensively [30]–[32], and a number of analytic results are available for it. We first describe the linear filter and static non-linearity for the LIF model, and then assess the accuracy of the LN approximation for the input-output mapping. The exponential integrate-and-fire (EIF) model [21] is a generalized integrate-and-fire model in which the action potentials are generated by an exponential current instead of a fixed threshold. It is the simplest model capable of reproducing membrane potential dynamics and action potential times of cortical neurons in presence of a noisy input, and it has been used to successfully fit data from layer pyramidal cells [39], interneurons [40] as well as cerebellar Purkinje cells [41]. In this section, we examine an EIF neuron with parameters corresponding to pyramidal cells (see Materials and Methods). So far we examined only models of the integrate-and-fire type, which are one-dimensional in the sense that action potential generation is controlled by a single variable, the membrane potential. In contrast, in biophysically more detailed models, the dynamics of the membrane potential are coupled to the dynamics of a number of ionic conductances, so that these models have higher dimensionality. In spite of this additional complexity, we will show that our results can be easily extended to a standard conductance-based model of Hodgkin-Huxley type, the Wang-Buzsáki model [42]. Studying the dynamics of conductance-based models in the presence of noise is in general very challenging, and the transfer and linear response functions are in general not known analytically. It has however been found that the exponential integrate-and-fire model with appropriately chosen threshold, reset, spike sharpness and refractory period closely reproduces the transfer and linear response functions of the Wang-Buzsáki model [21]. Although the values of these four parameters were chosen so as to reproduce the scaling of the transfer function around threshold and at strong inputs, the linear response functions of the Wang-Buzsáki and EIF models also match for any value of input parameters. In the original study this was observed with the help of direct numerical simulations of both models. We evaluated the transfer function and linear response function of the EIF model using the direct integration of Fokker-Planck equation [26] and comparing with simulations of the Wang-Buzsáki models we confirm the previously observed close match. The linear filter and static non-linearity for the Wang-Buzsáki model can thus be directly obtained from the transfer function and linear response function of the EIF model with appropriate parameters (see Materials and Methods). Fig. 7 illustrates the comparison between the firing rates obtained from the LN approximation, and numerical simulations of the full conductance-based model. As for the EIF model, the match between the LN estimate and the numerical PSTH is good. Moreover, the simplified rate models developed for the EIF model carry over to the Wang-Buzsáki model, and provide simple approximations for the rate dynamics using the transfer function alone. In particular, the adaptive-timescale rate model leads to very high correlation coefficients between the estimated firing rate and the numerical PSTH (cf. Fig. 7). Note that such a high value of can be somewhat misleading: in Fig. 7 A the correlation coefficient of corresponds to a root-mean square distance of Hz between the predicted firing rate and numerical PSTH. This is significantly larger than the error in the PSTH, which is for these parameters of the order of Hz (see Materials and Methods).
Phenomenological firing-rate models (and the closely related neural field models) are basic tools of theoretical neuroscience, and several earlier studies have looked for quantitative mappings between such models and more biophysically detailed, spiking neuron models. To our knowledge, our study is the first to compare extensively across parameter values the output of a phenomenological rate model to the firing rate dynamics of spiking neurons. The question of how to reduce the firing rate dynamics of populations of spiking neurons to simplified ‘firing rate’ models has been the subject of numerous previous studies. Most reductions however ignore the single cell dynamics and eventually end up with rate equations in which the only time scale is a synaptic time scale (see e. g. [45]). Such a reduction can be performed rigorously in the limit in which the dynamics of the synapses are slow [9], [46]. Another approach is to use another slow variable, e. g. an adaptation variable, as the only dynamical variable (see e. g. [15]). The drawback of this type of approach is that one can only capture the dynamics on the slow time scales, and all the fast time scales related to spiking dynamics are lost. Shriki et al (2003) reduced the dynamics of a network of specific conductance-based model neurons to firing rate dynamics, but their approach is based on numerical fits of both the static non-linearity and of the dynamical firing rate response. The correspondence between linear-nonlinear cascade models and spiking neuron models has been examined in several earlier works. In [47], [48], techniques were developed for computing the linear filter and static non-linearity for integrate-and-fire models, while similar questions for the Hodgkin-Huxley model were addressed in [49], [50]. In these works, the authors consider the situation in which background noise is absent, so that the neuron does not fire spontaneously in absence of input signal. In our framework, this corresponds to the double limit of and. The limit of periodically firing neurons, i. e. vanishing noise but non-zero firing rate, was investigated in [51]. In [16], linear-nonlinear cascade models were used to approximate the firing rate dynamics of a spike response model with escape noise [7]. In contrast, here we examined integrate-and-fire and conductance-based models in presence of more biophysically realistic diffusion noise. Similarly to our case, in [7] the linear filter was determined analytically, however the static nonlinearity was obtained by fits to the data. To produce trains of action potentials, the linear-nonlinear cascade model is often supplemented by a third step, a stochastic Poisson process which at every time step generates an action potential with a probability given by the instantaneous firing rate obtained from the cascade. In this study, we have not attempted to compare the full statistics of spike trains produced by such a linear-nonlinear-Poisson model with the statistics of spike trains of integrate-and-fire neurons. Instead we have concentrated on the instantaneous firing rate, which is equivalent to the first-order statistics of spike trains. The instantaneous firing rate provides information about the timing of individual spikes, but does not specify the correlations between successive spikes in a given train. It has been argued that the refractory period and other post-spike effects play an important role in determining precise spike timing [5], [52]–[54]. To reproduce faithfully the full statistics of spike trains of spiking neurons, the linear-nonlinear cascade would have to be supplemented with post-spike history filters leading to correct higher order statistics. Several modeling approaches have been developed to include post-spike filters [54], [55], most prominently generalized linear models (GLMs) [3], [56] and spike-response models (SRMs) [7]. The main difference between these two classes of models is that in SRMs, the quantity obtained after applying the linear filters to the inputs and previous spikes is interpreted as the membrane potential, while no such assumption is made in GLMs. In consequence SRMs are usually fitted to intra-cellular recordings [57], while GLMs are more often applied to extra-cellular recordings [4], [56]. In both classes of models, because of post-spike filters, the firing rate is an implicit function of the input signal, while in conventional LN models as used here the firing rate is an explicit function of the input signal, a very desirable property (for details see [16], [58]). It is not clear how to generalize an LN cascade to take into account correlations in the spike train while preserving this property. It should be noted that the linear filter we use incorporates effects of refractoriness - this is most noticeable at low noise, where the filter exhibits oscillations due to effective refractoriness (see also [55]). While additional effects would need to be incorporated in post-spike filters, these filters should affect only the higher order statistics of spike trains, and not the instantaneous firing rates. A large number of studies have exploited linear-nonlinear models to fit experimentally measured data. In the majority of these studies [2], [4]–[6], [54], the linear-nonlinear model represents the mapping between the stimulus and neuronal firing, and therefore typically encapsulates several processing stages that transform the stimulus into a direct input to the neuron. In contrast, here we considered the mapping between the direct input to the neuron and its output. Such direct mappings have been studied experimentally in vitro [55], [59]–[63]. In these studies, the input-output mapping of cortical neuron was investigated in absence of background noise (note in particular that the spike-triggered average inputs display oscillations as in our low noise case). In vivo recordings indicate that background synaptic noise is a fundamental component of cortical processing [64], [65], as ongoing neural activity in higher cortical areas implies that only a part of the total input to a neuron can controlled by a sensory stimulus. More recent in vitro studies [22], [23], [29], [61] have therefore injected artificial background activity on top of the repeating signal. These studies have however mostly explored the linear regime, and it seems important to further examine the non-linear regime, varying systematically signal and noise parameters. Exponential integrate-and-fire models have been used to predict individual action-potentials of cortical neurons, however post-spike adaptation currents had to be taken into account [39], [66], [67]. We therefore expect that the linear-nonlinear and rate models we developed here for the eIF model will have to be supplemented with additional adaptation components to reproduce accurately the firing rate dynamics of cortical neurons.
In integrate-and-fire models, action potentials are generated solely from the underlying dynamics of the membrane potential [7]. These dynamics are given by [21]: (17) where the membrane potential is determined with respect to the resting potential of the cell, ms is the membrane time constant, is a spike generating current, and is the total current (expressed in mV) elicited by synaptic inputs to the neuron. We studied two different versions of the integrate-and-fire model: Once the membrane potential crosses the threshold, it diverges to infinity in finite time. This divergence represents the firing of an action potential. Following the divergence, the membrane potential is reset to a value after a refractory time. The parameter quantifies the sharpness of the AP initiation. The parameter values used in most of this study were, (a typical value for pyramidal cells [39]), , and ms. We used the Wang-Buzsáki model [42], which is a modified Hodgkin-Huxley model. The dynamics of the membrane potential are given by (19) where is the membrane capacitance (), is the leak current (; mV), is the sodium current, is the delayed rectifier potassium current, and is the total synaptic input current. The activation of the sodium current is assumed instantaneous: (20) while the kinetics of the gating variables and are given by: (21) with. The functions and are given by: (22) (23) (24) The maximum conductance densities and reversal potentials are: , mV; , mV. As explained in the main text, in this study we assume that the synaptic inputs to the neuron are separated into two groups: (i) inputs that are identical across trials, and which we call the “signal” inputs; (ii) inputs that are uncorrelated from trial to trial, which we call the background noise. In consequence, the total synaptic input can be written as (25) We further assume that both signal and noise inputs consists of a sum of large number of synaptic inputs, each individual synaptic input being of small amplitude. We therefore use the diffusion approximation [30], and represent both signal and noise inputs as Gaussian random processes. Within the diffusion approximation, the difference between a conductance-based input and a current-based input is merely a rescaling of the membrane time constant [68], hence here we consider only a current-based input. For convenience, the mean of the input signal is taken to be zero. The correlation time and the standard deviation of are parameters which we systematically vary in this study. Realizations of are generated using (26) where is a Gaussian process, with zero mean, unit variance and vanishing correlation time. The same realization of is used in all trials. The background noise is a Gaussian process of mean, standard deviation and vanishing correlation time, uncorrelated from trial to trial. The parameters and were systematically varied. Here we provide the summary of definitions and expressions for the transfer function and linear response functions of integrate-and-fire neurons. For completeness full derivations are provided below. The transfer function determines the average firing rate of a neuron in response to a steady input of the form, where is a random process of zero mean and unit variance representing background noise, and is the membrane time constant. For the leaky integrate-and-fire neuron receiving background noise uncorrelated in time, the transfer function is given by [27] (27) For the exponential integrate-and-fire neuron receiving background noise uncorrelated in time, the transfer function can be expressed as [21] (28) A convenient method of evaluating Eq. (28) is to integrate the steady state Fokker-Planck equation [26]. The rate response function specifies the trial-averaged firing rate of the neuron in response to a time-varying input of small amplitude [12]. More precisely, in response to an input of the form, with a Gaussian random process independent from trial to trial and identical in all trials, at the linear order the firing rate of the neuron is given by (29) where. Taking the Fourier transform, Eq. (29) becomes (30) where, and are the Fourier transforms of, and. For the LIF receiving a background noise uncorrelated in time, the response function in frequency can be calculated from the Fokker-Planck equation [12], [24]. The result reads: (31) where, , and is given in terms of a combination of hypergeometric functions, or equivalently as the solution of the differential equation (32) with the condition that is bounded as. For the EIF, no explicit expression is available for, but a method has been developed for computing numerically from the Fokker-Planck equation [26]. That procedure is essentially equivalent to integrating Eq. (32) for the LIF model. This is the method we use here. For details, see [26]. To assess the precision of the firing rates predicted by various models, we have systematically compared the predicted firing rates with results of simulations of the LIF, EIF and Wang-Buzsáki neurons. The membrane potential dynamics of the neuronal models were simulated using a standard second-order Runge-Kutta algorithm with a time step of. For each set of parameters, a fixed, s long random instance of the input signal was applied in trials. The output firing rate was then computed by averaging over trials the spike trains binned in windows of ms. To obtain the predicted firing rates, the original input signal was sampled at intervals of ms, convolved with the linear filter (determined with a precision of ms) and/or fed through the static non-linearity. To compare quantitatively the prediction with the numerical firing rate, we computed the Pearson' s correlation coefficient: (33) where is the numerical firing rate and is the firing rate predicted by the model. The average is taken over time bins (34) where is the total number of bins. The value of the Pearson correlation coefficient, which lies between −1 and 1, represents the fraction of variance of accounted for by a linear, instantaneous transformation of. A caveat of the correlation coefficient is that its value can be high even if the means and variances of and are very different: quantifies only the match between the temporal variations of the two time series. We have therefore checked separately that, when values of are high, the LN and rate models provide accurate predictions of the mean and variance of the firing rate (data not shown). An alternative standard measure of the similarity between and is the root mean square distance defined by (35) If the means and variances of the two time series are identical, there is a simple relationship between and: (36) An advantage of the RMS distance over the correlation coefficient is that takes into account the match between the means and variances of and. A disadvantage is that the scale of varies when the input parameters are varied. To compare the predictions of the models as the parameters are varied, we have therefore chosen the dimensionless measure provided by the correlation coefficient. For a fixed set of parameters, the RMS distance is a very useful comparison between the PSTH and various models, as it provides the mean error in units of Hz: we therefore display these values to the graphs comparing the predictions of different models for fixed parameter values (Figs. 3 A, 6 and 7). The value of is bounded from below by the error induced by the finite number of trials used to estimate the PSTH. For a given instantaneous firing rate, the error can be estimated to be where is the number of trials and ms is the size of the time bin. As the instantaneous firing rates vary from to Hz in Figs. 3 A, 6 and 7, the error on the PSTH is of the order of Hz. Note that this precision is far superior to the one that can be reached in experiments, where the number of trials is typically smaller by several orders of magnitude. For the leaky integrate-and-fire neuron receiving background noise uncorrelated in time, the rate response in frequency can be obtained from the first-order perturbation of the steady-state Fokker-Planck equation [24]. The original perturbation study [24] was done in the context of a recurrent inhibitory network, but the rate response response function can be deduced in exactly the same way [12], [25]. Here we provide the direct derivation of the rate response response function, following the same steps as in [24]. We consider a leaky integrate-and-fire neuron with membrane potential dynamics defined by Eq. (17), receiving an input current of the form (37) where is a Gaussian white noise process of zero mean and unit variance, uncorrelated from trial to trial. To study the stochastic dynamics of the membrane potential, we look at the probability distribution of the membrane potential as function of time. The dynamics of the corresponding probability density obey the Fokker-Planck equation [69]: (38) This equation expresses the conservation of probability in time, and can also be written as (39) where the current of probability density is given by (40) The instantaneous firing rate is given by the flux of probability density through the threshold membrane potential: (41) The membrane potential is reset to after crossing the threshold, hence is discontinuous at and obeys: (42) As the membrane potential cannot exceed the threshold, for. Since the probability density current depends on the derivative of with respect to, must be a continuous function of, hence (43) (44) Eqs. (41–44) are the four boundary conditions for the Fokker-Planck Equation. In addition we will require that be integrable as.
Deciphering the encoding of information in the brain implies understanding how individual neurons emit action potentials (APs) in response to time-varying stimuli. This task is made difficult by two facts: (i) although the biophysics of AP generation are well understood, the dynamics of the membrane potential in response to a time-varying input are highly complex; (ii) the firing of APs in response to a given stimulus is inherently stochastic as only a fraction of the inputs to a neuron are directly controlled by the stimulus, the remaining being due to the fluctuating activity of the surrounding network. As a result, the input-output transform of individual neurons is often represented with the help of simplified phenomenological models that do not take into account the biophysical details. In this study, we directly relate a class of such phenomenological models, the so called linear-nonlinear models, with more biophysically detailed spiking neuron models. We provide a quantitative mapping between the two classes of models, and show that the linear-nonlinear models provide a good approximation of the input-output transform of spiking neurons, as long as the fluctuating inputs from the surrounding network are not exceedingly weak.
neuroscience/theoretical neuroscience
From Spiking Neuron Models to Linear-Nonlinear Models