Abstract
This study examines the ability of neurons to track temporally varying inputs, namely by investigating how the instantaneous firing rate of a neuron is modulated by a noisy input with a small sinusoidal component with frequency (f). Using numerical simulations of conductance-based neurons and analytical calculations of one-variable nonlinear integrate-and-fire neurons, we characterized the dependence of this modulation on f. For sufficiently high noise, the neuron acts as a low-pass filter. The modulation amplitude is approximately constant for frequencies up to a cutoff frequency, f_{c}, after which it decays. The cutoff frequency increases almost linearly with the firing rate. For higher frequencies, the modulation amplitude decays as C/f^{α}, where the power α depends on the spike initiation mechanism. For conductance-based models, α = 1, and the prefactor C depends solely on the average firing rate and a spike “slope factor,” which determines the sharpness of the spike initiation. These results are attributable to the fact that near threshold, the sodium activation variable can be approximated by an exponential function. Using this feature, we propose a simplified one-variable model, the “exponential integrate-and-fire neuron,” as an approximation of a conductance-based model. We show that this model reproduces the dynamics of a simple conductance-based model extremely well. Our study shows how an intrinsic neuronal property (the characteristics of fast sodium channels) determines the speed with which neurons can track changes in input.
- populations of spiking neurons
- noise
- dynamics
- sodium channel
- integrate-and-fire model
- conductance-based model
Introduction
The input-output transformation performed by a neuron is classically characterized by its frequency-current (f-I) relationship (McCormick et al., 1985; Powers and Binder, 2001). However, knowing of this relationship is not enough to predict neuronal responses to transients or time-dependent inputs, a problem that is, in general, particularly tricky because of the highly nonlinear nature of action potential generation. Nevertheless, if the temporal variations of the input are small enough, the dynamics are dominated by the linear response properties of the neuron. Therefore, the response of the neuron to an arbitrary input with weak temporal variations can be predicted if we can determine how its instantaneous firing rate is modulated by sinusoidal inputs.
The response of neurons to sinusoidal inputs has been investigated in various in vitro preparations, including the horseshoe crab Limulus polyphemus (Knight, 1972b), visual cortex (Carandini et al., 1996; Nowak et al., 1997), and vestibular nuclei (Ris et al., 2001; Sekirnjak and du Lac, 2002). Most of these studies restricted their analysis to low-frequency inputs (<20 Hz). Carandini et al. (1996) found that the response of regular-spiking cortical cells to injection of a broadband noisy input is linear and can be flat up to 100 Hz. Bair and Koch (1996), using in vivo recordings in middle temporal (MT) cortex of anesthetized monkeys, showed that the power spectra of the responses of neurons to randomly moving dots are low-pass with a broad range of cutoffs up to 150 Hz.
The linear response of neurons to noisy fluctuating inputs has also been investigated in theoretical studies relying on the leaky integrate-and-fire (LIF) model (Knight, 1972a; Gerstner, 2000; Brunel et al., 2001; Fourcaud and Brunel, 2002; Mazurek and Shadlen, 2002). It was found that in the presence of white noise, LIF neurons behave like low-pass filters, with a cutoff depending on the passive membrane time constant and the average firing rate of the neuron. The gain of the filter decays as 1/√f, where f is the input frequency, and its phase shift reaches the value of 45°, at sufficiently large f. It was also found that temporal correlations in the input noise improve the “accuracy” of the LIF response, because they suppress the decay at high frequencies and reduce the phase shift (Gerstner, 2000; Brunel et al., 2001; Fourcaud and Brunel, 2002). In the present study, we show that simple conductance-based (CB) models behave differently. Using a combined analytical and numerical approach, we investigated theoretically how neuronal properties affect the response of neurons to fluctuating inputs. In particular, we demonstrate that the suppression of the AC part of the spike response to fast varying inputs depends on the DC response of the neuron and on the characteristics of the fast sodium currents responsible for spike initiation.
Materials and Methods
Linear response of the instantaneous firing rate. We investigated the response of a neuron to a time-varying input. Our goal is to determine the instantaneous firing rate ν(t) (i.e., the probability that the neuron will fire between time t and time t ± dt, divided by dt, in the limit dt → 0). We assume that the neuron receives an input 1 where I_{0} is the DC component of the input, δI(t) is its deterministic temporal variation around I_{0}, and I_{noise} is a noisy component that is described in more detail below. The firing rate ν(t) is determined as a function of I_{0} ± δI(t) by averaging the response of the neuron over different realizations of I_{noise}. Equivalently, ν(t) represents the instantaneous firing rate of a population of noninteracting neurons that receive the same input I_{0} + δI(t), but in which each neuron receives an independent realization of the noise.
In the limit of small δI(t), a linear relationship can be written between ν(t) and δI(t): 2 where Φ(I_{0}) is the average firing rate of the neuron in response to a current I_{0} in the presence of noise I_{noise} (the f-I curve of the neuron), and the kernel, κ, depends on I_{0} and the properties of the noise. In particular, the response to a sinusoidally and weakly modulated I_{syn}(t): 3 where I_{1} << I_{0} is: 4 where ν_{0} = Φ(I_{0}), and ν_{1}(f) and ϕ(f) are related to the modulus and the argument of the Fourier transform of κ, , by and .
The gain of the response at frequency f, ν_{1}(f)/I_{1}, and the phase response, ϕ(f), completely characterize the linear filtering properties of the neuron. We compute these quantities as a function of the frequency f. Note that at sufficiently low frequencies, it is straightforward to relate the gain ν_{LF}/I_{1} to the f-I curve of the neuron. Indeed, for small f, one can write ν(t) = Φ[I_{0} + I_{1} cos(2πft)], so that because I_{1} << I_{0}: 5 whereas the phase shift at small f is ϕ_{LF} = 0.
Noise. When a neuron receives a large number of synaptic inputs per membrane time constant through synapses with small amplitude compared with threshold, the resulting currents can be well approximated by a mean input current with random Gaussian variations around the mean (diffusion approximation) (Tuckwell, 1988); thus, we consider Gaussian noise. If one assumes that postsynaptic currents are instantaneous, the Gaussian noise is white. For exponentially decaying synaptic currents, the Gaussian noise is colored with a correlation time equal to the synaptic decay time constant (see Appendix B2).
Conductance-based models
Two single-compartment neuronal models are used in this study. The first was proposed by Wang and Buzsáki (1996), which is a modified version of the original Hodgkin-Huxley model (Hodgkin and Huxley, 1952). The Wang-Buzsáki (WB) model has a leak current, a fast sodium current with instantaneous activation dynamics, m(V), and a delayed-rectifier potassium current. The firing current threshold of the neuron is I_{T} = 0.16 μA/cm^{2}. This corresponds to a voltage threshold V_{T} = -59.9 mV. Note that the WB model is type I (zero firing frequency at threshold) and not type II (nonzero firing frequency at threshold), like the Hodgkin-Huxley model with standard parameters (Hodgkin and Huxley, 1952) because of the different choice of parameters, in particular for the sodium and potassium conductance (Wang and Buzsáki, 1996).
The second model was proposed by Hansel and van Vreeswijk (2002), which has the fast sodium current and the delayed rectifier current of the WB model but also additional intrinsic currents: A-type current, slow potassium current, and persistent sodium current. Details of the models are given in Appendix A.
Numerical integration of CB models. Model equations in the presence of noise have been integrated using a stochastic second-order Runge Kutta (RK2) method (Honeycutt, 1992), with a time step <0.02 msec in all simulations. The instantaneous firing rate was computed by taking the average spike count in bins of one-thirtieth of the input period over a minimum of 2000 periods. The obtained instantaneous firing rate was fitted by a sinusoid (Eq. 4) using a least-square method or by computing directly the first Fourier component of the output spike train. Both methods yielded the same results. The input modulation amplitude I_{1} was adjusted depending on the input frequency to get an instantaneous firing rate modulation equal to 25% of the baseline. We checked that, to this level of input modulation, the amplitude of the firing rate modulation ν_{1} varied linearly with I_{1}.
Nonlinear integrate-and-fire neuron models
To understand the factors that determine the response to time-varying inputs, we investigated simpler models that are to some extent analytically tractable. We introduced a family of one-variable neuronal models whose membrane potential V dynamics are given by: 6 where C is the membrane capacitance, g_{L} is the leak conductance, V_{L} is the leak potential, I_{syn} is the external synaptic current (Eq. 3), and ψ(V) is a function of voltage that describes the spike-generating currents. The passive membrane time constant is τ_{m} = C/g_{L}.
The LIF (Knight, 1972a; Tuckwell, 1988) is the special case of this family of models for which ψ(V) = 0. In this model, a spike threshold V_{th} has to be imposed to obtain spike generation. The resulting spike is instantaneous, and the neuron is reset to a voltage V_{r} after a spike.
When the function ψ(V) is supralinear, the membrane potential diverges to infinity in finite time if the input current exceeds some threshold. This divergence can be identified with the firing of a spike, provided that one supplements Equation 6 with a reset condition for the membrane potential to a value V_{r}, as in the LIF model. A threshold voltage V_{T} can be defined as the voltage at which the slope of the I-V curve vanishes, which is given by: 7 This is also the largest steady voltage at which the neuron can be maintained by a constant input current. The corresponding current I_{T} = g_{L}(V_{T} - V_{L}) - ψ(V_{T}) is the threshold current above which tonic firing occurs. We also defined the spike slope factor Δ_{T} (in mV) as: 8 This parameter, which is inversely proportional to the curvature of the I-V curve at the threshold V_{T}, measures the sharpness of spike initiation. The quadratic integrate-and-fire (QIF) neuron corresponds to: 9 where I_{T} is the current threshold, V_{T} and Δ_{T} are the spike threshold and slope factor, respectively, as defined above. If one chooses V_{r} = -∞, this model represents the normal form of type I neurons (Ermentrout and Kopell, 1986; Ermentrout, 1996).
Other spike-generating currents, ψ(V), were considered in this study. Of particular interest is the exponential integrate-and-fire (EIF) neuron for which the spike-generating current is exponential: 10 It is easy to check that V_{T} and Δ_{T} satisfy Equations 7 and 8, respectively. Note that in the limit Δ_{T} → 0 (spike with very sharp initiation), the EIF neuron becomes equivalent to the LIF model with V_{th} = V_{T}.
Fitting IF to the Wang-Buszáki model. In all of the models considered in this study, we take C = 1 μF/cm^{2}, g_{L} = 0.1 mS/cm^{2}, and V_{L} = -65 mV. The spike threshold V_{th} of the LIF model was chosen such that the threshold current is equal to the threshold current I_{T} of the WB model. In the WB model, the f-I curve behaves near the onset of firing (I ≳ I_{T}) as β√I-I_{T}. This is also the case for the QIF and EIF models. Therefore, one can determine the parameters β and I_{T} to match the behavior of their f-I curves with one of the WB models. This gives β = 0.038 msec^{-1} μA^{-1/2} cm, I_{T} = 0.16 μA/cm^{2}. This also determines Δ_{T} = 3.48 mV and V_{T} = -59.9 mV. For larger I, the f--I curves are no longer well described by the square root behavior. The reset voltage, V_{r}, and the refractory period, τ_{ref}, were chosen to minimize the difference of its f-I curve with one of the WB models for input currents well above I_{T} (Fig. 1). Matching the f-I curves in this range of currents determines V_{r} = -63.8 mV and τ_{ref} = 0 (QIF model), V_{r} = -68 mV and τ_{ref} = 1.7 msec (EIF model).
Fokker-Planck equation. To analyze the response of an IF neuron to a time-dependent input, we studied how the distribution of its voltage, P(V, t), evolves over time. We assumed Gaussian white noise, I_{noise}(t) = σ √C_{gL} η(t), where η(t) has unit variance and σ has a dimension of a voltage. The distribution, P(V, t), at time t obeys the Fokker-Planck equation (Risken, 1984; Abbott and van Vreeswijk, 1993; Brunel and Hakim 1999; Brunel, 2000; Knight et al., 2000; Nykamp and Tranchina, 2000; Brunel et al., 2001; Fourcaud and Brunel, 2002): 11 where is the probability flux through the potential V (i.e., the probability per unit of time that the membrane potential of the neuron crosses the value V from below at time t). The firing rate, ν, of the neuron is the flux through V = ∞: 12 For a complete analytical description, one also needs the boundary conditions, which take into account that, after the spike, the neuron is reset to V_{r}.
All of our calculations are based on the above equations, which are valid for white input noise, or on their generalizations for more realistic noise models. The details of the calculations can be found in Appendix B.
Simulations of the EIF model. In the EIF model, an action potential is defined as a divergence of the voltage. In simulations, one has to introduce a cutoff at a finite voltage V_{th}. However, if one assumes that the neuron immediately spikes after reaching V_{th}, one ends up with a neuron that is effectively an LIF for some inputs. This can be avoided by treating the spike more carefully. In our simulations, we chose a threshold V_{th} large enough so that: 13 by a factor ≥100. For V < V_{th}, we integrated the dynamics using a stochastic RK2 algorithm (Honeycutt, 1992). For V > V_{th}, we neglected all of the currents but ψ(V) and analytically integrated Equation 6. In particular, the time it takes for the voltage to increase from V = V_{th} to V = +∞ is t_{sp} = τ_{m} exp[(V_{T} - V_{th})/Δ_{T}]. If V_{th} is too large, the RK2 procedure at large V underestimates the increase in the membrane potential at each time step. This leads to a systematic overestimate in the spike time and causes an additional phase shift in the neuronal response. In practice, we must make sure the time t_{sp} is large compared with the integration time step dt. Thus, one needs to check that the two conditions, (Eq. 13) and t_{sp} >> dt, are both satisfied. In the simulations with Δ_{T} = 3.48, a value for which both conditions are satisfied with V_{th} = -30 mV for the range of inputs that we investigated.
Results
Filter of a simple conductance-based model
We investigated how the instantaneous firing rate of the WB model is modulated by a sinusoidal input at frequency f in the presence of a noisy background input using numerical simulations. Figure 2A shows a particular realization of such an input current. The instantaneous firing rate is computed by averaging the response over many realizations of the noise, as shown in Figure 2, B and C. Figure 2, D and E, shows how the amplitude and phase shift of the firing rate modulation of the neuron depends on the frequency of the sinusoidal input when the noisy component is a Gaussian white noise with an SD of σ = 6.3 mV and its DC part is such that the average firing rate is ν_{0} = 40 Hz.
The amplitude versus frequency curve (Fig. 2D) shows a marked attenuation of the response above a frequency of ∼50 Hz. A fit of our numerical simulation data in the range 50 Hz < f < 500 Hz, assuming a power law decay, reveals that the firing rate modulation at large f decays as 1/f (Fig. 2D). Moreover, we find that in this frequency range, the phase shift does not vary much and is slightly larger than 90°, as shown in Figure 2E. We performed extensive numerical simulations varying the level of noise and the DC external input as well as the temporally correlated noise. In all of these simulations, the same behavior of the firing rate modulation was found at sufficiently large f, although the frequency above which it becomes apparent depends on the values of the varied parameters (Fig. 9).
In all of these simulations, the instantaneous spike rate was computed using a definition of spike time as the time the membrane potential crosses a voltage, V_{s} = 20 mV, from below. The response modulation amplitude does not depend on the precise value of V_{s}, provided it is sufficiently high (more than -30 mV). In contrast, the phase shift is more sensitive to V_{s}, especially at very high input frequencies where the increase of the phase shift is linear with a slope that depends on V_{s} (Fig. 2E).
Integrate-and-fire models with intrinsic spike generation mechanisms
Recent studies have shown that the high-frequency response of the LIF model depends on the temporal correlations of the noise. For white noise, the response decreases as 1/√f with a phase lag of 45°; for colored noise, the response stays finite with no phase lag (Brunel et al., 2001; Fourcaud and Brunel, 2002). These results are clearly at odds with what we found in the simulations of the WB model.
To understand the origin of the behavior of the WB model at high frequencies and why it is different from the LIF, we studied nonlinear integrate-and-fire neurons, in which the membrane potential obeys Equation 6. The advantage of using this simplified model is that it allows us to compute analytically both the f-I curve and the firing rate modulation at high frequencies, and hence to understand in detail the factors that affect it.
With a supralinear spike-generating current ψ(V), the neuron is endowed with an intrinsic mechanism for spike initiation. We considered two cases in particular, the QIF neuron defined by Equation 9 and the EIF neuron defined by Equation 10. These two neuronal models are type I, like the WB model, and their f-I curve increases proportionally to √^{I-I}_{T}, near the firing onset at current threshold, I_{T}. Hence, for both the QIF and EIF, the parameters can be chosen so that the onset of firing in response to a constant input reproduces the parameters of the WB model sufficiently close to the current threshold (Fig. 1) (see Materials and Methods for details). For the EIF model, this constraint determines the value of the “spike slope factor,” Δ_{T}. To determine the remaining free parameters, V_{r} and τ_{ref}, we also required a good fitting of the WB model f-I curve outside the region of bifurcation.
The f-I curves of the QIF and EIF, with parameters so defined, and the WB model are compared in Figure 1A. It shows that the firing rate response to constant current of the EIF and WB models are very close in a broad range of firing rates (up to 200 Hz). The QIF neuron reproduces the f-I curve of the WB neuron well at low firing rates but not at high rates. Finally, the LIF neuron cannot reproduce the f-I curve of type I neurons well at low rates, because the firing rate has a logarithmic, not square root, dependence on I - I_{T}. However, the LIF neuron can be made to reproduce the f-I curve reasonably well at high firing rates by a suitable adjustment of parameters.
We then compared the response of the QIF, EIF, LIF, and WB models to noisy input currents (Eq. 3, with I_{1} = 0). The corresponding f-I curves are plotted in Figure 1B for the four models with the same parameters as in Figure 1A. As in the noiseless case, the f-I curves of the QIF and EIF models match the WB model f-I curve very well in a broad range of input currents, and this range is larger for the EIF model than for the QIF model.
For a more in-depth look at the spiking dynamics of these various models, Figure 3 plots the voltage traces of the four models for the same realization of the input current. As shown in Figure 3A, the QIF, EIF, and LIF models behave on large time scales in a very similar way as the WB model. However, differences between these models are found on shorter time scales, as shown in Figure 3B. The best match with the WB model is obtained with the EIF model. The subthreshold potential traces are essentially indistinguishable. This is because the same leak current is present in both models and because the I-V curves of the two models are very similar near threshold voltage, V_{T}. The latter property stems from the fact that the activation curve of the fast sodium current in the WB model can be well approximated near threshold by an exponential function, as shown in Figure 3C. Thus, the membrane potential behavior at the spike onset is very similar for the WB and EIF models. After initiation, the dynamics of the spike in the WB model is controlled by the interplay between Na^{+} and K^{+} channel dynamics and is essentially independent of synaptic inputs. Thus, its spike shape is approximately invariant. Consequently, the only significant difference between the EIF and WB models, the precise shape of the spike after the onset, leads to an input-independent fixed short delay (of order 0.1-0.2 msec) between the spikes of the two models.
The other two models failed to reproduce the dynamics of spike initiation of the WB model. When a fluctuation in the input current induces a spike in all four models, the LIF neuron, which has an instantaneous spike, usually spikes first. The spike in the QIF neuron usually occurs later because of the long time spent by the neuron in the voltage threshold.
The high-frequency response of nonlinear integrate-and-fire neurons
We analytically calculated the firing rate response of a general IF model to a sinusoidally modulated input in the presence of Gaussian white noise (Eq. 3). The details of the calculations, which rely on the Fokker-Planck equation and an expansion of both instantaneous firing rate and membrane potential distribution in powers of 1/f, are described in Appendix B. There, we prove that, depending on the function ψ(V), there are three classes of possible behaviors for the firing rate modulation at high frequencies:
If ψ′(V)/ψ(V) does not vanish and is finite at large V, the phase lag ϕ(f) is 90°, whereas the response modulation is given by: 14 where ψ′(V) is the derivative of ψ. This happens if the spike-generating current has an exponential dependence on the voltage, as in the EIF model. For this model, one finds that at high frequencies: 15 and the phase lag is ϕ = 90°.
If ψ′(V)/ψ(V) goes to zero at large V, the rate modulation decays faster than 1/f at high frequencies. This is, for instance, the case of the QIF model for which the calculation of the order 1/f^{2} in the expansion of ν_{1} shows that at high frequencies, ν_{1}(f) ≈ (ν_{0}I_{1})/[g_{L}ΔT(2πfτ_{m})^{2}], and the phase lag is ϕ = 180° (see Appendix B).
If ψ′(V)/ψ(V) diverges when V goes to infinity, the response is less attenuated than 1/f. This happens in particular for the LIF neuron.
We also show in Appendix B that in cases 1 and 2, the high-input frequency behavior is independent of the temporal correlations of the noise. In contrast, in case 3, the properties of the noise matter, as found in previous work in which the LIF model was studied (Brunel et al., 2001; Fourcaud and Brunel 2002). Furthermore, it can be shown that the results hold for current-based as well as conductance-based noise, because in the case of suprathreshold membrane potential, the spike-generating current ψ(V) dominates the neuronal dynamics, and the synaptic input has little effect on what happens after spike initiation.
In conclusion, the amplitude of the firing rate modulation behaves at high frequencies like a power law ν_{1}(f) ≈ 1/f ^{α}, where the exponent α depends on the nonlinearity of the spike-generating current ψ(V). The phase lag at high-input frequency is απ/2. The exponent α and asymptotic phase shift of the various models are summarized in Table 1. The high frequencies behaviors of the rate modulation in WB and EIF models are extremely consistent. Other IF models failed to reproduce the behavior of the WB model.
Properties of the EIF neuron at intermediate frequencies
The simulation and analytical results described so far show that in the limit of low- and high-input oscillation frequencies, the behaviors of the EIF and WB models are very similar. We consider now the behavior of the EIF model at intermediate frequencies using numerical simulations.
When the noise level is weak, resonances appear at frequencies that are integer multiples of the average firing rate ν_{0}, f = kν_{0} (k = 1, 2,...) (Fig. 4, insets). This is because at low noise levels, the neuron behaves like an oscillator. When the noise increases, the neuronal firing process becomes less regular. The resonant peaks become less pronounced and they disappear, one after the other. For a sufficiently high noise level, the last resonance peak (k = 1) disappears, and the firing rate modulation monotonically decreases with the input frequency, f. The dashed line in Figure 4, computed by numerical simulations, represents the boundary in the plane (ν_{0}, σ) between the weak noise regime (below the line) where the firing rate modulation displays resonances and the strong noise regime (above the line) where the resonances are completely suppressed. In the latter region, the spike response modulation, ν_{1}(f), is approximately constant in some range of input frequencies, beyond which it decreases with f with an asymptotic decay in 1/f. Therefore, when the noise is sufficiently large, the EIF neuron behaves like a low-pass filter.
The gain and phase shifts of the filter in the presence of white noise are plotted in Figure 5, A and B, respectively, for fixed ν_{0} and three values of the spike slope factor, Δ_{T}. In all cases, the gain of the filter decreases as 1/f at large enough frequencies. This asymptotic regime is reached earlier for larger Δ_{T}. There is an overall reduction of the gain as Δ_{T} increases. This effect is in agreement with our analytical calculation (Eq. 15), which predicts that in the limit f → ∞, the firing rate modulation, ν_{1}, decreases with Δ_{T}. This is also consistent with Equation 5, because one can show (using Eq. 29) that the derivative of the f-I curve at ν_{0} = 20 Hz is a decreasing function of Δ_{T}. The phase shift is an increasing function of Δ_{T}, as shown in Figure 5B. This corresponds to what is expected intuitively, because a decrease of Δ_{T} increases the sharpness of spike initiation. At large frequencies, the phase shift of the filter goes to -π/2 for any Δ_{T}, but the convergence is slower when Δ_{T} is small.
Figure 6 shows the effect of temporal correlations in the input noise on the shape of the filter. The response is plotted for three values of the temporal correlation of the noise, τ_{s}, over a broad range of input frequencies, f. In these simulations, the noise level was adjusted so that the low-frequency response stays constant as τ_{s} was varied. All three curves overlap for sufficiently large frequencies, f > 100 Hz, in agreement with our analytical results. In contrast, for intermediate frequencies, 30 < f < 100 Hz, the response modulation increases with τ_{s}. This effect is reminiscent of the behavior of the LIF neuron (Brunel et al., 2001) for which correlations in the input noise suppresses the attenuation of the response modulation (Brunel et al., 2001; Fourcaud and Brunel, 2002).
To further characterize the filtering properties of the neuron, we defined the filter cutoff as the frequency, f_{c}, for which the filter gain decreases by a factor 1/√2 compared with its low-frequency limit, ν_{1}(f_{c}) = ν_{LF}/√2. We used numerical simulations to study how f_{c} depends on the average firing rate, ν_{0}, the spike slope factor, Δ_{T}, and the temporal correlations in the noise. Figure 7 shows the results for a noise level σ = 8 mV. For this value of σ, the EIF neuron is in the strong noise regime in the entire range of average firing rates and the spike slope factor explored.
Figure 7A shows that for a fixed value of the slope factor, f_{c} is finite in the limit of low average firing rate, ν_{0}, and that it is an increasing function of ν_{0}. For ν_{0} > 10 Hz, f_{c} varies linearly with ν_{0}. In Figure 7B, the cutoff frequency is displayed as a function of Δ_{T} for fixed ν_{0} and different values of τ_{s}, the noise correlation time. In general, f_{c} is a decreasing function of Δ_{T} but, at large enough Δ_{T}, it saturates toward a value that depends only weakly on τ_{s}. In contrast, the smaller the τ_{s}, the weaker the variations of f_{c} at small Δ_{T}. In particular, for white input noise, f_{c} is almost constant in all of the explored range of Δ_{T}. Temporal correlations in the noise increase the cutoff frequency, f_{c}, and this effect is more pronounced at smaller values of Δ_{T} (i.e., when spikes have sharper initiation).
Domain of validity of the linear response
In the previous sections, we analytically computed the linear response of the EIF neuron to modulated input in the limit of small or large modulation frequency, f. These results should provide a correct quantitative description of the neuron dynamics in these limits for sufficiently weak amplitude modulation, I_{1}. The numerical simulations presented above, performed for a finite but small value of I_{1}, confirm this expectation. They also extend our study by investigating the neuronal linear filter at intermediate values of the input frequency.
How large can the input modulation be for the response of the neuron to still be linear? To address this issue, we performed simulations at various levels of input modulation and measured the amplitude of the Fourier components of the response up to order 3. Results of these simulations are displayed in Figure 8 for σ = 6.3 mV, ν_{0} = 20 Hz, and an input modulation frequency of f = 20 Hz. For these parameters, ν_{1} varies linearly with I_{1}, and higher order Fourier components are negligible up to input modulations that induce temporal variations of the firing rate as large as 75% of the average firing rate. This shows that the neuron responds in a linear manner in this range of firing rate variations. A similar range of validity was obtained for all values of f tested (0-100 Hz) in the strong noise regime (the boundaries of strong and low-noise regimes are indicated in Fig. 4). In the low-noise regime, the domain of validity is on the same order of magnitude, except for frequencies close to resonances.
Comparison of the linear responses of EIF and conductance-based neurons
In Figure 9, we compare the gain and phase shifts of the EIF and WB models for two levels of noise, σ, and average firing rates, ν_{0}. We found a remarkable agreement for the gain of the two models in all ranges of frequencies studied. For the phase shift, the agreement is also excellent but in a smaller range of frequencies, up to 100 Hz. Above 100 Hz, the phase lag is significantly smaller in the EIF model than in the WB model. This is a consequence of the delay δ ≈ 0.2 msec between the spikes in the EIF and WB models, which is clearly visible in Figure 3B. The corresponding additional phase lag is 2πfδ. This phase shift depends on the voltage V_{s}, which is used for the definition of the spike times in the WB model.
In Figure 10, we show how the instantaneous firing rate of WB and EIF models responds to a current step that brings the neuron from 20 to 30 Hz. Again, the figure shows an excellent agreement of EIF and WB model responses, as expected from Figures 8 and 9.
Effect of sodium channel activation kinetics on high-frequency behavior
Sodium channels have a voltage-dependent activation time constant on the order of 0.1 msec (Martina and Jonas, 1997). The WB model makes the simplifying assumption that the sodium activation is instantaneous. Because ψ(V) depends on the voltage directly in the nonlinear IF models, the effective spike initiation is also instantaneous in these models. In most situations, this assumption will not affect the results much, but there is clearly a frequency above which the finite activation time constant will become important.
To measure at which frequency the activation kinetics affects the firing rate modulation, we introduced the activation kinetics in the WB model as in the Hodgkin-Huxley formalism. For the EIF model, we replaced the spike-generating current, ψ(V), by a dynamical variable with an activation time constant, τ_{act} = 0.1 msec: τ_{act}dψ/dt = -ψ + g_{L}ΔT exp[(V - V_{T})/Δ_{T}]. After each spike, the potential was reset to V_{r}, and the spike-generating current was reset to ψ(V_{r}). Results are shown in Figure 11. In both models, the firing rate modulation was not affected at low and intermediate frequencies. It was only slightly more attenuated and had a larger phase lag at input frequency beyond several hundred Hertz. In particular, the cutoff frequency was not affected, and the instantaneous model was sufficient to describe the neuron dynamics in a wide range of input frequencies (0-500 Hz). Slower activation kinetics decreased the frequency at which the deviation from the 1/f behavior occurred. For example, for time constants of 1 msec, the response is already significantly more attenuated than 1/f below 100 Hz.
High-frequency behavior of other conductance-based neurons
So far, our results suggest that the high-frequency response depends only on the properties of the current leading to spike generation. To test this prediction, we performed numerical simulations of a conductance-based model with the same fast sodium current as in the WB model but with several additional ionic currents (Hansel and van Vreeswijk, 2002) (see Materials and Methods). These additional currents modify the subthreshold behavior of the neuron and act on time scales slower than the sodium activation time. They have a significant effect on the response at low-input frequencies because of, in particular, the adaptation current. However, Figure 12 shows that, above the cutoff frequency, the amplitude of the firing rate modulation decreases again as 1/f as in the EIF model. The dashed line in Figure 12 shows the amplitude in the asymptotic regime predicted by the EIF neuron (Eq. 15) with the same values of C and Δ_{T} as for the regularly firing and WB neurons. This confirms that high-input frequency behavior is determined only by intrinsic currents leading to spike emission and the average firing rate.
Discussion
We studied how the firing rate response of a neuron to an oscillating input depends on the modulation frequency, f. The models that we investigated include conductance-based models as well as generalized integrate-and-fire models. For sufficiently strong noise, we found that these neurons behave like a low-pass filter. For f below a cutoff frequency f_{c}, the response modulation is weakly dependent on f. For f > f_{c}, it decreases rapidly with f and behaves at large f like a power law C/f ^{α}, where C is a prefactor independent of f. We analytically calculated α and C for generalized integrate-and-fire models and have shown that α depends on the nonlinearity of spike initiation.
Of particular interest is the case of neurons with an exponential spike nonlinearity (EIF model). In this case, α = 1, and C depends on both the average firing rate of the neuron and the sharpness of spike initiation characterized by the spike slope factor, Δ_{T}. This high-frequency behavior is independent of the properties of the noise; it holds for white noise as well as colored noise and for current-based and conductance-based noise. In contrast, all other types of nonlinearity lead to exponents α different from 1.
Fast sodium currents involved in action potentials increase exponentially near firing onset. Therefore, we expect that our analytical results derived in the framework of the EIF are relevant to predicting the response of real neurons to rapidly varying inputs. This is further supported by our simulations of conductance-based models, the quantitative behavior of which is identical to that found for the EIF.
We have also shown that the cutoff frequency of the EIF neuron increases approximately linearly with the average firing rate, and that it is a decreasing function of the spike slope factor Δ_{T}. The latter dependency is very mild for white noise but much more pronounced for colored noise. We expect this result to be generic, as confirmed by our simulations of simple conductance-based models.
Implications for modeling
We have shown how a simple conductance-based neuronal model can be reduced in a systematic way to a new type of integrate-and-fire model, the EIF model, in which the dynamics involve a current that increases exponentially with the membrane potential. Remarkably, the responses to rapidly fluctuating inputs of the full and reduced models are indistinguishable. More generally, the filter characteristics are quantitatively similar for the two models in the entire range of input frequencies. This is in contrast to the LIF model, the response properties (Gerstner, 2000; Brunel et al., 2001; Fourcaud and Brunel, 2002; Mazurek and Shadlen, 2002; van Rossum et al., 2002) of which differ qualitatively from those of conductance-based models at all frequencies, as shown here. Also, although QIF and conductance-based neurons behave similarly at low frequencies, qualitative discrepancies remain at high frequencies.
The stability of asynchronous firing in large networks of interacting neurons has been studied extensively in LIF (Abbott and van Vreeswijk, 1993; Treves, 1993; Brunel and Hakim, 1999; Brunel, 2000) and QIF (Hansel and Mato, 2001, 2003) neuronal networks. The approach developed in these studies can be generalized to some extent for neurons with more realistic dynamics, provided that one knows how they respond to weakly modulated oscillating inputs. Hence, the present study paves the way for a greater understanding of the synchronization properties of conductance-based neurons. In particular, the presence of a cutoff frequency in the filtering properties of conductance-based neurons indicates that there is an upper bound to the frequency of population oscillations, which can emerge in large networks (Geisler et al., 2002).
Functional implications
Our results can be used to determine how many neurons are needed to detect a transient signal in noisy conditions, given single-neuron characteristics. Equivalently, one can calculate the minimal duration a given signal should last to be detected. Consider, for example, a population of 100 neurons emitting at a background rate of 10 Hz and subjected to a transient signal that, if applied for a time long enough, would increase the firing rate to 20 Hz. An ideal detector must be able to distinguish the transient increase of the population rate induced by the signal from rate fluctuations attributable to the finite number of neurons. This can be done using a receiver operator characteristics analysis (Green and Swets, 1966). Briefly, one computes the distributions of the spike counts in a given interval under two conditions (signal or no signal). One then computes the probability of correct detection of the signal (i.e., the probability that the spike count is larger when the signal is present than when it is absent). If neurons have very sharp spikes (Δ_{T} close to zero) in the presence of correlated noise, the signal can be detected in ∼5 msec with 90% accuracy. If the spike slope factor is Δ_{T} = 3 mV, the signal can be detected only in ∼12 msec with the same accuracy. The detection time increases further for larger values of Δ_{T}.
Experimental implications
How sharp are spikes in cortical neurons?
Our work shows that the spike slope factor, Δ_{T}, is one of the main parameters on which the response of a neuron to fluctuating inputs depends. Activation curves of Na^{+} channels have been measured in several preparations, including neocortical pyramidal cells (Fleidervish et al., 1996), hippocampal pyramidal cells, granule cells, and basket cells (Martina and Jonas, 1997; Fricker et al., 1999; Ellerkmann et al., 2001). These authors used Boltzmann functions to fit the observed data. Using their best-fit parameters, one finds Δ_{T} in the range of 3-6 mV for these types of cells. However, in all cases, there are few data points in the region of the threshold, leading to a considerable uncertainty in the estimate of this parameter. Therefore, more experiments are needed to determine the spike slope factor of cortical neurons.
Experimental measurements of the linear response of neurons
The instantaneous firing rate of neurons responding to sinusoidal currents has been measured in slice preparations (Knight, 1972b; Carandini et al., 1996; Chance, 2000). Carandini et al. (1996) and Chance (2000) found cutoff frequencies in the range of 10-100 Hz that increase with the average firing rate. This is consistent with our findings. However, none of these studies have systematically explored high-frequency behavior. Bair and Koch (1996) measured the poststimulus time histogram power spectra of MT neurons in response to visual stimuli consisting of randomly moving dots. In the examples shown in this study, the output power spectrum is rapidly attenuated in the 30-100 Hz frequency range. However, it is difficult to interpret these results, because the frequency content of the synaptic input to MT neurons is unknown. Thus, additional experiments are needed to test our prediction for the 1/f attenuation of the response at high frequencies.
Testing our experimental predictions
The main predictions of our work are: (1) the neuronal gain decays as 1/f at high frequency, independently of the characteristics of the input, (2) the cutoff frequency increases with the average firing rate, and (3) the cutoff frequency increases when the spike slope factor decreases in the presence of temporally correlated noise. One possible experimental test of these results would be to study the response of neurons in vitro, in which sodium channels underlying spike initiation are blocked with TTX and replaced, using dynamic clamp techniques by an artificial sodium current with known properties. Note that because the dynamics of sodium currents underlying spike initiation are fast, a dedicated analog circuit may be needed to emulate them with the dynamic clamp method. A similar experiment can be devised in vivo (e.g., in cortical neurons). Using QX314 (lidocaine N-ethyl bromide quaternary salt) intracellularly, sodium currents in a specific cell can be blocked without affecting the activity of its neighboring cells. This would preserve the background noise that the cell receives from the rest of the network and allow for the study of its effect on the response.
Appendix
A. Conductance-based models
The Wang-Buszáki model was introduced by Wang and Buzsáki (1996). The membrane potential V(t) is governed by the following equation: where C_{M} is the membrane capacitance (C_{M} = 1 μF/cm^{2}), I_{L} = g_{L}(V - V_{L}) is the leak current (g_{L} = 0.1 mS/cm^{2}; V_{L} = -65 mV; τ_{m} = C/g_{L} = 10 msec), I_{Na} = g_{Na}m^{3}h(V - V_{Na}) is the sodium Hodgkin-Huxley current with an instantaneous activation variable, I_{K} = g_{K}n^{4}(V - V_{K}) is the delayed rectifier potassium current, and I_{syn}(t) is the input current defined in Materials and Methods.
The dynamical equations for the gating variables are: 16 where x = m, h, n. All functions x_{∞}(V), τ_{x} are given in Table 2. The respective maximum conductance densities and the reversal potentials of the ionic currents are given in Table 3.
Another conductance-based model used in this work is taken from Hansel and van Vreeswijk (2002). The dynamics of the model are described by: 17 where the currents I_{L}, I_{Na}, and I_{K} are identical to the WB model (except V_{L} = -70 mV), but three additional ionic currents are present: (1) a persistent sodium current, I_{NaP} = g_{NaP}s_{∞}(V) (V - V_{Na}), (2) an A-type potassium current, I_{A} = g_{A}a_{∞}(V)^{3}b(V - V_{K}), and (3) a slow potassium current responsible for spike adaptation, I_{Ks} = g_{Ks}z(V - V_{K}).
The dynamical equations for the gating variables m, h, n, b, and z are given by Equation 16. All gating variable parameters are given in Table 2. The respective maximum conductance densities and the reversal potentials of the ionic currents are given in Table 3.
In most of this study, we neglect the sodium activation kinetics, and we take m = m_{∞}.
B. Analytical calculations of the neuronal response to sinusoidal modulated noisy inputs
B1. Gaussian white noise
We rewrite Equation 6 as: 18 with F(V)=-V+V_{L}+ψ(V)/g_{L} and: 19 20 where I(t) is the deterministic part of the input current, and η(t) is a Gaussian white noise with zero mean and unitary SD. The frequency of the modulation is f and ω ≡ 2πf. For convenience, we also define μ(t) = I(t)/g_{L}, which can be written as μ(t) = μ_{0} + μ_{1} cos(ωt), with μ_{0} = I_{0}/g_{L}, μ_{1} = I_{1}/g_{L}. We denote the membrane time constant τ_{m} = C/g_{L}. Equation 18 can then be rewritten as: 21 The distribution of voltages at time t obeys the Fokker-Planck equation (Risken, 1984): 22 with the boundary conditions: The reset of the voltage to V_{r} after the spike is taken into account by requiring that: 23 We solve this equation assuming that μ_{1} << μ_{0}. For the sake of simplicity, we use complex notations and write to first order: 24 where P_{1}(V, ω), ν_{1}(ω), and J_{V,1}(V, ω) are complex quantities,|P_{1}| << P_{0},|ν_{1}| << ν_{0},|J_{V,1}| << J_{V,0}. 25
Substituting Equation 24 in Equation 22 and keeping only the terms at leading order, one finds that P_{0}(V) satisfies the equation: 26 which, after one integration over V, gives: where K_{+} and K_{-} are constants. The boundary condition (Eq. 12) and reset (Eq. 23) determine K_{+} and K_{-}: 27 The solution of Equation 26 is: 28 The normalization of the distribution P_{0}(V) determines the average firing rate, ν_{0}. One finds: 29 We now consider the first order contributions in the expansions of Equation 22. One finds: 30 The boundary condition at the same order gives: 31 In the last equality, we used the fact that F(V) goes to infinity at large V.
The solution of Equations 30 and 31 simplifies in two limits:
Low-input frequency
At low frequencies (ωτ_{m} << 1), P(V, t) follows adiabatically the slowly changing input, and the left-hand side of Equation 30 can be neglected. It is then easy to check that is a solution of Equation 30. Using P_{0}(V) ≈_{V→∞} ν_{0}/F(V) (Eqs. 26, 31), we can compute the firing rate modulation in the low-input frequency limit: 32
High-input frequency
Expanding Equation 30 in powers of 1/ω and keeping only the leading order gives: 33 and ν_{1} is determined by the large V behavior of P_{1}(V) (Eq. 29). Using Equation 26 and the fact that F(V) diverges at large V, one finds: 34 and thus 35 Inserting Equation 35 into Equation 33, we find: 36 Finally, inserting Equation 36 into Equation 31, we find: 37 This formula can be applied to the EIF model to show that ν_{1} decays like 1/ω at large ω. In contrast, for any polynomial F(V), the 1/ω term vanishes, and it is necessary to calculate P_{1} and ν_{1} at the next order (i.e., 1/ω ^{2}). One finds that ν_{1} = A/ω^{2}, where the prefactor A is proportional to lim_{V→∞}(F′^{2} - FF″)/F. In particular, if F ≈ V^{2} at large V, A is finite, and ν_{1} decays like 1/ω^{2} at large ω. This is the case for the QIF model for which one finds: 38 If F ≈ V^{k} > 2 at large V, A diverges. In this case, one expects the firing rate modulation to decrease like 1/ω^{α} with 1 < α < 2.
B2. Temporally correlated noise
In the previous section, we assumed white Gaussian input noise. Here, we show that the presence of correlations does not affect the high-frequency behavior of the spike response modulation of the EIF neurons. Assuming that synaptic current decays exponentially with a synaptic time constant, τ_{s}, the equations for the dynamics of the IF neurons are: The probability density function, P(V, μ_{syn},t) (Brunel et al., 2001; Haskell et al., 2001; Fourcaud and Brunel, 2002), satisfies a two-dimensional Fokker-Planck equation: 39 The probability flux is now a two-dimensional vector with components: 40 41 In particular, in the large V limit: 42 because F(V) diverges at large V. As before, the instantaneous firing rate is given by: 43 To solve these equations, one expands P(V, μ_{syn},t), J_{V}(V, μ_{syn},t) and ν(t) in a way similar to the case of white noise. One finds: 44 and 45 for k = 0, 1. A similar analysis, as in the white noise case, shows that at first order in 1/ω: 46 Combining Equations 44-45, one can show that: 47 48 which is the same expression as in the case of white noise (Eq. 37). Therefore, the asymptotic behavior at high frequencies of the instantaneous rate modulation is not affected by the presence of temporal correlations in the noise.
Footnotes
We thank Claude Meunier for a critical reading of a previous version of this manuscript.
Correspondence should be addressed to Nicolas Brunel, Centre National de la Recherche Scientifique, Neurophysique et Physiologie du Système Moteur, Université Paris René Descartes, 45 rue des Saints Pères, 75270 Paris Cedex 06, France. E-mail: nicolas.brunel{at}biomedicale.univ-paris5.fr.
Copyright © 2003 Society for Neuroscience 0270-6474/03/2311628-13$15.00/0