From Spiking Neuron Models to Linear-Nonlinear Models

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.


Introduction
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][3][4]. The most basic version of these models is the linearnonlinear (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][9][10][11][12][13][14][15][16][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 linearnonlinear 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 linearnonlinear 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 parameterfree 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-andfire (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.

Results
We model a typical setup in which a given stimulus is repeatedly applied to a preparation, and action potentials of a neuron are recorded over many trials. We represent this neuron as a spiking neuron (either integrate-and-fire or conductance based) receiving a time-varying input. Here we consider only the case of input current, but our results could be easily extended to an input conductance. This current is assumed to consist of a sum of two components: an input signal, i.e. a time-varying input that is identical in all trials, and a background noise that is uncorrelated from trial to trial. The input signal can be interpreted as a feed-forward input from sensory pathways responding reliably to a stimulus which is identically repeated over trials, while the noise component represents the background activity of the surrounding network and inputs from other areas not directly controlled by the stimulus. Because of the noise, the spiking output varies from trial to trial. We therefore represent the output by its Peri-Stimulus Time Histogram (PSTH), i.e. the time dependent firing rate of the neuron [1], where the instantaneous firing rate is obtained by averaging over trials (see Fig. 1 and Materials and Methods). Equivalently, the PSTH can be interpreted as the firing rate of a large population of uncoupled neurons that all receive an identical input signal as well as background noise that is uncorrelated from neuron to neuron.
Our aim is to examine the extent to which the mapping between the input signal and the output firing rate can be approximated by a linear-nonlinear (LN) cascade consisting of two steps: (i) a linear temporal filter applied to the input signal; (ii) a static non-linear function applied to obtain the instantaneous firing rate (see Fig. 1). A standard method for determining the elements of a LN cascade is to use reverse correlations [19,20], i.e. apply a signal with whitenoise temporal statistics, compute the spike triggered average of the signal, which corresponds to the linear filter, and then determine the associated static non-linear function. Here we use a different approach: we exploit the known analytic results for integrate-and-fire neurons to infer the linear filter and the static non-linearity for particular limits of parameter values. Extending these expressions to the whole parameter space, we obtain the linear filter and static non-linearity in a parameter-free form. We then systematically assess the accuracy of the corresponding LN cascade by comparing its predictions for the firing rate with numerical simulations of spiking neurons.
To limit the available parameter space, we assume that the temporal statistics of the noise input are Gaussian with mean I 0 , standard deviation s and no temporal correlation (white noise), while the temporal statistics of the signal input are Gaussian with zero mean, standard deviation I s (which we refer to as the amplitude of the signal), and correlation time t s (colored noise). Because of background noise, in absence of signal (I s~0 ), the neuron is not silent, but fires at a baseline rate r 0 determined by I 0 and s. In this study, instead of varying I 0 , we vary r 0 , which is equivalent since all model neurons considered in this paper have continuous and monotonically increasing f -I curves.. The parameter space that we explore is thus four-dimensional and consists of r 0 , s, I s and t s .

Determining the elements of the LN cascade
We wish to approximate the trial-averaged firing rate r(t) of the neuron using a linear-nonlinear cascade: where s is the signal input, D is a temporal linear filter, F is a static non-linearity, and D Ã s(t)~Ð ? 0 dtD(t)s(t{t) is the convolution between D and s. Moreover, s(t)~I s n(t), where n(t) is a Gaussian process of zero mean, unit variance and correlation time t s .
The LN approximation of firing rate dynamics becomes exact in two extreme cases: (i) the linear limit of vanishing signal amplitude I s ?0; (ii) the adiabatic limit of very long signal correlation time, t s ??. We first determine the linear filter D and the static non-linearity F in these two limits. To extend the obtained expressions to the full parameter space, we simply set the linear filter D and the static non-linearity F to be independent of the input parameters I s and t s , in which case the linear and adiabatic limits fully determine D and F. In this section we describe the general approach, which we later apply to specific neural models.
Linear filter. For a signal of vanishing amplitude I s , the variation of the firing rate r(t) around the baseline firing rate r 0 is small. Expanding Eq. (1) to linear order, the LN model reduces to where F (0)~r 0 . On the other hand, in the linear limit the firing rate of the spiking neuron is given by where R n (t) is the so-called rate response function of the neuron in presence of white noise [12,[21][22][23]. The Fourier transform of R n (t) can be computed analytically for the leaky integrate-and-fire model [12,24,25], and for the exponential integrate-and-fire model an efficient numerical method has been designed to determine R n (t) from the Fokker-Planck equation [26]. For completeness, the analytic derivation for the LIF is included in the Appendix.

Author Summary
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.
Comparing Eqs. (2) and (3), it is straightforward to identify F '(0)D(t) and R n (t). Note that R n (t) depends on the properties of background noise through r 0 and s.
Static non-linearity. In the limit t s ??, the input signal varies on a timescale that is much longer than the timescale of the linear filter, so that D Ã n(t)~D 0 n(t), where D 0~Ð ? 0 dtD(t), and in the LN approximation In the same limit, as the input signal varies slowly, at each point in time the neuron effectively receives a white noise input of mean I 0 zs(t) and standard deviation s. The response to such an input is given by the so-called f {I or transfer function W : The transfer function for the LIF and EIF models receiving white noise is known analytically [27]. For the EIF an efficient numerical method has been designed to determine W from the Fokker-Planck equation [26]. In both cases, the shape of W depends on the standard deviation s of background noise.
Comparing Eqs. (4) and (5) leads to the following identification: Extension to the full parameter space. For finite signal amplitude I s and correlation time t s , the LN cascade does not provide an exact description of firing rate dynamics, but only an approximation. Here we choose the linear filter D and the static non-linearity F to be independent of the parameters of the signal I s and t s . In that case, the two limits I s ?0 and t s ?? determine D and F uniquely, up to multiplicative constants D 0 and F '(0). From Eq. (6), these two constants are constrained by F '(0)D 0~W '(I 0 ). As the two constants enter the LN model only as a product, without loss of generality we set F '(0)~1, so that D 0~W '(I 0 ). We therefore get: F (L)~W(I 0 zL=W' where W is the transfer function and R n is the rate response function of the neuron. The units of D(t) are Hz/mV, so that D Ã s(t) is the linear estimate of the firing rate in Hertz. Leaky integrate-and-fire neurons 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][31][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.
Linear filter. As we set the linear filter D of our LN cascade to be independent of input parameters, D is given by the rateresponse function of the model neuron. For the LIF model, the rate response function to oscillatory inputs in presence of white noise has been studied in [12,24]. Its analytical expression is known in frequency (see Materials and Methods Eq. (31)), and we use the Fast Fourier Transform to obtain the values in time, as in [33]. The derivation of the rate response function is included in the Appendix. Some analytic results for D(t) have been obtained for the limit t?0 [34,35]. In particular, in that limit, D(t) diverges as t {1=2 , so that the LIF model is capable of responding very quickly to changes in the input signal.
Background noise strongly modulates the response of the neuron [36,37], so that both the timecourse and amplitude of the linear filter D(t) for tw0 depend on the mean I 0 and standard deviation s of the background noise, as shown in Fig. 2.
For strong background noise, D(t) is a monotonically decreasing function, the decay time being of the order of a couple of milliseconds. This decay time depends on the mean I 0 of the background noise, or equivalently the baseline firing rate r 0 of the neuron [33]. For r 0 ?0, D(t) becomes proportional to the subthreshold filter exp ({t=t m ) (this result was not reported in previous works, see Materials and Methods for the derivation), the decay time of D(t) is thus given by the membrane time constant (10 ms in our case). As r 0 is increased, the decay time decreases, so that it is always faster than t m .
For weak background noise, D(t) displays oscillations, at a frequency approximately given by the baseline firing rate r 0 (see Fig. 2). The decay time of D(t) in this case is much longer (of the order of tens to hundreds of milliseconds), and it increases as the standard deviation s of background noise is decreased.
The amplitude of the linear filter approximately scales as the inverse of the standard deviation s of background noise, so that the amplitude of the firing rate modulation in response to an input signal of a given amplitude I s strongly depends on the amplitude of background noise.
It is interesting to compare the analytic linear filter D(t) with linear filters obtained numerically using reverse correlations analysis. In the simulations, we inject an input signal s(t) with very short correlation time (t s~1 ms), and various amplitudes I s , as well as a background noise of mean I 0 and standard deviation s. The numerical linear filter is then obtained by computing the spike-triggered average of s(t) only. As shown in Fig. 2, our analytic filter matches well the numerical STAs, in absence of any fitting parameter. While this was expected for small I s , the match is good for any value of I s . This fact, observed previously in [38], supports our choice of a linear filter independent of I s .
Non-linearity. The analytic derivation of the transfer function W of the LIF in presence of white noise can be found in [27], and the expression is reproduced in Materials and Methods Eq. (27). For weak noise, W exhibits a sharp threshold, above which it is concave. As the standard deviation of noise is increased, the threshold is replaced by an inflexion point.
The static non-linearity F is obtained by rescaling the transfer function W of the neuron, in such a way that the firing rate at the origin is equal to the baseline firing rate r 0 , and the slope is unity (see Eq. (8)). As shown in Fig. 2, depending on whether r 0 lies in the convex or concave region of W, the non-linear function F is either super-linear or sub-linear. Depending on the parameters of background noise, F can thus display qualitatively different behaviors.
To compare our static non-linear function F with the static nonlinearity that would be obtained using a reverse correlation analysis, we follow the standard method [18]: we first compute the linear estimate L(t) of the instantaneous firing rate (in bins of 1 ms); for each bin, we then compare the linear estimate with the actual firing rate obtained from simulations. For a given value of L, there is a distribution of actual firing rates. The numerical nonlinearity is given by the average of this distribution as function of L. The width of the distribution, which can be significant, is not due to statistical noise in the data, but rather to the limitations of the LN model, as this distribution corresponds to variations in the output firing rate that cannot be accounted for by a linearnonlinear transformation of the signal. In Fig. 2, the numerical non-linearities are compared with the non-linear function F obtained from the transfer function, for different values of the signal correlation time t s . While a good match is expected for very large t s , F is close to the reverse-correlation estimates even for small t s . This observation supports our choice of F independent of the signal correlation time t s .
Accuracy of the LN approximation for the rate dynamics. Once the linear filter and the static non-linearity are determined, we are in position to compare the estimate of the instantaneous firing rate provided by the LN cascade with the actual, numerically determined firing rates for different points in parameter space. Fig. 3 A illustrates the comparison between the numerical PSTH and three different approximations: (i) the linear approximation, in which only the linear filter is applied to the signal; (ii) the nonlinear approximation, in which only the static non-linearity F is applied to the signal; (iii) the full LN cascade. The main qualitative differences between the three approximations can be clearly seen in Fig. 3 A. In the linear approximation, the firing rates can take negative values; applying subsequently the nonlinearity corrects for this, and further amplifies or attenuates the firing rates depending on whether the non-linearity is superthreshold or sub-threshold. In the purely non-linear approximation, the firing rate depends only on the instantaneous value of the signal, so that the firing rate fluctuates too fast compared to the PSTH. Applying a linear filter corrects for that. Note that in this example, I s is large enough, and t s short enough, so that the instantaneous firing rate has very large fluctuations (between 0 and 100Hz) on fast time scales. Thus, we are far from the parameter ranges where one would expect the approximations made to derive filter and non-linearity to hold.
The degree to which various approximations match the numerical PSTH clearly depends on the parameters of the input signal and background noise. To get a quantitative comparison, we computed the Pearson's correlation coefficient r between the estimates and the numerical PSTH, as well as the root mean square (RMS) distance d between both, for various values of the parameters (see Methods for pros and cons of the two measures). The correlation coefficient is 1 in case the two traces are identical, and 0 in case the traces are fully uncorrelated, the precise value depending only on the temporal dynamics of the traces and not on their mean and standard deviation. To interpret the measured values of r, we also compute the correlation between the signal and the PSTH, which corresponds to the performance of an estimate of the firing rate obtained by simply rescaling the mean and variance of the signal. For example, in Fig. 3A, the correlation between signal and PSTH is 0.78 (corresponding to d*14Hz) while the correlation between the LN model and the PSTH is 0.92 (d*8Hz). This means that in this example, values of r of order 0.8 correspond to poor approximations of the PSTH, while values of order 0.9 or more represent significant improvements compared to a trivial rescaling of the signal. Fig. 3 displays the performance of different models as the amplitude I s of the signal is increased, the other parameters being held constant. For I s ?0, the linear and LN models become identical. In that limit, the correlation coefficient r between the models and the PSTH approaches 1, but does not strictly reach this limit because of statistical fluctuations in the data: for I s small, the signal-induced variations in the output firing rate are small, and become comparable to the fluctuations in the numerical PSTH that are due to the finite number of trials used for averaging.
As I s is increased, large values of the input signal lead to increasingly precise and reliable emission of action potentials. As a consequence, the variations of the output firing rate become larger, and the non-linearity plays an increasingly important role. The accuracy of the LN model progressively deteriorates, but r remains above 0:9 far into the non-linear regime. As an illustration, Figs. 3 A and C display the comparison between the LN estimates and the numerical PSTHS for two values of input signal amplitude I s , all other parameters being identical. The accuracy of the purely non-linear model is approximately independent of I s , but relatively low (r&0:8). It is only in the limit of very strong signal amplitude that the accuracy of the nonlinear model becomes comparable to the LN model.
The correlation between the signal and the output increases as the correlation time of the signal is increased (Fig. 3). In parallel, the accuracy of the LN model also increases, but to a relatively smaller degree.
Although the non-linear filter and static non-linearity depend on the parameters of background noise (see Fig. 2), the performance of the model varies weakly with the baseline firing rate r 0 and noise amplitude s. The main exception is the limit of small s in which none of the models significantly improves over the correlation between the signal and the output (see Fig. 3 B for an illustration). The reason for this is that for weak noise D(t) displays oscillations at long timescales: the input signal at a given time influences the output rate on a range of different timescales, and the non-linear interferences between these timescales are not well described by a linear-nonlinear cascade. However, the LN approximation performs well already for sw2 mV, the correlation coefficient between the LN estimate and the PSTH being larger than 0:9. This correlation coefficient further increases with increasing s, while it appears to be independent of the baseline firing rate r 0 .
In summary, the linear-nonlinear model of input-output mapping provides a good approximation of the firing rate dynamics for most of the parameter space, two notable exceptions being the limit of weak background noise (s?0) and the limit of very large input signal (I s &1). Illustration for a given set of parameters (r 0~5 Hz, s~6 mV, I s~3 :3 mV, t s~5 ms). From top to bottom: input signal trace; raster of action potentials in a subset of 200 trials; comparison between the numerical PSTH (black) and the linear, LN, and nonlinear estimates. B: Comparison between PSTH and LN estimate at low noise (r 0~3 0 Hz, s~0:5 mV, I s~0 :4 mV, t s~5 mV). C: Comparison between PSTH and LN estimate at very strong input signal (r 0~5 Hz, s~6 mV, I s~6 :5 mV, t s~5 mV). In panels A-C, the values of r and d indicate respectively the correlation coefficiant and root-mean-square distance between the predicted firing rate and the numerical PSTH. D: Correlation coefficient between the numerical PSTH and various estimates as the amplitude I s and correlation time t s of the signal are varied, for fixed background noise parameters (r 0~5 Hz and s~6 mV). From left to right: (i) I s is varied for t s~5 ms; (ii) t s is varied for I s~5 mV; (iii) I s is varied for three increasing values of t s (t s~1 ,5 and 10 ms), shown with curves of increasing darkness. E: Correlation coefficient between the numerical PSTH and various estimates as the baseline firing rate r 0 and the amplitude s of the noise input are varied, for input signal parameters t s~5 ms and I s~0 :5s. From left to right: (i) s is varied for r 0~1 0 Hz; (ii) r 0 is varied for s~6 mV; (iii) data for four increasing values of r 0 (r 0~5 , 15 Exponential integrate-and-fire neurons 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 5 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).
LN cascade. The linear filter D(t) and the static non-linearity F for the EIF model are qualitatively similar to the LIF model, a crucial difference being that for the EIF model D(t) remains finite in the limit t?0, and therefore exhibits slower timescales. The influence of background noise on D(t) is qualitatively comparable to the LIF model: for large background noise D(t) decays monotonically, while for weak background noise it displays oscillations (see Fig. 4). Depending on the parameters of background noise, the static nonlinearity displays the super-linear and sub-linear regimes, as in Fig. 2. Rate model. In the case of the EIF model, as the linear filter D(t) remains finite in the limit t?0, it is possible to exploit the known asymptotic results to derive a simpler analytic approximation for the firing rate dynamics, in which the linear filter is exponential with a single, effective timescale t eff : With such an exponential filter, the linear non-linear cascade of Eq. (1) can be rewritten as so that the dynamics of the firing rate are given by a rate model (see [18] chapter 7.2 and [7] chapter 6).
To derive an analytic expression for the timescale t eff of the equivalent exponential filter, we note that the amplitude of the Fourier transform of D decays as r 0 =(D T 2pt m f ) for large frequencies f , while for vanishing frequencies it reaches a finite value equal to the slope W' of the transfer function W, evaluated at the baseline firing rate r 0 and background noise amplitude s [21]. Matching these asymptotics to those of the Fourier transform A=(1zi2pf t eff ) of an exponential filter, we get: To compare quantitatively the full linear filter D with the approximate exponential filter D eff , in Fig. 4 we display the root mean square distance d between the Fourier transforms of D and D eff , as a function of the baseline firing rates r 0 and noise amplitude s. Independently of r 0 , d displays a minimum as function of s, at a location s min close to 4 mV for small r 0 . Around that minimum, the exponential filter provides an excellent approximation of the full filter. For very small s, the exponential approximation clearly fails to reproduce the oscillations in D(t).
As shown in Fig. 4 C, the value of t eff is a decreasing function of the baseline firing rate r 0 , and a decreasing function of the noise amplitude s (Fig. 4). The effective timescale t eff is faster than the membrane time constant t m (10 ms in our case) in most of parameter space. The only exception is the limit of small s, in which the exponential filter is a bad approximation of the full linear filter, so that values of t eff larger than t m are an artifact of the approximation. Fig. 6 illustrates the comparison between a numerical PSTH and the estimate obtained from the rate model. Fig. 5 displays the accuracy of the rate model approximation, quantified by the correlation coefficient between the estimate and the PSTH, as a function of the parameters of the input signal and background noise. In most of the parameter space, the accuracy of the rate model is comparable to the performance of the LN approximation in which the full filter is used. For some values of the parameters, the rate model appears to outperform the full LN cascade. This happens when the exponential filter is faster than the full linear filter, in which case the rate model predicts large firing rate transients better than the full LN model, but overestimates small transients in contrast to the full LN cascade. When the input signal amplitude is large, large transients dominate in the value of r, so that the rate model leads to larger values of r than the full cascade.
Similarly to the full LN model, the performance of the rate model degrades in the two limits of weak background noise and very large input signal amplitude. The advantage of the rate model over the full LN cascade is its simplicity, which allows for a very efficient and robust implementation.
Adaptive timescale rate model. While the LN and rate models provide good estimates of firing rate dynamics in most of parameter space, their accuracy deteriorates as the amplitude I s of input signal increases. To improve the firing rate estimates in that parameter regime, we introduce adaptive models, in which the linear filter is allowed to change over time.
So far, the linear filter D and static non-linearity F used in various approximations were evaluated at the baseline firing rate r 0 and held constant in time. However, for large input signal amplitude, the instantaneous firing rates deviate strongly from r 0 . As noted previously, the timescale of D decreases with increasing r 0 , hence at high r 0 the neurons respond faster to inputs, but the LN model does not incorporate this effect. To improve the estimates of firing rate at large I s it seems therefore natural to replace at every timestep the static linear filter D(r 0 ) by the linear filter corresponding to the instantaneous firing rate r (i.e. the estimate of the firing rate at the previous time step). We call such a model an adaptive linearnonlinear model as at every timestep the linear filter is matched to the instantaneous firing rate. With such a modification, the firing rate response becomes faster for larger input currents.
In the adaptive LN model, the linear filter has to be computed in principle at every timestep by integrating the Fokker-Planck equation, which is computationally cumbersome. Instead, for the EIF model, at every timestep we approximate the instantaneous filter by the corresponding exponential filter (see Eq. 9). We thus obtain an adaptive timescale rate model: In this model the timescale t eff depends implicitly on time through the instantaneous firing rate, and can be obtained from where r is the instantaneous firing rate (i.e. the estimate of the firing rate at the previous time step). Fig. 6 illustrates the comparison between a numerical PSTH and the estimate obtained from the rate model. When the instantaneous firing rates are high, the adaptive-timescale rate model provides a better estimate than the non-adaptive models. For small firing rates, the adaptive-timescale model overestimates some transients; this effect is due to the fact the the effective timescale t eff at low rates is faster than the timescale of the full linear filter D.
As shown in Fig. 5 B, in contrast to non-adaptive models, the accuracy of the adaptive-timescale rate model does not degrade with increasing input signal amplitude, but remains constant, the correlation coefficient r between the numerical PSTH and the estimate being greater than 0:95. Overall, the adaptive-timescale rate model thus provides an excellent estimate of the firing rate dynamics in most of parameter space. In particular, as function of background noise amplitude s, r displays a broad maximum where it reaches values of 0:99 (cf. Fig. 6). The accuracy of the adaptive-timescale rate model deteriorates in the limit of very small noise amplitude.

Conductance based model
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 conductancebased 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 r can be somewhat misleading: in Fig. 7 A the correlation coefficient of 0:99 corresponds to a root-mean square distance of 5 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 1 Hz (see Materials and Methods).

Discussion
In this study, we examined the ability of phenomenological models to describe the firing rate output of spiking neurons in response to a time-varying input signal that the neurons receive on top of background synaptic noise. The phenomenological models we considered belong to the class of linear-nonlinear cascade models: the firing rate is estimated by first applying a linear filter to the input signal and then correcting for deviations from linearity using a static non-linear function. Instead of using a fitting procedure, the linear filter and static non-linearity were obtained in a parameter-free form by exploiting analytic results valid for particular limits of input signal parameters. This approach allowed us to systematically quantify the accuracy of the phenomenological models by comparing their predictions with results of numerical simulations of spiking neurons. We found that linear non-linear models provide a quantitatively accurate description of firing rate dynamics of leaky integrate-andfire, exponential integrate-and-fire and conductance-based models, as long as the background noise is not excessively weak. In the limit of vanishing variance of background noise, the spiking of neurons exhibits locking to the input signal [43] that cannot be accounted for by the linear-nonlinear cascade. The accuracy of the cascade models also decreases as the amplitude of the input signal is increased, but this effect becomes quantitatively important only for very strong input signals that lead to instantaneous firing rates in the range of hundreds of hertz. As methods for computing analytically the linear filter and static non-linearity are available only for white noise background inputs, we have not systematically studied the situation in which the background noise is colored. However, one could in principle compute numerically both linear filter and static non-linearity, and then use the same approach as introduced here.
For the exponential integrate-and-fire and conductance-based models, the linear filter can be accurately approximated by a single exponential in a large range of noise amplitudes, so that the linearnonlinear model can be reduced to a firing rate model. We obtained a simple analytic expression for the time constant of the rate model, directly relating it to the biophysical parameters of the neuron. The value of the time constant in particular depends on the sharpness of action potential initiation and the baseline firing rate of the neuron.
Interestingly, the EIF model is essentially the only non-linear integrate-and-fire model that can be described by such a simple rate model, since it is the only model in this class whose firing rate response decays as 1=f in the high frequency limit, independently of whether the background noise is white or colored [21]. The firing rate response of LIF neurons decay as 1= ffiffi ffi f p in the case of white noise [12,24], which makes it impossible to reduce the LIF model to a simple rate model. In the case of colored noise, the response stays finite in the high frequency limit [12]. This fact makes it possible to approximate the LIF model with colored noise by an even simpler rate model, in which the firing rate depends instantaneously on the input through the f -I curve (the 'nonlinear' model described here, see also [15,44]). A comparison of the predictions of the nonlinear model with simulations of the LIF neuron with colored noise however showed that the accuracy of the nonlinear model did not improve significantly with respect to the white noise case (data not shown), presumably because the response function of LIF model does not become totally flat as the noise correlation time is increased [44]. Finally, the firing rate response of QIF neurons decays as 1=f 2 . This model could therefore be approximated by a two-variable rate model, in order to reproduce this high frequency behavior.
Finally, we introduced a simple generalization of the rate model in which the time constant depends on the instantaneous firing rate of the neuron. This adaptive-timescale rate model reproduces with a high precision the firing rate dynamics of exponential integrateand-fire and conductance-based models, even for input signals of very large amplitude. This model can be extended to the case of colored noise, and its accuracy degrades gracefully as the correlation time of the background noise is increased (data not shown).

Comparison with previous studies
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 s?0 and r 0 ?0. The limit of periodically firing neurons, i.e. vanishing noise but nonzero 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.

Predicting full spike trains
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 linearnonlinear-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][53][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.

Relationship with experimental data
A large number of studies have exploited linear-nonlinear models to fit experimentally measured data. In the majority of these studies [2,[4][5][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][60][61][62][63]. In these studies, the inputoutput 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.

Integrate-and-fire models
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]: where the membrane potential V is determined with respect to the resting potential of the cell, t m~1 0 ms is the membrane time constant, y(V ) is a spike generating current, and I 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: -Leaky integrate-and-fire (LIF) modelin this model, y(V )~0, there is no spike-generation current, and an action potential (AP) is emitted when the membrane potential crosses a fixed threshold value V T . The membrane potential is subsequently reset to a value V R after a refractory period t rp . The values used in the simulations were V T~2 0 mV, V R~1 0 mV and t rp~2 ms. -Exponential integrate-and-fire (EIF) modelin this model, the spike generation current is exponential: Once the membrane potential crosses the threshold V T , 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 V R after a refractory time t rp . The parameter D T quantifies the sharpness of the AP initiation. The parameter values used in most of this study were D T~1 mV, (a typical value for pyramidal cells [39]), V T~1 0 mV, V R~3 mV and t rp~2 ms.

Conductance based model
We used the Wang-Buzsáki model [42], which is a modified Hodgkin-Huxley model. The dynamics of the membrane potential are given by where c m is the membrane capacitance (c m~1 mF=cm 2 ), I L~gL (V {V L ) is the leak current (g L~0 :1 mS=cm 2 ;V L{ 65 mV), I Na~gNa m 3 h(V {V Na ) is the sodium current, I Kg K n 4 (V {V K ) is the delayed rectifier potassium current, and I(t) is the total synaptic input current.
The activation of the sodium current is assumed instantaneous: while the kinetics of the gating variables h and n are given by: with x~h,n. The functions a x and b x are given by: a n (V )~0 The maximum conductance densities and reversal potentials are: g Na~3 5mS=cm 2 , V Na~5 5 mV; g K~9 mS=cm 2 , V K~{ 90 mV.

Inputs to the neurons
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 I(t) can be written as 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 currentbased 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 I signal (t) is taken to be zero. The correlation time t s and the standard deviation I s of I signal (t) are parameters which we systematically vary in this study. Realizations of I signal (t) are generated using where j (t) is a Gaussian process, with zero mean, unit variance and vanishing correlation time. The same realization of j (t) is used in all trials. The background noise I noise (t) is a Gaussian process of mean I 0 , standard deviation s and vanishing correlation time, uncorrelated from trial to trial. The parameters I 0 and s were systematically varied.

Transfer function and linear response of integrate-andfire neurons
Here we provide the summary of definitions and expressions for the transfer function and linear response functions of integrate-andfire neurons. For completeness full derivations are provided below.
The transfer function W determines the average firing rate of a neuron in response to a steady input of the form I 0 zs ffiffiffiffiffi t m p g(t), where g(t) is a random process of zero mean and unit variance representing background noise, and t m 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] For the exponential integrate-and-fire neuron receiving background noise uncorrelated in time, the transfer function can be expressed as [21] A convenient method of evaluating Eq. (28) is to integrate the steady state Fokker-Planck equation [26].
The rate response function R n (t) 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 I 0 zI 1 (t)zs ffiffiffiffiffi t m p g(t), with g(t) a Gaussian random process independent from trial to trial and I 1 (t) identical in all trials, at the linear order the firing rate of the neuron is given by where r 0~W (I 0 ). Taking the Fourier transform, Eq. (29) becomeŝ wherer r,Î I 1 andR R n are the Fourier transforms of r, I 1 and R n . For the LIF receiving a background noise uncorrelated in time, the response function in frequency R n can be calculated from the Fokker-Planck equation [12,24]. The result reads: where y T~( V T {I 0 )=s, y R~( V r {I 0 )=s, and u(y,v) is given in terms of a combination of hypergeometric functions, or equivalently as the solution of the differential equation with the condition that u is bounded as y?{?.
For the EIF, no explicit expression is available forR R n , but a method has been developed for computingR R n 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].

Comparison between model predictions and simulations
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 10ms. For each set of parameters, a fixed, 5 s long random instance of the input signal was applied in 50,000 trials. The output firing rate was then computed by averaging over trials the spike trains binned in windows of 1 ms.
To obtain the predicted firing rates, the original input signal was sampled at intervals of 1 ms, convolved with the linear filter (determined with a precision of 1 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: where r(t) is the numerical firing rate andr r is the firing rate predicted by the model. The average is taken over time bins where T is the total number of bins. The value of the Pearson correlation coefficient r(r,r r), which lies between 21 and 1, represents the fraction of variance of r accounted for by a linear, instantaneous transformation ofr r. A caveat of the correlation coefficient is that its value can be high even if the means and variances of r(t) andr r(t) are very different: r quantifies only the match between the temporal variations of the two time series. We have therefore checked separately that, when values of r 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 r(t) andr r(t) is the root mean square distance d defined by If the means and variances of the two time series are identical, there is a simple relationship between d(r,r r) and r(r,r r): An advantage of the RMS distance d(r,r r) over the correlation coefficient r(r,r r) is that d takes into account the match between the means and variances of r(t) andr r(t). A disadvantage is that the scale of d(r,r r) 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 r.
For a fixed set of parameters, the RMS distance d 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 d 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 r, the error can be estimated to be ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi n=(Ndt) p where N~50000 is the number of trials and dt~1 ms is the size of the time bin. As the instantaneous firing rates vary from 0 to 100{200 Hz in Figs. 3 A, 6 and 7, the error on the PSTH is of the order of 1{2 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.
Calculating the rate response function of the leaky integrate-and-fire neuron For the leaky integrate-and-fire neuron receiving background noise uncorrelated in time, the rate response in frequencyR R n 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 where g(t) 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 V as function of time. The dynamics of the corresponding probability density P(V,t) obey the Fokker-Planck equation [69]: This equation expresses the conservation of probability in time, and can also be written as where the current of probability density J is given by The instantaneous firing rate is given by the flux of probability density through the threshold membrane potential V T : The membrane potential is reset to V R after crossing the threshold, hence J is discontinuous at V R and obeys: As the membrane potential cannot exceed the threshold, for V wV T P(V ,t)~0. Since the probability density current J depends on the derivative of P with respect to V, P must be a continuous function of V, hence P(V T ,t)~0, Vt, ð43Þ Eqs. (41)(42)(43)(44) are the four boundary conditions for the Fokker-Planck Equation.
In addition we will require that P(V ,t) be integrable as V ?{?.
Steady state. If I 1 (t)~0, the neuron receives a steady input current, and its output rate is constant in time, r(t)~r 0 . Eqs. (41) and (42), then imply that J(V ,t)~r 0 for V R vV vV T and J(V,t)~0 for V vV R . From Eq. (40) the steady-state probability density P 0 (V ) is proportional to r 0 . Integrating Eq. (40) determines P 0 (V ) up to a multiplicative constant, which is obtained from the normalization condition for P. The final result reads To calculate the linear perturbation of the firing rate arising from a time-varying input current I 1 (t) we write r(t)~r 0 (1zr 1 (t)) ð48Þ Q(y,t)~Q 0 (y)zQ 1 (y,t): Keeping only first-order terms, the Fokker-Planck equation becomes where Lf~d dy yf z d 2 dy 2 f , and the boundary conditions read wherer r 1 is the Fourier transform of r 1 .
The solution of Eq. (53) can be expressed aŝ Q Q 1 (y,v)~a z 1 (v)w 1 (y,v)zb z 1 (v)w 2 (y,v)zQ Q p 1 (y,v) ywy r a { 1 (v)w 1 (y,v)zb { 1 (v)w 2 (y,v)zQ Q p 1 (y,v) yvy r ( ð56Þ whereQ Q p 1 is a particular solution of the non-homogeneous equation and w 1 and w 2 are two independent solutions of the homogeneous equation.
As shown in [24], a particular solution is given bŷ The homogeneous equation reads Two independent solutions of Eq. (58) can be expressed as [24] w 1 (y,v)~M½(1{vt m )=2,1=2,{y 2 ð 59Þ w 2 (y,v)~ffi where M½a,b,c are confluent hypergeometric functions [70]. The full solution forQ Q 1 is obtained by determining the four unknown coefficients a { 1 , a z 1 , b { 1 and b z 1 from the four boundary conditions. The boundary conditions however depend onr r 1 , which is not known at this point. This discrepancy is resolved by noting that w 2 decays exponentially as y?{? while w 1 decays algebraically. Hence only w 2 is integrable, and ifQ Q 1 is to be integrable we must have a { 1~0 . Enforcing this condition leads tô where u(y,v)~exp y 2 w 2 (y,v): ð62Þ In conclusion, we havê Note that the function u is a solution of where L Ã f~{y d dy f z d 2 dy 2 f , i.e. L Ã is the operator adjoint to L. In practice, when evaluating Eq. (61), it is often preferable to evaluate u by integrating Eq. (64) rather than using available implementations of the confluent hyper-geometric functions.
Limit of vanishing firing rate. We consider here the limit I 0 ?{?, keeping s fixed. In that limit y th ?? and r 0 ?0. In the limit y??, we have ( [70], Eq.13.