Frontiers in Computational Neuroscience

Once the impulse response h 0 (t) is known, one can predict the output rate r(t) for any arbitrary, but sufficiently small, stimulus perturbation ∆x(t) by a convolution, 4 r t x h t r x h t () ()() ()(). = = + * * 0 0 0 ∆ (1) The output rate r 0 in response to a constant input x 0 is typically described by a nonlinear activation function r 0 = g(x 0) (" f-I curve "). Note, that in several studies the activation function is called " transfer function " (e.g., Destexhe and Sejnowski, 2009). Following the terminology of linear-systems theory, we use the term " transfer function " exclusively for the Fourier transform of the impulse response. The dynamical response properties of neuron populations, i.e., the properties of the transfer function H 0 (f), have been extensively studied both theoretically (e.

Once the impulse response h 0 (t) is known, one can predict the output rate r(t) for any arbitrary, but sufficiently small, stimulus perturbation ∆x(t) by a convolution, 4 r t x h t r x h t ( ) ( )( ) ( )( ). = = + * * The output rate r 0 in response to a constant input x 0 is typically described by a nonlinear activation function r 0 = g(x 0 ) ("f-I curve"). Note, that in several studies the activation function is called "transfer function" (e.g., Destexhe and Sejnowski, 2009). Following the terminology of linear-systems theory, we use the term "transfer function" exclusively for the Fourier transform of the impulse response.

IntroductIon
Instantaneous firing rates are commonly used to describe either the compound spiking activity of neuron ensembles (population rate) or the trial-averaged response of individual neurons to multiple repetitions of the same stimulus. The dynamics of firing rates is often studied by means of population or firing-rate models. The main purpose of such models is to reduce the dimensionality and complexity of the microscopic dynamics of neural networks in order to obtain tools which allow mathematical treatment, efficient simulation, and intuitive understanding. Since the seminal studies of Wilson and Cowan (1972) and Knight (1972a), a zoo of population models has emerged in the literature. Often, they have to be considered purely phenomenological descriptions. Under simplifying assumptions, they can be mathematically derived from the single-neuron dynamics. In many cases, however, the mathematics of these "reduced" models is still very complex.
A standard solution to this problem is to study the dynamics in a restricted region of the state space. By means of linear-response analysis, for example, one can describe the time-dependent output firing rate r(t) of a population of neurons (or an ensemble of trials) in response to some input signal x(t) in the vicinity of a stationary working point (x 0 , r 0 ) for small perturbations ∆x(t) = x(t) -x 0 . The outcome of this analysis is the system's impulse response 1 h 0 (t), or, in the frequency domain, the transfer function H 0 (f) =  [h 0 (t)](f ) 2,3 .
Rate dynamics of leaky integrate-and-fire neurons with strong synapses Institute of Mathematical Sciences and Technology, Norwegian University of Life Sciences, Ås, Norway Firing-rate models provide a practical tool for studying the dynamics of trial-or population-averaged neuronal signals. A wealth of theoretical and experimental studies has been dedicated to the derivation or extraction of such models by investigating the firing-rate response characteristics of ensembles of neurons. The majority of these studies assumes that neurons receive input spikes at a high rate through weak synapses (diffusion approximation). For many biological neural systems, however, this assumption cannot be justified. So far, it is unclear how time-varying presynaptic firing rates are transmitted by a population of neurons if the diffusion assumption is dropped. Here, we numerically investigate the stationary and non-stationary firing-rate response properties of leaky integrate-and-fire neurons receiving input spikes through excitatory synapses with alpha-function shaped postsynaptic currents for strong synaptic weights. Input spike trains are modeled by inhomogeneous Poisson point processes with sinusoidal rate. Average rates, modulation amplitudes, and phases of the period-averaged spike responses are measured for a broad range of stimulus, synapse, and neuron parameters. Across wide parameter regions, the resulting transfer functions can be approximated by a linear first-order low-pass filter. Below a critical synaptic weight, the cutoff frequencies are approximately constant and determined by the synaptic time constants. Only for synapses with unrealistically strong weights are the cutoff frequencies significantly increased. To account for stimuli with larger modulation depths, we combine the measured linear transfer function with the nonlinear response characteristics obtained for stationary inputs. The resulting linear-nonlinear model accurately predicts the population response for a variety of non-sinusoidal stimuli. et al., 2003;Naundorf et al., 2005;Richardson, 2007;Pressley and Troyer, 2009) and experimentally (Knight, 1972b;Silberberg et al., 2004;Köndgen et al., 2008;Blomquist et al., 2009;Boucsein et al., 2009). One of the major insights has been that populations of neurons can track time-varying input signals much faster than single cells. Theoretical studies have shown that the shape of the transfer function H 0 (f ) depends crucially on the properties of the singleneuron dynamics (Knight, 1972a;Fourcaud-Trocmé et al., 2003;Naundorf et al., 2005;Richardson, 2007) and on the filter properties of the synapses (Brunel et al., 2001).

Eilen Nordlie, Tom Tetzlaff* and Gaute T. Einevoll
Most previous studies have investigated the population rate in response to a time-varying noisy input current injected into the soma of the cells. At each instance in time t, the current amplitudes for the ensemble of neurons (trials) were assumed to be uncorrelated and to follow a normal distribution with mean μ(t) and variance σ 2 (t). Strictly speaking, this corresponds to a scenario where spikes from presynaptic sources arrive at the target cells independently with infinite rate, and where the impact of each spike on the postsynaptic neuron, i.e., the synaptic weight, is infinitesimal ("diffusion approximation"; Johannesma, 1968). For many neural systems, these assumptions are violated: In the mammalian neocortex, for example, a considerable fraction of excitatory synapses gives rise to postsynaptic potentials (PSPs) which cannot be considered small (PSP amplitudes larger than 1 mV; see, e.g., Fetz et al., 1991;Song et al., 2005;Bannister and Thomson, 2007;Lefort et al., 2009). The same holds for hippocampal neurons (Fetz et al., 1991). In the early visual pathway, synapses between retinal ganglion cells and cells in the lateral geniculate nucleus (LGN) are so strong that single retinal spikes can initiate action potentials in the thalamic targets (Cleland et al., 1971;Sirovich, 2008). Moreover, the assumption that the rate of incoming events is large is often questionable. Intracellular recordings in behaving rats, for example, have recently shown that cortical neurons may typically fire at much smaller rates (1 s -1 ) than previously thought (Lee et al., 2006). Even though such neurons receive inputs from several thousands of sources, the compound arrival rate of presynaptic events may therefore be relatively small. Finally, inputs to different neurons can be highly correlated, e.g., due to overlapping presynaptic populations (Bruno and Sakmann, 2006;DeWeese and Zador, 2006). Finite synaptic weights, low firing rates and correlated firing cause non-Gaussian input distributions (DeWeese and Zador, 2006). Under these conditions, the results obtained under the diffusion approximation have to be questioned.
In (1), the precise meaning of the "input" x(t) has so far been left open. In previous studies, x(t) has been interpreted either as the mean μ(t) or the variance σ 2 (t) of the input distribution. Depending on this choice, the meaning and properties of the transfer function H 0 (f ) can be very different. It has been shown, for example, that a modulation in the variance of a noisy synaptic input current (multiplicative noise) can be tracked much faster than a modulation in the mean input current (additive noise; Lindner and Schimansky-Geier, 2001;Silberberg et al., 2004;Richardson, 2007;Boucsein et al., 2009;Pressley and Troyer, 2009). In general, a modulation in the firing rate of a presynaptic neuron population affects both the mean and the variance of the ensemble of input currents. If presynaptic spike trains can be treated as realizations of independent Poisson point processes, the input variance is proportional to the input mean. For the perfect integrate-and-fire neuron model, Pressley and Troyer (2009) recently showed that the resulting combination of additive and multiplicative modulation of the input current enables the cell population to perfectly replicate the input signal, i.e., the time-varying input rate. So far it is unclear to what extent this can be generalized to other neuron models or to real neurons.
In the present study, we investigate the stationary and transient response properties, i.e., g(x 0 ) and H 0 (f ), of the leaky integrateand-fire (LIF) neuron in a regime beyond the diffusion limit. Here, neurons receives input spikes at relatively low rates through strong excitatory synapses. In contrast to former studies, the input signal is represented by the time-varying firing rate of a hypothetical presynaptic neuron or neuron population. Thus, we do not impose any assumption on the distribution of input currents. Moreover, with this setup we do not artificially separate the effects resulting from a timevarying input mean and a time-varying input variance. Instead, both parameters are modulated simultaneously. The transfer function is measured by injecting spike trains generated by an inhomogeneous Poisson process with sinusoidal rate profile for a broad range of modulation frequencies and modulation depths. By studying the frequency spectrum of the rate response, we investigate in which regime the cell (or trial) ensemble can be treated as a linear system. To account for nonlinearities, we combine the transfer function with the activation function measured for stationary input. We demonstrate that the resulting linear-nonlinear model can accurately predict the population-averaged spike response of an ensemble of independent LIF neurons for a variety of non-sinusoidal stimuli.

Model descrIptIon
We investigate the firing-rate response properties of the LIF neuron (Lapicque, 1907;Tuckwell, 1988) with membrane capacitance C m , membrane time constant τ m , spike threshold V th , refractory period τ ref , and reset potential V reset (see Table 1). The neuron receives inputs through a linear synapse modeled by an alpha-function kernel (postsynaptic current, PSC)  i.e., by first drawing spike times s i from a homogeneous Poisson process with unit rate and then, in a second step, applying the inverse t i = Λ -1 (s i ) of the cumulative rate function Λ( ) ( ) t as t = ∫ 0 ds to each event time s i of the realization (see e.g., Devroye, 1986;Nawrot et al., 2008;Cardanobile and Rotter, 2010).
A detailed description of the neuron, synapse and input model can be found in Table 1 (Nordlie et al., 2009). Simulations were performed using the exact-integration method presented in (Rotter and Diesmann, 1999) and implemented in NEST (Gewaltig and Diesmann, 2007) with a time resolution of dt = 0.1 ms. The simulation time was set to T = 2 20 dt ≈ 105 s. The parameters V th , V reset , τ ref , C m, and δ were kept constant throughout the whole study. Unless stated otherwise (Sections 3.4 and 3.5), default values τ s = 2 and τ m =10 ms were used for the synaptic and membrane time constants, respectively. Response spike trains were recorded for a wide range of synaptic weights w and stimulus parameters a 0 , a 1 , f stim . The parameter settings are summarized in Table 3.

characterIzatIon of response propertIes
For each parameter set, we constructed a discretized version of the instantaneous firing rate r(t) (spike density; t ∈{k dt|k = 0,1,2,…,T/ dt}) by counting the number of response spikes in each time interval [t,t + dt) and normalized it by dt = 0.1 ms. As dt < τ ref , the spike count in each bin was either 0 or 1. For each stimulus frequency f stim , we obtained the time averaged firing rate r 0 = T −1 R(0) and the response amplitudes For the principal harmonics (k = 1), we also measured the phase φ of R(f) at frequency f stim . The discrete Fourier transform was computed using the Fast Fourier Transform (FFT; Press et al., 1992) implemented in Python/Numpy 5 . To improve the efficiency of the FFT, the simulation time T was chosen such that the number of time steps T/dt = 2 20 is an integer power of 2. Further, to ensure that the stimulus frequencies f stim were represented in the discrete Fourier transforms, they were chosen such that T always contained an integer number of stimulus periods, i.e., T f stim ∈.
Assuming Poissonian spike-train statistics, the variance ∆ 2 r = 4r 0 /T of the amplitude r 1 is determined by the expected firing rate r 0 and the observation time T (Papoulis and Pillai, 2002). This holds for both homogeneous and inhomogeneous Poisson processes and for all frequencies f. To evaluate whether a measured modulation amplitude r 1 significantly differs from zero (homogeneous case), we computed the z-score z r r r T r = = for each amplitude r 1 . For a homogeneous Poisson process, the probability of observing an amplitude r 1 with z > 2 is less than 5% (values above the horizontal dashed lines in Figures 4A,C,E).
The stationary response to constant input (f stim = 0, a 1 = 0) is described by the activation function r 0 = g(a 0 ), measured on a grid a 0 ∈ {0,2,…,100} s −1 (see Figure 3). For the prediction of population-averaged responses (see Section 3.5), a continuous representation of the activation function was obtained by linear interpolation. To characterize the response to nonstationary input (f stim > 0, a 1 > 0), we computed the complex transfer function (Rotter and Diesmann, 1999). For a given synaptic weight w, the PSP maximum psp max (w) can be calculated numerically. We refer to the minimal weight w crit which causes the PSP to reach the threshold, i.e., for which psp max (w crit ) = V th , as the critical weight. The dependence of w crit on the synaptic and membrane time constants τ s and τ m , respectively, is shown in Figure 2. In this work, all synaptic weights are considered relative to w crit , i.e., w r = w/w crit . The present study is restricted to purely excitatory synapses with w > 0.
The input spike train to the cell is modeled as a realization of an inhomogeneous Poisson point process with sinusoidal rate a(t) with oscillation frequency f stim ∈[1,1000] Hz, mean rate a 0 ∈[0,100] s −1 and modulation depth a 1 ∈[0,100] s −1 (see Figure 1). For data analysis, only data for a 1 < a 0 were taken into account (non-rectified input). To test the predictive power of the rate model we also applied other, non-sinusoidal, stimuli, e.g., firing-rate steps (see Section 3.5). Realizations of inhomogeneous Poisson processes with rate function a(t) were generated by applying the timerescaling theorem (see Brown et al., 2001, and references therein),

Model overvIew
-Reset + refractoriness: Exact integration (Rotter and Diesmann, 1999) Critical weight w crit = min. weight w required for with a time constant τ = 1/(2π f c ). As a measure of the robustness of the fit, we computed the standard deviation of fit parameters f c , γ, and d (error bars in Figure 5) obtained from an ensemble (100 trials) of surrogate data, generated by drawing random numbers from a complex Gaussian distribution with mean  H f 0 ( ) stim and standard deviation σ(  for each stimulus frequency f stim . For further analysis, only parameter sets {w r , a 0 , a 1 , τ s ,τ m } with z-scores z(f stim ) ≥ 2 for all f stim ≤ 2f c were taken into account.
Across wide parameter regions, we found that (6) provides a reasonable fit of the measured amplitudes and phases. We do not imply, however, that the first-order low-pass filter is the optimal description. We mainly use it as a tool to extract transfer features like the cutoff frequency and the low-frequency gain, which in turn can be used to construct a firing-rate model (see Section 2.3). In Section 3.5, we demonstrate that such a model can be used to reliably predict population responses to non-sinusoidal stimuli.
A summary of the data analysis can be found in Table 2.

lInear-nonlInear Model
The linear description (1) is frequently extended by combining a state-independent linear kernel h(t) with a nonlinear activation function g(·), resulting in a linear-nonlinear integral-equation model An alternative approach based on Wiener kernels is presented, for example, in (Sakai, 1992). Firing-rate models are often expressed in terms of differential equations (e.g., Wilson and Cowan, 1972). Indeed, for an exponential kernel h(t) (which corresponds to the first-order low-pass filter (6)) and with u(t) = (a *h)(t) and from the measured response amplitude r 1 (f stim ) and phase φ(f stim ) (see Figure 4). To obtain a more compact description of the response characteristics, we fitted the transfer function of a first-order low-pass filter with cutoff frequency f c , low-frequency gain γ and constant (frequency independent) delay d to the measured transfer function H 0 ( f ) by minimizing the Euclidean distance err = Re Im Im FIgure 1 | sketch of the setup. A leaky integrate-and-fire (LIF) neuron receives spike trains generated by a Poisson point process with sinusoidal rate of mean a 0 and modulation amplitude a 1 . The synapse is modeled as a current kernel (PSC) with amplitude w. The response firing rate is characterized by its mean r 0 , amplitude r 1 and phase φ.
Average response rate (first and second harmonics r 1 , r 2 )

Model and sIMulatIon paraMeters
(see Nordbø et al., 2007, and references therein) 6 . Note, that from this point of view the dynamical variable u(t) corresponds to the linear response of the firing rate. It describes neither the average synaptic input nor the membrane potential. A priori, the kernel h(t) can therefore not be interpreted as the synaptic or membrane filter. For the mathematical and numerical analysis, the differential form (10) of the rate model is often more convenient. The integral form (9) represents a more general description for arbitrary linear kernels h(t).
For small perturbations ∆a(t) of a stationary working point a 0 , (9) should be to leading order identical to (1) (for x(t) = a(t)). As a consequence, the dependence of the transfer kernel h 0 (t) on the choice of the working point a 0 is determined by the slope g′(a 0 ) = dg(a 0 )/da 0 of the activation function at this point, i.e., As we will show in Section 3.2, for a wide range of parameters the dependence of the filter kernel h 0 (t) on the working point a 0 is indeed given by the slope g′(a 0 ) of the activation function alone, i.e., the extracted kernel h(t) = h 0 (t)/g′(a 0 ) is independent of the applied stimulus. We therefore pool our results in this region of the parameter space by averaging h(t) over all a 0 and all a 1 < a 0 . Based on the kernel h(t) obtained this way and the measured activation function g(a 0 ), we set up a linear-nonlinear model (9) which allows us to compute the responses r(t) to arbitrary stimuli a(t) (see Section 3.5).
The model (9) was extracted from the Fourier transform of the spike response of a single neuron. This corresponds to an analysis of the period-averaged signal. Provided the responses in consecutive stimulus periods are independent, period averaging is equivalent to averaging across independent trials or across a population of independent neurons. In Section 3.5, we test to what extent the rate model extracted from a single neuron by period averaging can predict the population-averaged spike response of an ensemble of LIF neurons receiving input spike trains with non-sinusoidal firing rate.

statIonary response
In the present study, the synaptic input is purely excitatory. For stationary inputs, the output rate r 0 therefore increases monotonously with the input rate a 0 ( Figure 3A). For subcritical synaptic weights w r < 1, the postsynaptic neuron has to integrate a number of presynaptic events to reach the threshold (cf. Usrey et al., 1998, in the context of the retina-LGN system). In consequence, the output rate r 0 is smaller than the input rate a 0 . As with increasing input 6 In fact, Eq. (9) can be transformed to a set of differential equations for the whole family of exponential kernels rate a 0 more and more charge is accumulated, the transfer ratio r 0 /a 0 increases with a 0 (Figure 3B). The activation function r 0 = g(a 0 ) is therefore concave in this regime. At the critical weight w r = 1, an isolated input spike causes a depolarization of the postsynaptic cell which just reaches the threshold V th . For w r ≥ 1, one could therefore expect that each input spike is transferred, such that r 0 /a 0 = 1. Due to refractoriness, however, a fraction of input spikes is lost. This fraction increases with the input rate a 0 because the number of short inter-spike intervals increases. In consequence, the activation function becomes convex for larger a 0 , i.e., the transfer ratio r 0 /a 0 decreases with a 0 . For strong weights w r > 1, fewer input spikes are lost because the postsynaptic currents are extended in time (alpha-PSCs). Even if an input spike arrives within the refractory period, the resulting PSC can still be suprathreshold after the neuron has recovered. For the same reason, the transfer ratio can even become larger than one for larger synaptic time constants t s (see for example Figures 8L,O). Figure 3B, the transfer ratio r 0 /a 0 is non-continuous at the critical weight w r = 1. This is most apparent for small and moderate input rates a 0 . Even if the synaptic weight is only slightly smaller than the critical weight, the transfer ratio is considerably smaller than 1. At the critical weight, it jumps to values close to 1.

nonstatIonary response
As illustrated in Figure 4, the response to a time-varying input rate can be well described by a linear first-order low-pass filter (6) with cutoff frequency f c , low-frequency gain γ and fixed delay d. Both the cutoff frequencies f c and the low-frequency gain γ are larger for supercritical synaptic weights w > w crit as compared to subcritical ones. Note that the cutoff frequency for an intermediate weight w r = 0.4 (A) is almost identical to the cutoff frequency for a weight w r = 0.95 which is close to the critical weight (C). At the critical weight, the cutoff frequency jumps to higher values (E, w r = 1.05). The dependence of f c and γ on the synaptic weight is explicitly shown in Figures 5A,B (for a 0 = 40 and a 1 = 30 s −1 ). For subcritical weights w < w crit , the rate-transfer cutoff frequency f c is smaller than or close to the cutoff frequency  f s ≈ 50 Hz of the synaptic filter (alpha-kernel with time constant τ s = 2 ms; see Section 2.1). At the critical weight w r = 1, f c jumps to high values above 100 Hz. In contrast to this sharp transition in the cutoff frequency, the low-frequency gain γ increases continuously with w r . For large weights, it approaches values close to one. Here, the response modulation at small frequencies is almost as large as the modulation in the input rate. In agreement with (11), the low-frequency gain γ is well explained by the slope g′(a 0 ) of the activation function at the working point a 0 (thin solid line in Figure 5B).
Figures 5C-F depict the dependence of the filter characteristics f c and γ on the mean a 0 and modulation depth a 1 of the input signal. For subcritical weights w r < 1, the cutoff frequency f c is approximately independent of both a 0 and a 1 (Figures 5C,E). Only for small synaptic time constants τ s , the cutoff frequency f c slightly increases with a 0 (data not shown). For supercritical weights w r > 1, we observe a decrease in f c with a 0 (light gray curve in Figure 5C). The low-frequency gain γ (Figures 5D,F) is in good approximation determined by the slope g′(a 0 ) of the activation function evaluated at the input rate a 0 (thin solid curves in Figures 5D,F).
For a linear-nonlinear model (9), the kernel h(t) is given by h 0 (t)/g′(a 0 ) (see Section 2.3). For subcritical weights, we observe that h(t) is approximately independent of both a 0 and a 1 . An important prerequisite for a model which does not depend on the input a(t) is thereby fulfilled. For the construction of the firing-rate model (9), we therefore pooled the data by averaging h(t) over all a 0 and all a 1 < a 0 . Although the cutoff frequency depends on the choice of

lInearIty of the response
Strictly speaking, measuring the linear response H 0 (f) of a nonlinear system at a working point a 0 requires infinitesimally small the working point a 0 for supercritical weights, we also attempt the same strategy here. In Section 3.5, we will test the performance of this model for non-sinusoidal stimuli. frequency f stim will give rise to a single peak in the response amplitude spectrum at the same frequency f stim . Leaving the linear regime is in general accompanied by the occurrence of higher harmonics.
To map the linear regime, we therefore compared the amplitudes perturbations. In practice, it is a priori unclear which perturbation amplitudes can be considered sufficiently small. In the framework used here, we can quantify this by analyzing the frequency contents of the response r(t). In the linear regime, a sinusoidal stimulus of Nordlie et al. LIF: strong synapses indeed bounded from above by the synaptic cutoff frequency  f s . For supercritical weights w r > 1, f c depends on  f s in a nonlinear (but still monotonous) fashion. For small τ s , the cutoff frequency of the rate transfer can considerably exceed the synaptic cutoff frequency. In other words, the time constant τ = 1/(2πf c ) can be much smaller than the synaptic time constant τ s .
As shown in Figure 7B, the low-frequency gain γ is, as a good approximation, given by the slope g′(a 0 ) of the activation function, independently of the choice of the time constants τ s and τ m , and the synaptic weight w r .

test of the lInear-nonlInear Model
By systematically varying the parameters f stim , a 0 and a 1 of the sinusoidal input rate, we obtained an activation function g(·) (Figure 3) and a linear state-independent transfer kernel h(t) (Figures 4, 5, and 7) for each parameter set {w,τ s ,τ m }. With these ingredients we set up firing-rate models of the form (9). The model extraction is based on the period-averaged responses of a single neuron to sinusoidal inputs, a procedure which is often applied in experimental neuroscience. Provided the spike responses in consecutive stimulus periods are independent the resulting firing-rate model should also yield accurate predictions for the instantaneous population-averaged response of an ensemble of independent (unconnected) neurons. Figure 8 illustrates that this is in fact the case for subcritical synaptic weights w r < 1 (Figures 8A-I): here, the population-averaged step response of an ensemble of 50000 independent LIF neurons (black curves) is well described by the linear-nonlinear model (gray curves). Both the time constants τ = 1/(2πf c ) of the rate dynamics (9) and the stationary response rates are correctly predicted. Similarly good agreement between the rate-model predictions and the results of direct simulation was obtained for a variety of other test stimuli.
In a few cases, we observe some subtle differences between model predictions and simulation results: For small τ s and subcritical weights, the time constants of the model are slightly underestimated (see e.g., of the principal and second harmonics r 1 and r 2 , respectively, for all neuron and stimulus parameters. Figure 6 shows the results for a mean input rate a 0 = 40 s −1 and a modulation frequency f stim = 10 Hz (similar results were obtained for other stimulus frequencies and input rates). In Figure 6A, the amplitudes of the first two harmonics are shown for three different synaptic weights. For sufficiently small input amplitudes a 1 , the amplitude r 2 of the second harmonics does not exceed the noise level. Note, that r 2 is never exactly zero because of a finite observation time T (see Section 2.2). The ratio r 1 /(r 1 + r 2 ) shown in Figure 6B is thus never exactly 1. For weights w r > 0.6 and intermediate input amplitudes a 1 , the rate transfer is approximately linear. For synaptic weights w r < 0.6, we observe an increase in r 2 for larger modulation depths ( Figure 6B) which is due to a rectification of the response modulation (r 1 ≥ r 0 ; see dotted curves in Figure 6B).

dependence on synaptIc and MeMbrane tIMe constants
At first sight, one might expect that the tracking speed of a population of neurons is limited by the low-pass filter properties of the neuronal membrane, i.e., the membrane time constant τ m . Previous studies (e.g., Gerstner, 2000) have however shown, that the cutoff frequency of the rate transfer is rather determined by the synaptic time constant τ s . Here, we test this for a number of time constants τ s ∈ {0,5,1,2,5} and τ m ∈ {10,20} ms. In Figure 7A, the cutoff frequency f c of the rate transfer normalized by the cutoff frequency of the synaptic kernel  f s is plotted as a function of the relative synaptic weight w r for different synaptic and membrane time constants τ s and τ m , respectively. For the whole range of tested synaptic weights, the cutoff frequencies do not significantly differ for the two considered membrane time constants τ m = 10 ms (solid curves) and 20 ms (dashed curves). In the subcritical regime w r < 1, the firing-rate cutoff frequency f c is close to or slightly smaller than the synaptic cutoff   Figure 8A). For supercritical weights w r > 1, the population-averaged step response sometimes exhibits an overshoot immediately after the jump in the input rate which is not captured by the model (Figures  8J,M). Further, the model overestimates the rate transmission delay d in some cases (Figures 8K,L,N). It is likely that these deviations for small synaptic time constants or supercritical weights are a result of the dependence of the cutoff frequency on the choice of the working point a 0 (see Figure 5C). The averaging strategy leading to the state-independent kernel h(t) ignores this dependence.

dIscussIon
In the present study, we investigated how LIF neurons transmit presynaptic time-varying firing rates to the output rate in a period-or population-averaged sense. In contrast to earlier studies, we focused on a regime where the synapses have moderate or large synaptic weights and where the average input rates are relatively low. The response properties were characterized by measuring the activation function g(a 0 ), i.e., the response r 0 to constant input rates a 0 , and the transfer function H 0 ( f ) describing the transmission of time-varying rate perturbations in the vicinity of a working point a 0 . Both concepts, the activation function and the transfer function, were used to set up a firing-rate model which was evaluated by comparing its predictions with simulation results for non-sinusoidal test stimuli.
In agreement with previous work (Sirovich, 2008), we found that even for strong synapses only a fraction of the input spikes is transmitted. Only for (presumably) unrealistically large synaptic weights above a critical value w crit , the transfer ratio r 0 /a 0 approaches or even exceeds one. We have observed this behavior for both stationary and nonstationary inputs.
In the diffusive regime, i.e., for weak synapses, different transfer scenarios have been described in the literature for different neuron and stimulus models. All previous studies predict some sort of low-pass behavior (at least if one takes into account that incoming spikes are low-pass filtered by the synapses). We observe similar behavior for large synaptic weights. If the synapses are so strong that a single isolated input spike can reliably cause a response spike in the postsynaptic neuron, i.e., for w > w crit , one would expect that any signal carried by the presynaptic spikes can be transmittedno matter what its frequency contents is and independently of the synaptic filter properties. Both weights larger than w crit and infinitesimally small weights are unrealistic for most biological neural systems. It is tempting to interpolate between these two extreme cases and assume that the ability of a neuron population to track fast stimuli gradually increases with the synaptic strength w. In fact, Helias et al. (2010) have recently shown that for finite synaptic weights the population response of an ensemble of integrate-andfire neurons exhibits an instantaneous component which gradually increases with the synaptic weight. The analysis presented in (Helias et al., 2010) is however still restricted to small synaptic weights. For strong but subcritical weights w < w crit , we found that the cutoff frequency f c of the rate transfer is bounded from above by the synaptic cutoff frequency  f s . Only for supercritical weights w > w crit , the tracking speed of the population can be significantly increased. In agreement with earlier work (e.g., Gerstner, 2000), our data suggest that the cutoff frequency f c of the rate transfer does not depend on the choice of the membrane time constant.
In most previous studies, the stimulus is treated as a modulation of the distribution of postsynaptic somatic currents. In our study, in contrast, the stimulus is represented by a modulation in the presynaptic firing rate. Here, the effective transfer function is a combination of the neuronal (the current-to-rate transfer) and the synaptic filtering (rate-to-current transfer). For correlated noise, it has been shown that ensembles of LIF neurons can faithfully transmit fast modulations in the mean or the variance of the input current (see e.g., Brunel et al., 2001;Lindner and Schimansky-Geier, 2001): the response amplitudes remain finite even at high frequencies. In this case, the net filtering at high frequencies is governed by the synaptic kernel h s (t) alone. For alpha-function shaped synaptic currents (as used in this study), one would therefore expect second-order low-pass characteristics with response amplitudes decaying as 1/f 2 at high frequencies (see (3)). Note, however, that this result refers to the diffusion limit, i.e., infinitesimal synaptic weights. For non-diffusive regimes, the shape of the neuronal transfer function is unknown. Here, we for simplicity describe the data by a first-order low-pass filter. For most parameters, this results in reasonable fits (see, for example, Figure 4). In practice, however, the high-frequency response can easily be affected by noise: due to the stochastic nature of the input spike trains and finite kernel explicitly into account. As shown in Figure 8, the predictive power of the reduced model with a state-independent kernel is however acceptable even for fast synapses and low firing rates. For small synaptic weights, previous studies have shown that the characteristics of the transmission of time-varying input currents to output rates crucially depends on the properties of the neuron model (Knight, 1972a;Fourcaud-Trocmé et al., 2003;Naundorf et al., 2005). Our study is restricted to the LIF model. Future work needs to show how more (or less) biologically realistic neuron models alter our findings on the transmission of input firing rates (rather than input currents) with finite synaptic weights. In our and in most previous studies on the transfer characteristics of neuron ensembles, synapses are modeled as static kernels describing the PSC. The effect of conductance-based (see e.g., Kuhn et al., 2004, and references therein) or dynamical synapses (Tsodyks et al., 1998), for example, on the firing-rate filter properties of a neuron ensemble is so far unclear. Another restriction in our work is the neglect of inhibitory inputs. Pressley and Troyer (2009) recently pointed out that the rate transfer of a population of perfect integrate-and-fire neurons depends on the interplay between presynaptic excitation and inhibition.
A number of experimental studies has systematically characterized the stationary (e.g., Arsiero et al., 2007) and nonstationary response properties (Köndgen et al., 2008;Boucsein et al., 2009) of populations of neurons in vitro. As we have suggested here, the combination of both transfer characteristics can -in principlelead to powerful firing-rate models which allow the description of neuronal activity on a macroscopic scale. We (and others) pointed out that the transfer characteristics are to a large extent determined by the properties of the synapses. The above-mentioned experimental studies are based on direct somatic current injections. The role of the synapses has thus been ignored so far. The controlled stimulation of a large number of individual synapses in in vitro (and in vivo) experiments is difficult. Recent developments in experimental techniques (e.g., Boucsein et al., 2005), however, give hope that this problem may be overcome in the near future.

acknowledgMents
We thank Hans Ekkehard Plesser and the two reviewers for constructive discussions and valuable comments on the manuscript. We acknowledge partial support by the Research Council of Norway (eNEURO, NOTUR) and the Honda Research Institute Europe. Simulations and data analysis were carried out using the NEST simulation tool 9 and Python 10 .
observation times, small amplitudes are in general overestimated. This may lead to a flattening of the transfer function. Here, we do not address the exact shape of the transfer function at high frequencies. Instead, we restrict this study to a description of the cutoff frequency and the low-frequency gain in order to design a simple firing-rate model. We tested that these characteristics are rather robust against the choice of the transfer function 7 . Since the work of Wilson and Cowan (1972), firing-rate models with exponential kernels (i.e., first-order low-pass filters) have become a standard tool in neuroscience 8 . In the present study, we show that the presence of strong synapses does not necessarily disqualify such models. For an ensemble of unconnected LIF neurons, we demonstrated that a simple 1st-order model yields accurate predictions of the population-averaged response for a wide range of stimulus, neuron, and synapse parameters. For recurrent systems, however, the exact shape of the rate transfer function at high frequencies can be crucial: Nordbø et al. (2007), for example, have shown that the shape of the rate-transfer kernel determines the stability of persistent states in neural-field models. Solutions which are stable for exponential kernels can become unstable for alpha-function kernels. In future studies, the performance of the rate model extracted here therefore needs to be tested under recurrent conditions.
The approach presented in this study can in principle be applied to any nonlinear system. In general, this leads to transfer functions which depend on the choice of the working point (a 0 , r 0 ) in a non-trivial way. For large but subcritical synaptic weights and not too small synaptic time constants, we observed that the normalized rate-to-rate transfer function H(f ) = H 0 (f )/g′(a 0 ) is approximately independent of the working point and mainly determined by the synaptic filter. Here, a simple linear-nonlinear model r(t) = g((a *h (t)) with a state-independent kernel h(t) yields accurate predictions of the population-averaged firing rate. In accordance with previous studies (e.g., Knight, 1972a;Brunel et al., 2001;Lindner and Schimansky-Geier, 2001;Mattia and Del Giudice, 2002), we found that for small synaptic time constants τ s the cutoff frequencies f c increase with the firing rate. If the synaptic cutoff frequency  f s is small compared to the firing rate r 0 (i.e., for larger synaptic time constants), the synaptic filter seems to dominate. If  f s is large, a precise prediction of the response rate would require a model which takes the state-dependence of the transfer 7 Fitting of a second-order filter to the data led to similar result as the first-order model (data not shown). 8 For exponential kernels h(t), the linear-nonlinear model (9) can be transformed into a one-dimensional differential Eq. (10). The mathematical and numerical analysis of the model is thereby considerably simplified. 9 www.nest-initiative.org 10 www.python.org