Abstract
We demonstrate that singlevariable integrateandfire models can quantitatively capture the dynamics of a physiologically detailed model for fastspiking cortical neurons. Through a systematic set of approximations, we reduce the conductancebased model to 2 variants of integrateandfire models. In the first variant (nonlinear integrateandfire model), parameters depend on the instantaneous membrane potential, whereas in the second variant, they depend on the time elapsed since the last spike [Spike Response Model (SRM)]. The direct reduction links features of the simple models to biophysical features of the full conductancebased model. To quantitatively test the predictive power of the SRM and of the nonlinear integrateandfire model, we compare spike trains in the simple models to those in the full conductancebased model when the models are subjected to identical randomly fluctuating input. For random current input, the simple models reproduce 70–80 percent of the spikes in the full model (with temporal precision of ±2 ms) over a wide range of firing frequencies. For random conductance injection, up to 73 percent of spikes are coincident. We also present a technique for numerically optimizing parameters in the SRM and the nonlinear integrateandfire model based on spike trains in the full conductancebased model. This technique can be used to tune simple models to reproduce spike trains of real neurons.
INTRODUCTION
Detailed conductancebased neuron models, sometimes termed Hodgkin–Huxleytype models (Hodgkin and Huxley 1952), can reproduce electrophysiological measurements to a high degree of accuracy (for a review, see Bower and Beeman 1995; Koch and Segev 1989). Unfortunately, because of their intrinsic complexity, these models are usually difficult to analyze and are computationally expensive in numerical implementations. For this reason, simple phenomenological spiking neuron models such as integrateandfire (IF) models (Geisler and Goldberg 1966; Hill 1936; Stein 1965; Tuckwell 1988) are highly popular and have been used to discuss aspects of neural coding, memory, or network dynamics (for a review, see Gerstner and Kistler 2002; Maass and Bishop 1998).
By replacing the rich dynamics of Hodgkin–Huxleytype models by an essentially onedimensional fireandreset process, many details of the electrophysiology of neurons will be missed. In particular, standard leaky IF models do not correctly reproduce neuronal dynamics close to the firing threshold. A systematic reduction of the neuronal dynamics of type I models [i.e., neurons with a smooth frequency–current curve (Hodgkin 1948)] in the limit of very low firing rates yields a canonical type I model (Ermentrout 1996; Ermentrout and Kopell 1986; Hoppensteadt and Izhikevich 1997) that is equivalent to a quadratic IF model (Hansel and Mato 2001; Latham et al. 2000). Although quadratic IF models give, by construction, a correct description of neuronal dynamics close to the firing threshold, it is unclear how well quadratic and other generalized IF models such as the exponential IF model (FourcaudTrocmé et al. 2003) or the Spike Response Model (SRM) (Gerstner and Kistler 2002; Kistler et al. 1997) perform for a realistic, timedependent input scenario where the neuron could spend a significant amount of time far away from the firing threshold. Keat and colleagues (2001) have shown that a phenomenological model of neuronal activity can predict every spike of lateral geniculate nucleus (LGN) neurons with a millisecond precision. However, these neurons produce very stereotyped spike trains with short periods of intense activity followed by long periods of silence (Reinagel and Reid 2002) so that their approach could be limited to LGN neurons only. Recently, it was shown that an IF model can predict the mean rate of pyramidal cells recorded in in vitro experiments (Rauch et al. 2003) over a broad range of different timedependent inputs. Moreover, the experimental distributions of membrane potentials in the subthreshold regime are well reproduced by leaky IF models (Destexhe et al. 2001). Quadratic IF neurons can approximate the frequency–current curve of a detailed conductancebased model (Hansel and Mato 2003). Finally, it was shown that generalized IF models can approximate the dynamics of the classic Hodgkin–Huxley model of the squid giant axon with high accuracy (Feng 2001; Kistler et al. 1997). The aim of the present paper is 2fold.
First, we attempt to illustrate the relation of phenomenological spiking neuron models to conductancebased models (Abbott and Kepler 1990; Destexhe 1997; Ermentrout 1996; Ermentrout and Kopell 1986; Kistler et al. 1997; Latham et al. 2000). To do so, we will use a stepbystep analytical derivation of 2 formal spiking neuron models starting from a detailed conductancebased model of a fastspiking cortical interneuron. In both cases, we compare the behavior of the reduced model to that of the detailed conductancebased model for 3 different input scenarios: a strong isolated current pulse, constant superthreshold input, and random current. Although isolated current pulses and constant drive are standard experimental paradigms, a random current is thought to be more realistic in that it reflects an approximation to the random conductance background caused by input spikes of cortical neurons in vivo (see, e.g., Calvin and Stevens 1968; Destexhe and Paré 1999).
Second, we apply a numerical technique to map the reduced models to the detailed conductancebased model of an interneuron. The optimal parameters characterizing the reduced models are extracted from a sample spike train generated with the conductancebased neuron model by a procedure which generalizes previous approaches (Brillinger 1988; Brillinger and Segundo 1979; Wiener 1958). We compare quantitatively the predictions of the reduced models with that of the full conductancebased model in the random input scenario. We explore the regime of validity of the reduced models by varying the mean and the standard deviation of the input current in a biologically realistic range, and we test the robustness of the reduced models to further simplifications. A generalization of our approach to adapting neurons is indicated. Finally, we extend our numerical technique to a conductance injection scenario, which has become increasingly popular in recent years (Destexhe and Paré 1999; Destexhe et al. 2001; Harsch and Robinson 2000; Reyes 2003).
METHODS
Integrateandfire models
In formal spiking neuron models, action potentials are generated by a threshold process. The neuron fires whenever the variable u reaches a threshold ϑ from below (1) where t̂ is called the firing time of the neuron. We focus in this article on models that are fully described by a single variable u. Wellknown idealized spiking neuron models can differ from one another in the specific way that the dynamics of the variable u are defined. In the standard leaky integrateandfire model (LIF model), the evolution of u is given by a linear differential equation (2) If we identify u with the membrane potential of the neuron, and I^{ext} with the external driving current, we may interpret τ_{m} as the membrane time constant, R as the input resistance, and u_{eq} as the equilibrium potential of the leakage conductance. Integration of Eq. 2 yields the membrane potential as a function of time. Firing is defined by the threshold condition (Eq. 1). After firing, the membrane potential is reset to a value u_{reset}. The LIF neuron may also incorporate an absolute refractory period, in which case we proceed as follows. If u reaches the threshold at time t = t̂, we interrupt the dynamics (Eq. 2) during an absolute refractory time δ_{refr} and restart the integration at time t̂ + δ_{refr} with the new initial condition u_{reset}. We emphasize that firing is a formal event defined by the threshold condition (Eq. 1). The time course of the action potential is disregarded in standard IF models and only the firing time t̂ is recorded. In a general nonlinear integrateandfire model (NLIF model), Eq. 2 is replaced by (Abbott and van Vreeswijk 1993) (3) As before, the dynamics is stopped if u reaches the threshold ϑ and reinitialized at u = u_{reset}. A specific instance of a NLIF model is the quadratic model (Feng 2001; Hansel and Mato 2001; Latham et al. 2000) (4) with parameters a_{0} > 0 and u_{eq} < u_{c} < ϑ. For I^{ext} = 0 and initial conditions u < u_{c}, the voltage decays to the resting potential u_{eq}. For u > u_{c}, the voltage increases up to ϑ where an action potential is triggered; u_{c} can therefore be interpreted as the critical voltage for spike initiation by a short current pulse. The quadratic IF model is closely related to the Θneuron, a canonical type I neuron model (Ermentrout 1996; Hoppensteadt and Izhikevich 1997).
Spike response model
In contrast to the NLIF model, the LIF model defined in Eq. 2 can be analytically integrated for arbitrary timedependent input I^{ext}(t). Let us denote the last firing time by t̂. For t > t̂ + δ_{refr}, the membrane potential is (5) where H(x) is the Heaviside step function defined by H(x) = 1 for x ≥ 0 and zero otherwise; u_{reset} is the initial condition of the integration at time t̂ + δ_{refr}. If we replace the exponentials multiplied by a step function by “response kernels” η and κ, we may rewrite Eq. 5 in the form (6) Equation 6 is a generalization of Eq. 5 and has been termed the Spike Response Model (SRM; Gerstner 1995; Kistler et al. 1997). The kernel η models, in the general case, the spike itself and the afterhyperpolarization that follows the spike. In the specific case of Eq. 5, η describes the hard reset of u(t) to u_{reset} at time t = t̂ + δ_{refr}. The kernel κ describes the response of the membrane potential to an input current pulse. In case of the LIF neuron, both kernels are characterized by an exponential decay with time constant τ_{m} (compare Eq. 6 and 5).
Although both the NLIF model (Eq. 3) and the SRM (Eq. 6) contain the LIF model (Eq. 2) as a special case, the direction of the generalizations is somewhat different. In the NLIF model, parameters are made voltage dependent, whereas in the SRM they depend on t − t̂ (i.e., the time since the last spike). To illustrate the relation, compare the NLIF models in Eq. 3 and 4 to an IF model with a timedependent time constant (i.e., a SRM; Stevens and Zador 1998; Wehmeier et al. 1989) (7) with C = τ_{m}/R. Starting the integration at u(t̂ + δ_{refr}) = u_{reset}, we find for t > t̂ + δ_{refr} (8) which is a special case of Eq. 6. A timedependent time constant takes care of the fact that during and shortly after a spike, many ion channels are open so that the resistance of the membrane and hence its time constant is reduced.
As a further generalization, we replace in the SRM the fixed threshold by a dynamic one (Fuortes and Mantegazzini 1962; Geisler and Goldberg 1966; Holden 1976; Stein 1967; Weiss 1966) (9) During an absolute refractory period δ_{refr}, we may, for example, set ϑ to a large and positive value to avoid firing and let it relax back to its equilibrium value for t > t̂ + δ_{refr}.
If input is provided by the activation of a conductancebased synapse (e.g., attributed to presynaptic spike arrival), rather than current injection, a different formulation of the state of the neuron is provided. The convolution with the response kernel κ is replaced by a sum and a response kernel ε that takes into account the filtering ascribed to the synapses (10) where the sum runs over all firing times f of all presynaptic neurons j with w_{j} the strength of synapse j. For the sake of simplicity, we will assume throughout this article that all synapses have the same strength w_{j} = 1. The kernel ε can be interpreted as the postsynaptic potential generated by each spike and it also depends on the time elapsed since the last emitted spike (see above). This latter equation can be transformed back into a convolution product. Let us define the sequence of presynaptic spikes by S(t) = Σ_{j} Σ_{f} δ(t − t_{j}^{(f)}). Equation 10 can thus be restated (11) This latter formulation will be used for the numerical fitting of the SRM to the full conductancebased model in case of random spike arrival.
Full conductancebased model
As our reference model, we take the conductancebased model of a fastspiking cortical interneuron proposed by Erisir and coworkers (1999). We have chosen this specific model so that, because fastspiking neurons show little adaptation, we avoid most of the complications caused by slow ionic processes that would not be well captured by the class of idealized spiking neuron models reviewed above. In fact, because the model of a fastspiking neuron is comparatively simple, we can hope to illustrate the steps necessary for a reduction to spiking models in a transparent fashion.
The fastspiking model neuron (Erisir et al. 1999) follows the Hodgkin–Huxley formalism and consists of a single homogeneous compartment with a nonspecific leak current and 3 active ionic membrane currents (12) The sodium current I_{Na} is described by the equation (13) with two gating variables m and h. A slow potassium current is modeled by (14) with gating variable n_{1}. Because of its slow time constant, the variable n_{1} builds up over several spikes and contributes to a subtle form of adaptation. There is a further, much faster potassium current (15) with a conductivity g_{K2} that is much larger than g_{K1}. It is a strong outward current associated with the Kv3.1–Kv3.2 channel that shapes the downstroke of the action potential and contributes to a rapid reset of the membrane potential after a spike (Erisir et al. 1999). Finally, there is the linear leak current (16) which defines a passive membrane time constant τ_{m} = C/g_{I}. Each gating variable follows the equation (17) where x stands for m,h,n_{1}, or n_{2}. We note that both potassium channels are essentially closed at rest. The parameters of all currents (see Table 1) were chosen by Erisir et al. (1999) so as to reproduce the behavior of fastspiking neocortical interneurons. We adopt their set of parameters, except for the leak current where we have chosen g_{l} = 0.25 mS/cm^{2} (instead of their value of 1.25 mS/cm^{2}) so as to set the passive membrane time constant to about 4 ms. This larger value of the passive membrane time constant seems to be realistic for a broad class of interneurons (Wang et al. 2002). With this set of parameters, the Erisir model exhibits the typical behavior of a type I neuron model, that is, a continuous frequency–current curve that allows firing at very low frequencies and delayed action potentials for pulse input that is just superthreshold.
The model was simulated on a Sun Workstation using the forward Euler integration method with a time step of 0.01 ms. Because the time step is much shorter than all intrinsic time constants of the model (see Table 1), such a simple firstorder method is sufficient.
Current injection
In the current injection scenario, the neuron is driven with an uncorrelated Gaussiandistributed random current. We vary both the mean μ and the standard deviation σ of the Gaussian current. This kind of highly variable timedependent input is thought to be more realistic than other current injection stimulation protocols (Calvin and Stevens 1968; Destexhe and Paré 1999).
Conductance injection (stochastic presynaptic spike arrival)
An even more realistic input scenario is to consider stochastic spike arrival at excitatory and inhibitory synapses. Such a scenario is equivalent to driving the neuron with a highly variable random conductance multiplied by a driving force that depends on the instantaneous membrane voltage of the cell. More precisely, the total synaptic input current I_{syn} is given by (Robinson and Kawai 1993) (18) where g_{exc} (g_{inh}) is the total excitatory (inhibitory) conductance and E_{exc} (E_{inh}) is the corresponding reversal potential. The synaptic conductances g_{exc} and g_{inh} consist of the summed input from N_{exc} excitatory synapses and N_{inh} inhibitory synapses, respectively (19) (20) The dynamics of the jth excitatory synapse is described by a variable with (Dayan and Abbott 2001) (21) The value of increases by an amount A^{exc} for each presynaptic spike arriving at the synapse at time t_{j,k}, and decays with a time constant τ_{exc}. The result is a lowpass–filtered version of the presynaptic spike train {t_{j,k}}. The dynamics of inhibitory synapses are defined similarly.
Presynaptic spike trains are described by random homogeneous Poisson processes. At each time step, N̅ independent random variables are generated and distributed among the N > N̅ synapses to generate slightly correlated spike trains (Destexhe and Paré 1999). Numerical values are summarized in Table 2.
Coincidence factor Γ
The coincidence factor Γ between two spike trains (Kistler et al. 1997) is defined by (22) where N_{1} is the number of spikes in the reference spike train S_{1}, N_{2} is the number of spikes in the spike train S_{2} that is compared to the reference spike train, N_{coinc} is the number of coincidences with precision Δ between the 2 spike trains, and 〈N_{coinc}〉 = 2νΔN_{1} is the expected number of coincidences generated by a homogeneous Poisson process with the same rate ν as the spike train S_{2}. In this article, the reference spike train S_{1} is always generated by the full conductancebased model of a fastspiking interneuron, whereas S_{2} is a spike train generated by one of the reduced models (IF models). The factor N = 1 − 2νΔ normalizes Γ to a maximum value of one, which is reached if and only if the spike train of the reduced model reproduces exactly that of the full model. A homogeneous Poisson process with the same number of spikes as the reduced model model would yield Γ = 0. We compute the coincidence factor Γ by comparing the 2 complete spike trains: the spike train S_{1} generated by the full conductancebased model and the train S_{2} predicted by the reduced model.Therefore, in this article, Γ gives a measure of the ability of the reduced model to predict the spike train of the full model.
Analytical reduction to an IF model
To find a NLIF model that approximates the dynamics of the full fastspiking neuron model, we proceed in 2 steps. In a first step, we keep all variables, but introduce a threshold for spike initiation. We call this the multicurrent IF (MCIF) model. In a second step, we separate gating variables into fast and slow ones (Abbott and Kepler 1990; Kepler et al. 1992; Rinzel 1985). The fast variables are turned into instantaneous ones, whereas the slow variables are replaced by constants. The result is the desired NLIF model, which depends on only a single variable.
step 1.
For the first step, we make use of the observation that the shape of an action potential of the fastspiking neuron model is always roughly the same, independently of the way the spike is initiated. Instead of calculating the shape of an action potential again and again, we can therefore simply stop the costly numerical integration of the nonlinear differential equations as soon as a spike is triggered and restart the integration after the downstroke of the spike about 1.5–2 ms later. We call such a scheme a multicurrent IF model. The interval between the spike trigger time t̂ and the restart of the integration corresponds to an absolute refractory period δ_{refr}.
The MCIF model is defined by a voltage threshold ϑ, a refractory time δ_{refr}, and the reset values from which the integration is restarted. We fix the threshold at ϑ = −40 mV; the exact value is not critical, and we could take values of −20 to −45 mV without significantly changing the results. With ϑ = −40 mV, a refractory time of δ_{refr} = 1.7 ms is suitable and an appropriate reset value for the voltage variable is u_{reset} = −85 mV.
The first important approximation involves the reset of the gating variables m, h, n_{1}, and n_{2}, because their time courses are not as stereotyped as that of the membrane potential. To illustrate this point, let us focus on the variable h. At 1.7 ms after spike initiation, the variable h is at a value of 0.16 during periodic firing at about 40 Hz; however, it is at a value of 0.26 when only a single action potential is triggered (not shown). Thus, the optimal value of h_{reset} depends on the choice of input scenario. If the biological neuron under consideration has a mean firing rate of 40 Hz, then we should choose a reset value appropriate for the particular regime of periodic firing. If, on the other hand, the biological neuron is near rest most of the time and emits only an occasional spike, we should choose different reset values. In the following, we adjust the reset values based on a scenario with constant drive current I^{ext} = 5 μA/cm^{2} that leads to a mean firing rate of about 40 Hz. The reset values are m_{reset} = 0.0; h_{reset} = 0.16; n_{1,reset} = 0.874; n_{2,reset} = 0.2; and u_{reset} = −85 mV. This set of parameters leads to a nearperfect fit of the time course of the membrane potential during periodic firing at 40 Hz and approximates the gain function of the full fastspiking neuron model to a high degree of accuracy (Fig. 1, A and B).
For stimulation with isolated 2ms current pulses, the MCIF model reproduces both the threshold behavior and the firing times of the full fastspiking neuron model. The hyperpolarizing spikeafterpotential, however, shows significant differences (see Fig. 1C). If we took reset values adapted to pulse input from rest, approximation would be more precise. We will, however, stick to the previous set of reset parameters and keep it fixed throughout the rest of this section.
We now turn to random input. The amplitude σ_{I} of the fluctuations determines the mean firing rate of the neuron model. We see from Fig. 1D that the mean rate of the MCIF model follows closely that of the full model. For a more detailed comparison of the MCIF with the full model, we stimulate both models with exactly the same random current (i.e., identical initiation of the random number generators). In Fig. 1E, we see that the voltage time course of the 2 models is indistinguishable most of the time. Occasionally, the MCIF model misses a spike or adds an extra spike. For this specific input scenario (where the input fluctuations have strength 25 μA/cm^{2}), about 95% of the spike times or the full model are reproduced correctly (with a resolution of Δ = ±2 ms) by the MCIF model. We see from Fig. 1F that the coincidence rate Γ varies as a function of the fluctuation amplitude, but stays above 0.75 in the whole range that we considered. As expected, the value of Γ is highest when the neuron fires at an average rate between 35 and 55 Hz (i.e., in the regime for which parameters have been optimized).
step 2.
The above MCIF neuron is not yet the desired NLIF model because it still depends on all 5 variables, u, m, h, n_{1}, and n_{2}. To reduce it further to a singlevariable model that depends only on the membrane potential u, we need to consider a gating variable x either as fast compared to u, in which case we replace x by x_{∞}(u), or as slow compared to u, in which case we take x as constant (see also Abbott and Kepler 1990; Kepler et al. 1992; Rinzel 1985). In our case, m(t) is the only fast variable (in the subthreshold range [−100, −40] mV, we find 0.08 ≤ τ_{m} ≤ 0.25 ms, 4.28 ≤ τ_{h} ≤ 14.45 ms, 44.66 ≤ τ_{n1} ≤ 144.10 ms and 0.44 ≤ τ_{n2} ≤ 4.19 ms). We eliminate m by replacing m(t) by its equilibrium value m_{0}[u(t)] (see Fig. 2A). The treatment of the other gating variables deserves some extra discussion.
We start with the variable n_{2}. A thorough inspection of the time course of n_{2}(t) shows that n_{2} is close to its resting value most of the time, except for a short interval during and immediately after the downstroke of an action potential (see Fig. 2B). If we take a refractory time of δ_{refr} = 4 ms, most of the excursion trajectory of n_{2} falls within the refractory period. Between spikes we can therefore replace n_{2} by its equilibrium value at rest n_{2} ≡ n_{2,eq} = n_{0}(u_{eq}).
The gating variables h and n_{1} vary only slowly so that, for a given input scenario, the variables may be replaced by their average value h_{av} and n_{1,av}. The average, however, does depend on the input scenario. To illustrate this point, we consider 3 different situations. First, for a neuron at rest that receives only occasionally an isolated current pulse, the average values are close to the resting values. In other words, we may use the parameters at rest (i.e., h_{av} = 0.87 and n_{1,av} = 0.00057). If the neuron is stimulated by a constant bias current that shifts the membrane potential just below threshold, we should take h_{av} = 0.54 and n_{1,av} = 0.019. Finally, in a periodic firing regime (constant stimulation with I = 5 μA/cm^{2}) reasonable averages are h_{av} = 0.45 and n_{1,av} = 0.8 (see Fig. 2, C and D). The last pair of values will be used in the following as our standard set of parameters unless stated otherwise.
With m = m_{0}(u) and constant values for h, n_{1}, and n_{2}, the dynamics of the full fastspiking neuron model defined in Eqs. 12–17 reduces to (23) After division by C, we arrive at a single nonlinear equation (24) The function F depends on the choice of constants h_{av} and n_{1,av}, which can, as indicated above, be optimized differently for each of the 3 input scenarios discussed. In particular, the zero crossings F(u) = 0 that define the resting potential u_{eq} and the critical voltage u_{c} for spike initiation are shifted in a model adapted to periodic firing compared to a model adapted to rest. By definition, the passive membrane time constant of the model is inversely proportional to the slope of F at rest: τ = −1/(dF/du) where the derivative is to be evaluated at u = u_{eq}. In principle the function F could be further approximated by a linear function with slope −1/τ and then combined with a threshold at, say, ϑ = −45 mV. This would yield a standard linear LIF model. Alternatively, F could be approximated by a quadratic function that would lead us to the quadratic IF neuron. The canonical type I model would correspond to a quadratic IF neuron where parameters are optimized in the limit of a very low firing frequency (Ermentrout 1996; Ermentrout and Kopell 1986; Hansel and Mato 2001; Latham et al. 2000).
In the following, we do not make any further approximations. Instead we work directly with the nonlinear function F(u). The refractory period δ_{refr} = 4 has already been introduced above. To finish the definition of the model, we have to specify the threshold ϑ and the reset value. We take ϑ = −45 mV (the exact value is not critical). After firing, integration is restarted at u_{reset} = −85 mV.
Analytical reduction to a spike response model
For a fitting of the full fastspiking neuron model defined in Eqs. 12–17 to the SRM defined in Eqs. 6 and 9, we need to determine the kernels η(t − t̂) and κ(t − t̂, s), and furthermore the (timedependent) threshold ϑ(t − t̂) must be adjusted. As a first step, we stimulate the full model by a short suprathreshold current pulse to determine the time course of the action potential and afterhyperpolarization. Let us define t̂ as the time when the membrane potential crosses an (arbitrarily set) threshold ϑ (e.g., ϑ = −50 mV). The kernel η(t − t̂) is defined by the time course of the membrane potential for t > t̂ (i.e. during and after the action potential) in the absence of external input. If we were interested in a purely phenomenological model, we could simply record the numerical time course u(t) and define η(t − t̂) = u(t̂) for t > t̂ (see Numerical optimization method). Instead of such a purely numerical method, it is instructive to take a semianalytical approach and study the 4 gating variables m, h, n_{1}, and n_{2} immediately after action potential generation. About 2 ms after initiation of a spike, all 4 variables have passed their maximum or minimal values and are on their way back to equilibrium (cf. Fig. 3, A and B). We set δ_{refr} = 2 ms. For t ≥ t̂ + δ_{refr}, similar to the approach of Destexhe (1997), we fit the approach to equilibrium by an exponential (25) where x = m, h, n_{1}, n_{2} stands for each of the 4 gating variables. The parameter τ_{x} is a fixed time constant and x_{reset} is the initial condition at t = t̂ + δ_{refr}, and x_{eq} = x_{∞} (u_{eq}). In other words, the voltagedependent differential Eq. 17 for x is replaced by the linear voltageindependent equation (26) Numerical values for x_{eq}, τ_{x}, and x_{reset} are summarized in Table 3. We note that, for the gating variable m, the time constant is not related to the gating dynamics but reflects the time course of the membrane potential u toward the resting potential.
With the time course of the gating variables defined in Eq. 25, we know the conductance of each ion channel as a function of time. For example, the potassium current I_{K2} is (27) where g_{2}(t − t̂) is an exponential function with time constant τ_{n2}/2. We insert the timedependent conductance into Eq. 12 and find for t ≥ t̂ + δ_{refr} (28) where the sum runs over the 4 ion channels I_{Na}, I_{K1}, I_{K2}, and I_{1}. With the definitions τ(t − t̂) = C/Σ_{j} g_{j}(t − t̂) and I^{ion}(t − t̂) = Σ_{j}g_{j}(t − t̂) E_{j}, we arrive at (29) which is a linear differential equation with timedependent time constant (cf. Eq. 7). The effective time constant is shown in Fig. 3C as a function of t − t̂.
According to Eqs. 27 and 28, the ion currents have, at time t̂ + δ_{refr}, a value of I^{ion}(δ_{refr}) = Σ_{j} g_{j}(δ_{refr})E_{j}. For t − t̂ = δ_{refr}, the effective time constant τ(δ_{refr}) of the membrane voltage is extremely short (Fig. 3C) so that we may assume that it is at its steady state value given I^{ion}(δ_{refr}). For t ≥ t̂ + δ_{refr}, we therefore integrate Eq. 29 with the initial condition (30) This yields the SRM in the form of Eq. 6 for t > t̂ − δ_{refr} with the kernels given by (31) (32)
The kernels κ and η that we have constructed so far are limited to t > t̂ + δ_{refr}. During the absolute refractory period t̂ < t < t̂ + δ_{refr} the neuron is not responsive to external input. We therefore set κ to zero. The action potential itself that occurs in the interval t̂ < t < t̂ + δ_{refr} is, for the sake of simplicity, approximated by a triangular voltage pulse. Finally, we introduce a dynamic threshold (33) with constants ϑ_{0}, ϑ_{1}, ϑ^{refr}, and τ_{ϑ}. During the absolute refractory period δ_{refr}, the threshold is set to a value ϑ^{refr} = 100 mV that is sufficiently high to prevent the neuron from firing. After refractoriness, the threshold starts at ϑ_{0} + ϑ_{1} and relaxes with a time constant of τ_{ϑ} to an asymptotic value of ϑ_{0}. The initial value ϑ_{0} + ϑ_{1} was arbitrarily set at 0 mV. ϑ_{0} = −50 mV and τ_{ϑ} = 6 ms were then chosen so that for 2 different inputs (i.e., I_{0} = 5 μA/cm^{2} and I_{0} = 30 μA/cm^{2}) the mean firing rate of the SRM was approximately correct.
Numerical optimization method
In the SRM framework, 3 quantities are needed to define the model. These are the kernel η, describing the shape of a spike; the kernel κ, describing the input integration process and finally the threshold ϑ for spike initiation. The previous section described a semianalytical method to reduce the full conductancebased model to a SRM, obtaining approximations of the kernels. In this section, we describe a numerical method that finds kernels that optimize the fit of the SRM to a spike train. The method proceeds in 3 steps. First, we use the fact that spikes have always roughly the same shape to extract η from the course of the membrane voltage. Second, we use a linear approximation for the integration of the input to extract κ as the best linear filter between the driving current and the fluctuations of the membrane potential in the subthreshold regime. Third, we fit the parameters of the threshold ϑ to find the best reproduction of the specific spike train that we use. The resulting fit is finally evaluated on a set of new spike trains that have not been used during optimization. We now discuss the 3 steps in more detail.
Throughout this subsection, we suppose that we have at disposal a simultaneous recording of the discretized membrane voltage u_{t}^{data} and the discretized driving current I_{t}^{ext} of a cell. We assume that both u_{t}^{data} and I_{t}^{ext} are recorded at the same sampling frequency.
1) To extract the discretized version of the kernel η from the data, we align all the spikes. The spike triggered average of the membrane potential (i.e., the “mean shape” of an action potential) yields, after smoothing, the kernel η (see Fig. 4, B and C). Detection and alignment are realized relative to an arbitrary threshold condition on the first time derivative of the membrane voltage (see Fig. 4A).
2) After the extraction of η, we move to the extraction of κ. We start by discretizing the SRM equations. The equivalent of Eq. 6 in discrete time is (34) We suppose that the driving current is time dependent and has the form (35) where μ_{1} and σ_{1} are constants and ξ_{1} are independent random variables drawn from a normal distribution with mean 0 and unit variance. For the sake of simplicity, we set μ_{I} = 0 for the rest of this discussion. However, the method can easily be generalized. Equation 34 can therefore be rewritten (36) It is useful to keep in mind that the version of the SRM described in this article is a memoryless model (i.e., nothing is remembered of what happened before the last spike). By shifting the ηterm to the lefthand side, we subtract the effects of the action potentials and arrive at an expression for the subthreshold regime, which we denote as V_{t−t̂,t}^{SRM} (37) Similarly, for the data u_{t}^{data}, we subtract the time course of the average action potential that yields a subthreshold behavior V_{t−t̂,t}^{data} (38)
The righthand side of Eq. 37 can be interpreted as a numerical convolution between σ_{I}ξ_{t} and a family of filters parameterized by t − t̂. The problem of finding the best linear filter (κ_{t−t̂,s}) of some device given the output vector (V_{t−t̂,s}^{data}) and the input vector (σ_{I}ξ_{t}) is known in the literature as the Wiener–Hopf optimal filtering problem. However, the classic formulation of this problem is not well suited for the present question because of the explicit dependency of the filter on t − t̂. In the appendix, we show that the Wiener–Hopf approach can nevertheless be adapted to the present situation. The formulation we propose circumvents a couple of technical issues of the classic formulation. In particular, there is no need to choose a specific time window in which to perform the extraction of the filter. Instead, all segments are aligned for a given t − t̂. We find that κ_{t−t̂,s} is a solution of the following linear system (39) where s_{max} is the maximum time lag for which κ_{δt,s} is nonzero, δt = t − t̂, and X[·, ·] is defined by (40) where T is the number of spikes in the spike train and f and g denote arbitrary vectors. X[·, ·] can be interpreted as a generalization of the empirical crosscorrelation between f and g. It is therefore possible to find an expression for the family of filters κ_{δt,s} for each δt and s by solving Eq. 39. To smooth the results, we fit them with a suitable function. In the present article, we use a single exponential decay in the variable s (see Fig. 4D). It is then possible to plot the time constant τ of this function versus the delay δt and fit the dependency on δt by the function τ(δt) = α_{1}(1 + tanh (α_{2}(δt − α_{3})]) with free parameters α_{1}, α_{2}, and α_{3}. The SRM can then be turned back into a differential equation with a timedependent time constant (Eq. 7; see also discussion), which is very efficient for simulations.
3) Finally, we choose a specific threshold condition and optimize it in terms of spike train reproduction. We take a dynamic threshold defined by Eq. 33, where δ_{refr} is always kept at 2 ms and ϑ^{refr} at a value of 100 mV to prevent firing. ϑ_{0}, ϑ_{1}, and τ_{θ} are free parameters. To find the best fit for these 3 parameters, we simulate the full spike train with the SRM and compute 1 − Γ. For a definition of Γ, see Eq. 22. Then, we use the downhill simplex method (Nelder and Mead 1965) algorithm to find a set of parameters that minimize 1 − Γ.
For this article, we have used for the parameter optimization a single spike train of 10 s length generated with Gaussiandistributed random current (cf. Fig. 4F). Parameters of the current are μ_{I} = 0 μA/cm^{2} and σ_{1} = 25 μA/cm^{2}. The value of the driving current is changed every 0.2 ms. Both the membrane voltage and the driving current are sampled at 5kHz sampling frequency. For this set of parameters, the mean rate is 34.6 Hz.
In case of conductance injection protocol (stochastic spike arrival), steps 1 and 3 of the method are left unchanged, whereas step 2 is slightly modified. The state of the neuron in discrete time is given by (see Eq. 10) (41) with an εkernel for each type of synapses (excitatory and inhibitory synapses). Using the same stratagems as before (see Eqs. 11 and 37 and Fig. 9), we find (42) where S_{t}^{exc} and S_{t}^{inh} are the histograms of spike arrivals at excitatory and inhibitory synapses, respectively. More precisely, a value S_{t}^{exc} = n means that in a time step of duration 0.2 ms, a total of n excitatory synapses received spike inputs. At this stage, a formalism very similar to the one applied before can be used to find the best linear filters ε^{exc} and ε^{inh} both at the same time (see appendix). The filters ε^{exc} and ε^{inh} are fitted by double exponentials (see Fig. 9). The kernel η was maintained from current injection simulations. The rest of parameter optimization was done on a single spike train of 1 s length generated with discharge frequencies ν_{exc} = 0.3 Hz and ν_{inh} = 2.5 Hz. Both the membrane voltage and the driving current are sampled at 5 kHz sampling frequency. In this specific spike train, the mean rate is 29.0 Hz. Finally, S^{exc} (S^{inh}) is the poststimulus time histogram PSTH computed over all excitatory (inhibitory) synapses with a given time binning (here, each bin equals 0.2 ms). See also Eq. 11.
RESULTS
We have compared a detailed conductancebased model of a fastspiking neuron (Erisir et al. 1999), in the following called “full model” with several versions of generalized IF models (also referred to as “reduced models”). The details of the full model and the reduced models can be found in the methods section. The results of our work are organized as follows.
We first present results with a nonlinear integrateandfire (NLIF) model that has been optimized by the analytical and numerical methods explained in the methods section. A slightly different analytical approach leads to the Spike Response Model (SRM) that we discuss in the second subsection. For both the NLIF and the SRM, the quality of the reduced model is assessed for 3 different input paradigms: constant current of variable strength, pulse input of 2ms duration and variable amplitude; and random input current of adaptable variance. Extensions of the model to the case of random conductance input or adaptation are discussed separately.
Comparison of the NLIF model and the full model
We have analytically approximated (see methods) the 5 equations of the full model of a fastspiking neuron (Erisir et al. 1999) by a NLIF with a single variable u, given by the differential equation (43) The NLIF neuron fires when u hits a threshold of ϑ = −45 mV. Integration then restarts after a refractory period of 4 ms at u_{reset} = −85 mV. The function F(u) is shown in Fig. 5F and parameterized by 3 constants, h_{av}, n_{1,av}, and n_{2,eq}. When the analytical reduction is based on the assumption that the neuron is close to the resting state, F(u) is given by the solid line; if we assume that the neuron is close to but below threshold, F(u) is given by the shortdashed line; if we assume that the neuron fires repetitively at about 40 Hz, then the analytical reduction yields the function F(u) indicated by the longdashed line. Thus although it is always possible to approximate the dynamics of the full model by a onedimensional NLIF, the exact form of the function F depends on the regime for which the NLIF is optimized. To emphasize this fact, we now compare the behavior of the NLIF to that of the full fastspiking neuron model for 3 different input scenarios.
Let us focus on constant input first. With the function F that has been derived for repetitive firing at 40 Hz, the membrane trajectory during periodic firing at that frequency is well approximated (cf. Fig. 5A). The firing frequency as a function of input current (gain function) is reasonably well approximated at frequencies above 40 Hz, but not at lower frequencies. In particular, the threshold for periodic firing is not reproduced correctly (cf. Fig. 5B). If we take the function F(u) that is appropriate at threshold, the behavior at very low firing rates is well approximated, but the gain function increases rapidly to unplausible values.
In our scenario with pulse input, the neuron is at rest just before a current pulse of 2ms duration arrives. Thus the best choice of the function F(u) is the one that has been derived for the resting state. Indeed, with the function F(u) that is adapted to a neuron at rest, the NLIF model generates spikes whenever the amplitude of the 2ms pulse is above 8.5 μA/cm^{2}, which is close to the pulse initation threshold of 8.8 μA/cm^{2} of the full model. Moreover, for pulse input that is just superthreshold, the NLIF model responds with delayed action potential initiation, just as the full model (data not shown). If we take, however, the NLIF model with the function F that we found optimal for periodic firing at 40 Hz, the NLIF generates an action potential only for strong input pulses with amplitude I^{ext} ≥ 15 μA/cm^{2}, whereas the full model triggers spikes for I^{ext} ≥ 8.8 μA/cm^{2}. Moreover, a difference in the resting potential is evident in a comparison of the 2 models (cf. Fig. 5C). The resting potential of the NLIF model is given by the leftmost zero crossing of the function F(u) plotted in Fig. 5F. It is reduced by a value of about 5 mV compared to that of the full model.
In summary, the 2 previous artificial stimulation paradigms (i.e., constant input and pulse input) have shown that the choice of the funtion F is critical for a good approximation of the dynamics of the full model by the NLIF model. We now turn to stimulation by random input current, which we consider the most realistic paradigm, given that neurons in vivo are thought to receive highly fluctuating synaptic current (Calvin and Stevens 1968; Destexhe and Paré 1999). To choose a function F that yields a good performance of the NLIF model for random input, we take advantage of the fact that F is characterized by 3 constants h_{av}, n_{1,av}, and n_{2,eq} (see methods). We take these constants and the threshold as the 4 free parameters that we optimize by a standard numerical procedure (simplex algorithm; Nelder and Mead 1965) using a fluctuating input current of zero mean and variance σ_{I} = 25 μA/cm^{2} during 10 s of stimulation. Ideally, the 2 neuron models (i.e., full and reduced model) should produce exactly the same spike train if stimulated by the same timedependent input current (that is, the same realization of a random current). The quality of the reduced model can be assessed by a comparison on the level of individual spikes or on the coarser level of firing rates. We choose to optimize the 4 parameters by comparing the spike train of the full model and that of the reduced model on a spikebyspike basis. In particular, we maximize the coincidence rate Γ that measures the similarity of the 2 spike trains with a temporal resolution of 2 ms (see methods). The parameters found by the numerical optimization are given in data supplements, Table 1.A sample spike train is shown in Fig. 5D. Most of the spikes coincide, whereas some spikes are missed and other spikes are added. Action potentials of the NLIF model are indicated by triangular pulses that span a refractory period of δ_{refr} = 2 ms; the performance does not change significantly if we take δ_{refr} in the range of 1 to 4 ms.
We now turn to a systematic exploration of the performances of the numerically optimized NLIF for fluctuating input with different mean μ_{I} and variance σ_{I}. In Fig. 5E we plot the mean firing rate as a function of the fluctuation amplitude for both the full model and the NLIF model. The NLIF model gives a fair reproduction of the mean rates of the full model except for large variance [σ_{I} > 50 (μA/cm^{2})] of the input current. For a more detailed comparison we plot the coincidence rate Γ as a function of both the mean and the variance of the stimulating current (Fig. 5G). The reproduction of spike trains on a spikebyspike basis is good (coincidence rate Γ > 0.7) over a broad range of stimulation parameters. Thus, even though the function F(u) has been optimized using data for fixed μ_{I} = 0 and σ_{I} = 25 μA/cm^{2}, the same NLIF model can also approximate spike trains of the full model when the stimulus has nonzero mean or higher variance. Only for very small fluctuation amplitudes σ_{I} < 20 μA/cm^{2} the approximation of the full model by a NLIF model breaks down. Raw data of Fig. 5G are reported in data supplements, Table 2, which also reports, in addition, results of a quadratic IF neuron where the 5 parameters (i.e., ϑ, τ_{m}, a_{0}, u_{c}, and R) have been optimized numerically (see Eq. 4). As expected, the results are very close to those of the NLIF model, but the performance is more sensitive to the exact choice of input variance than that of the NLIF model.
Comparison of the SRM and the full model
As indicated in the methods section, the 5 equations of the full conductancebased model of a fastspiking neuron model (Erisir et al. 1999) may also be approximated by a single differential equation (44) where the ion currents I^{ion} are exponential functions of the time t − t̂ that has elapsed since the last spike at t̂ (see Eqs. 28 and 29). We emphasize, that, in contrast to a standard leaky IF model the “time constant” τ is a function of t − t̂. Integration of Eq. 44 yields the equation of the SRM (45) The kernels η and κ, are given by exponential functions (see methods for details). If u(t) reaches a dynamic threshold ϑ(t − t̂) from below, then the neuron fires and t̂ is set to a new value t̂ = t, thus effectively resetting u(t), the time constant, and the ion currents. The parameters ϑ_{0}, ϑ_{1}, and τ_{θ} of the dynamic threshold (see methods, Eq. 33) are free parameters that we now adjust so as to optimize the performance for a given stimulation paradigm.
For pulse input of 2ms duration, the model neuron emits at most a single spike so that only the asymptotic threshold ϑ_{0} matters. For a sufficiently strong suprathreshold current pulse, spike initiation and hyperpolarizing spike afterpotential of the full model are reproduced by the SRM to a high degree of accuracy (cf. Fig. 6C). With ϑ = −50 mV the SRM produces action potentials whenever the amplitude of the current pulse is above 12.4 μA/cm^{2}. Delayed spikes caused by current pulses of 2 ms with amplitude I^{ext} between 8.8 and 12.4 μA/cm^{2} are not reproduced. If we adjusted the threshold ϑ_{0} to a lower value, the SRM would generate action potentials in this critical regime; the timing of the action potentials, however, would not be correct because the SRM—just as any other model with a strict voltage threshold—cannot account for delayed action potentials. Subthreshold excitation with amplitudes I^{ext} ≤ 7 μA/cm^{2} is reproduced to a fair degree of accuracy by the SRM.
We now turn to constant current input. For constant stimuli, a SRM with fixed threshold ϑ_{0} = −50 mV fails to approximate the gain function of the full model to a satisfactory degree (cf. Fig. 6B). We therefore use a dynamic threshold (Eq. 33; see methods) which approaches a value of ϑ_{0} = −50 mV with a time constant τ_{ϑ} = 6 ms and starting at a value of ϑ_{0} + ϑ_{1} = 0 mV at time t̂ + δ_{refr}. With the dynamic threshold, we get a fair approximation of the gain function of the full fastspiking neuron model. The approximation for currents that are just suprathreshold is bad, but for I^{ext} ≥ 5 μA/cm^{2} the rates are not too different from those of the full model. In Fig. 6A, the time course of the membrane potential of the SRM is compared to that of the full model. Even though the approach to threshold is not reproduced correctly, the period is about correct.
Similar to the case of the NLIF model, the choice of model parameters of the SRM is essential and depends on the regime the neuron model is optimized for. In the following we will focus on timedependent fluctuating current input and numerically optimize the SRM so as to reproduce a 10 s spike train of the full model, whereas both the full and reduced model are stimulated with random current of zero mean and variance σ_{I} = 25 μA/cm^{2}. For the numerical optimization of the SRM defined in Eq. 45 we need to determine the spike shape η, the inputresponse filter κ, as well as the parameters ϑ_{0}, ϑ_{1}, and τ_{ϑ} of the dynamic threshold (see methods). The kernel η is the average spike shape and (see Fig. 4, B and C). From Fig. 4D, we see that the kernel κ can be approximated by an exponential κ(δt, s) = exp[−s/τ(δt)] with a time constant τ(δt) shown in Fig. 4E. The threshold parameters are optimized so that most spikes occur with the correct timing (Fig. 4F). The optimal set of threshold parameters is ϑ_{0} = −53 mV, ϑ_{1} = 12.7 mV, and τ_{ϑ} = 54 ms (see also data supplements, Table 3). We keep these parameters fixed in the following. For the specific spike train used for parameter optimization, we find a coincidence factor Γ = 0.83.
We now turn to a systematic exploration of the performances of the numerically optimized SRM. A typical spike train is shown in Fig. 6D. To quantify the performance for different inputs, we first keep the variance of the input fixed and change the mean input μ_{I}. Figure 6E shows that the mean firing rates as a function of μ_{I} of the full model is comparable to that of the SRM. Moreover, the SRM reproduces the spike timing of the full model to a high degree of accuracy over a broad range of different σ_{I} and μ_{I}, as quantified by the coincidence rate Γ (cf. Fig. 6G). In data supplements, Fig. 1, it is furthermore shown that both the interval distribution and the coefficient of variation of the full model are correctly predicted by the simple model indicating that missed or added spikes do not modify significantly the spike pattern. Finally, in the subthreshold regime, the membrane potential of the reduced model approximates that of the full model within ±2 mV (cf. Fig. 6F). To summarize, the SRM predicts the subthreshold fluctuations, the correct number of spikes, most of them (>70%) with the correct timing and with a firing pattern very similar to the one produced by the full conductancebased neuron model.
The amplitude of the fluctuations of the stimulus (here, σ_{I}) is a crucial factor for the quality of the predictions within the SRM framework. This can be easily understood. First, if the variance of the input signal is large, the amplitude of fluctuations of the membrane voltage is large, which facilitates the emission of spikes at a correct timing with the threshold condition (Eq. 33) (Brette and Guignon 2003; Bryant and Segundo 1976; Mainen and Sejnowski 1995; Tanabe and Pakdaman 2001). Second, when the constant part of the driving current dominates the fluctuations (μ_{I} ≫ σ_{I}), whenever a spike is missed or added by the SRM, errors propagate further in time and the coincidence factor Γ decreases dramatically.
Comparison with simple threshold models
In the previous subsections, we have seen that both the NLIF and the SRM show a good performance for random input current, if the fluctuation amplitude is large enough. We may therefore wonder whether this is a universal property that would hold for any onedimensional threshold model. We therefore studied 3 models that are even simpler than the SRM or NLIF model. As a first simplification, we turn the dynamic threshold (Eq. 33) of the SRM into a constant one (except for the absolute refractory period) (46) We used the same spike train as before for optimization of the threshold. The optimal value for the parameter ϑ_{0} is −45.550 ± 0.001 mV (mean ± SD). Figure 7A shows the quality of predictions with a constant threshold. The results are significantly worse than those for the full SRM with dynamic threshold or those of the NLIF model (compare with Fig. 5G and 6G).
The SRM with constant threshold is very close to the standard LIF neuron. We simply use κ(t − t̂,s) = exp(−s/τ) for t > t̂ + δ_{refr} with τ = 4 ms and δ_{refr} = 2 ms, and κ[(t − t̂, s)] = 0 otherwise. This is equivalent to Eq. 2 of the LIF model. After firing, the membrane potential is reset to a value of −85 mV. The optimal parameter ϑ_{0} for the LIF neuron is −47.33 ± 0.08. Figure 7B shows the quality of predictions for the LIF neuron and, as expected, the results are close to the results of the SRM with constant threshold.
Let us now simplify even further. We use neither reset nor spike afterpotential η (i.e., we neglect the interaction between action potentials). Moreover, the kernel κ is replaced by lim_{t−t̂→∞} κ_{t−t̂,s}. The membrane voltage is then (47) We refer to this as minimal model 1. Specifically, κ_{∞,s} = κ_{0} exp(−s/τ) with τ ≈ 4 ms (i.e., a simple lowpass filter). As before, spikes are triggered whenever the variable u reaches a constant threshold ϑ from below. However, in contrast to the LIF model, there is no reset after firing. Although the minimal model 1 still performs well in the neighborhood of the optimization point, Γ decreases quickly as the stimulation parameters change (see data supplements, Table 2). It even goes below 0 for high frequencies and thus performs worse than a homogeneous Poisson neuron model (see definition of Γ in Eq. 22). This is attributed to the fact that for larger stimulation, the firing frequencies of minimal model 1 are too high and spikes occur systematically too early.
We now consider a minimal model 2, which is even simpler: instead of a lowpass filter we take κ as instantaneous (i.e., it takes a fixed positive value 1 during one time step and is zero thereafter). This is an approximation to a neuron in the highconductance state where the voltage follows the current quasi instantaneously. In this case, the variable u is given by (48) Firing is defined by the usual threshold condition. For both minimal models the only free parameter is the threshold ϑ_{0}, which has been optimized on a spike train generated by random input current with zero mean and variance σ_{I} = 25 μA/cm^{2}. Minimal model 2 performs poorly (see data supplements, Table 2). In particular, even at the point for which the parameter ϑ_{0} has been optimized, Γ attains a value of only 0.16 ± 0.01 (mean ± SD). The results with the minimal models indicate that the SRM and the NLIF model add features beyond simple threshold models and these features (i.e., spike–spike interaction for the SRM and nonlinear voltage dependency and voltage reset for the NLIF) are necessary if we want a model that works over a broad range of different inputs.
Extending the framework to slow processes
Real neurons often exhibit a rich repertoire of ion channels, in particular, slow ion channels that contribute to frequency adaptation, a widespread phenomenon in biological neurons. To take into account such slow processes, we need to leave the framework of singlevariable models and introduce an adaptation variable A. We use a simple relaxation dynamics (49) Each time the neuron emits a spike, the variable A is increased by an amount A_{1}/τ_{adapt}. In the absence of spikes A relaxes with time constant τ_{adapt} to an equilibrium value A_{0}. The sum on the righthand side of Eq. 49 runs over all output spikes of the neuron. For periodic firing with frequency v, the variable A fluctuates around a mean value of A_{0} + A_{1}v. Equations such as Eq. 49 yield standard phenomological models of adaptation (Benda and Herz 2003; Rauch et al. 2003).
To include adaptation into the SRM framework, we make the time constant τ_{∞} of the exponential response kernel κ (Fig. 4E) depend on A. In the case of the Erisir model of fastspiking interneurons, the variable n_{1} of slow K^{+} channels accumulates over consecutive spikes, so that the total conductance is increased and the effective membrane time constant τ_{∞} is shortened. We computed the response kernel κ(δt, s) (see methods) for different stimulation parameters and plotted the parameter 1/τ_{∞} as a function of the output frequency (Fig. 8A). A linear fit 1/τ_{∞}(v) = A_{0} + A_{1}v gives us the parameters A_{0} and A_{1} and the identification of the adaptation variable A in Eq. 49 as A = 1/τ_{∞}. The free parameter τ_{adapt} is set to a value of 450 ms. Figure 8, B–D shows the comparison between nonadapting SRM (Fig. 8C) and adapting SRM (Fig. 8D) for current injection with Gaussian white noise. After some time, parameters of the input are switched, leading to an increase in output frequency (Fig. 8B). The nonadapting SRM maintains τ_{∞} as tuned for the first part of the stimulation paradigm and thus performs poorly when stimulation parameters are changed (Fig. 8C, from left to right). The performance of the adapting SRM, however, are much better once adaptation has taken place (Fig. 8D). In terms of the coincidence factor Γ, the nonadapting SRM yields Γ = 0.63 when the threshold is optimized for intermediate stimulation (as in the rest of the article; Fig. 8C), whereas the adapting SRM yields Γ = 0.76 (Fig. 8D). The same methodology can also be applied to strongly adapting neuron models, as illustrated in Figure 8, E–H for the case of a 2compartment neuron model with Ca^{2+} channels and Ca^{2+}gated K^{+} channels responsible for slow adaptation (Wang 1998). Injection of random current in the soma allows us to identify the time constant τ_{∞}; a test with changing input confirms that the extended SRM can follow the adaptation of the full model.
Conductance injection scenario (stochastic spike arrival)
A more realistic input scenario than random current injection would be a model of stochastic conductance changes (Destexhe and Paré 1999). We preferred in the previous sections to work with the more common random current input because this allows us to characterize the input in a transparent fashion. Furthermore, for any random input scenario with stationary white noise characteristics, a random conductance scenario can be replaced, to a high degree of accuracy, by a random current scenario (Richardson 2004). Nevertheless, the numerical method exploited in previous sections can be extended to the case of random conductance injection. Instead of optimizing a kernel κ(t − t̂, s) that describes the linear response of the membrane potential to an input current, we now optimize two kernels ε^{exc}(t − t̂, t − t_{j}^{(f)}) and ε^{inh} (t − t̂, t − t_{j}^{(f)}) that describe the response of the membrane potential to spike arrival at an excitatory or inhibitory synapse. Here t_{j}^{(f)} denotes the time of arrival of spike number f at synapse number j. In other words ε^{exc} as a function of t − t_{j}^{(f)} describes the excitatory postsynaptic potential (EPSP), whereas ε^{inh} as a function of t − t_{j}^{(f)} describes the inhibitory postsynaptic potential (IPSP). In the following, we model stochastic spike arrival by Poisson point processes (see methods section). Samples of presynaptic spike trains and a histogram of spike arrival times are plotted in Fig. 9, A and B. For parameter optimization (see methods) we used 10 s of voltage data from the full model of a fastspiking neuron (Erisir et al. 1999) while stimulated by spike arrival at 8,000 excitatory model synapses (rate v_{exc} = 0.3 Hz each) and 2,000 inhibitory model synapses (rate ν_{inh} = 2.5 Hz each). Best fits for the dynamic threshold parameters are indicated in data supplements, Table 4.
As a result of the optimization procedure we get the shape of the kernels ε^{exc} and ε^{inh} shown in Fig. 9, D and E. Both EPSP and IPSP can be well approximated by double exponentials. For this specific spike train, the coincidence factor is Γ = 0.82.
We then keep the kernels ε^{exc} and ε^{inh} fixed and turn to a systematic exploration of the performances. Figure 10 shows the quality of predictions for different values of the discharge frequency of inhibitory synapses when holding the discharge frequency of excitatory synapses constant. The mean firing rates do not match as well as for current injection (Fig. 10A) except at the point at which parameters have been optimized. For this set of input parameters, the coincidence rate is good (Γ = 0.7; see Fig. 10B). The subthreshold behavior of the membrane voltage is also nicely predicted (Fig. 10, C and D). However, conductance injection is known to produce a suppression of subthreshold membrane voltage fluctuations (Monier et al. 2003) and a reduction of the membrane time constant (Destexhe and Paré 1999; Destexhe et al. 2001). We already mentioned that large membrane voltage fluctuations are a key point for correct prediction of the timing of the spikes. The fact that the effective membrane time constant depends on the stimulation paradigm implies that the EPSPs and IPSPs that we compute are optimal only in a limited range of input characteristics. This is an important difference to the current injection scenario where changes of the kernel κ with the input characteristics are less important. It is clear then that, in the case of conductance injection, our approach is valid only over a limited range of input characteristics with mean conductances close to those used to compute the εkernels, which explains the rapid decrease of Γ for low and high values of ν_{inh}. In particular, for reduced input rates (e.g., ν_{inh} = 1 Hz), the input conductance of the full model is reduced and hence, the membrane time constant is increased. Because the εkernels of the SRM do not change their time constant, the mean membrane potential of the SRM and thus the predicted mean firing rate are too low (see Fig. 10A).
DISCUSSION
In this paper, we mainly focused on deterministic singlevariable IF models (i.e., a spike is elicited if a single variable u crosses a threshold ϑ from below. If the evolution of u depends on the instantaneous voltage, we arrive at the NLIF model; if it depends on the time since the last spike, we arrive at the SRM. Models with 2 (Arcas et al. 2003; FitzHugh 1961; Nagumo et al. 1962) dimensions or more offer a wide spectrum of behavior but are often more difficult to optimize numerically. We proposed a straightforward extension of the SRM to the case of adapting neurons allowing a stepbystep procedure for numerical fitting. An extension toward neurons with hyperpolarizationactivated ion currents is conceivable but has not been explored. Another direction of research goes from deterministic to stochastic models of spike triggering (Arcas et al. 2003; Brillinger 1988; Brillinger and Segundo 1979; Gerstner 2001; Tuckwell 1988). This is important in regimes where the neuron shows a low degree of spike train reproducibility under repeated stimulation with the same signal (Bryant and Segundo 1976; Mainen and Sejnowski 1995).
Reduction of neuronal complexity to a singlevariable model implies that these neuron models are valid only in a limited regime. For example, the quadratic IF model as a canonical type I neuron model is optimal at very low firing frequencies; there is no a priori reason why it should perform well far from threshold, but it does so for some instances of models (Hansel and Mato 2003). Similarly, the SRM is based on a linearization about a reference state and there is no a priori reason why it should work well far from that reference point. Nevertheless, we found in this paper that the combination of linear summation with a dynamic threshold works for randomcurrent injection over a broad range of parameters. Because the SRM is based on the combination of linear summation with a sharp threshold, any input stimulus that probes specifically properties close to threshold will highlight the limits of the approach. In particular, we have seen that stimulation with a constant current just above threshold yields a frequency–current curve that is not captured by the SRM neuron, but well represented by a quadratic or nonlinear IF model. Moreover, stimulation by short current pulses with amplitudes close to the threshold amplitude will cause delayed responses of the Erisir model and other type I neuron model that cannot be captured by the sharp threshold process of the SRM, but which are perfectly accounted for by nonlinear IF models.
Given that all simplified models have a restricted regime of applicability, the big question is whether a model performs well in the biologically relevant regime. In our opinion, random input is the most realistic scenario and we focus our discussion therefore on this scenario. For random input, we measure the quality of simplified models by a coincidence factor Γ defined in Eq. 22. An alternative approach could be to derive a quality measure from information theory as proposed recently (Arcas et al. 2003). Although information theory provides a systematic theoretical framework, we think that our coincidence measure offers several practical advantages: in particular, it is easy to evaluate and allows a straightforward interpretation.
On our random input scenario, the MCIF model clearly yields the best performance. Although it is straightforward to implement and rapid to simulate, it is difficult to analyze mathematically. Strictly speaking, it does not fall in the class of singlevariable IF models. It is interesting to realize, however, that even the MCIF model, which is based on a seemingly innocent approximation, has a coincidence rate Γ significantly below one. The performance of the MCIF model could be improved by replacing the fixed reset values m_{reset}, h_{reset}, n_{1,reset}, n_{2,reset}, and u_{reset} by values that depend on the gating variables m, h, n_{1}, and n_{2} at time t̂ (i.e., at the moment of spike firing). We decided not to implement such a scheme because the main focus of this study has been on singlevariable IF models.
The onedimensional NLIF model that we used in this paper is very similar to the exponential IF model proposed recently (FourcaudTrocmé et al. 2003). A direct comparison of the NLIF model with the SRM shows that both yield a similar performance. The timedependent threshold that has been included in our definition of the SRM is an important component to achieve generalization over a broad range of different inputs. In fact, turning the dynamic threshold into a constant one reduces the stimulation regime where good predictions (Γ ≥ 0.7) are achieved. Further simplifications such as neglecting the reset have a dramatic effect when trying to generalize the predictions of the SRM over a broad range of different inputs.
Singlevariable IF models such as the NLIF or SRM models are easy and efficient in simulations. For numerical implementation, the SRM is simply reformulated as the equivalent IF model with timedependent time constant. The 2 models, SRM and NLIF, run at equal speed and with about the same performance. In summary, generalized IF models correctly predict up to 80% spike times of a much more complicated detailed neuron model driven by in vivo–like timedependent input. Our numerical methods for optimizing the parameters of generalized IF models can be directly applied to experimental data of real neurons. The only requirement is that a few seconds of intracellular recordings of membrane voltage during stimulation with a known timedependent input current are available.
APPENDIX
Optimal filtering problem
In this appendix, we adapt the classic Wiener–Hopf optimal filtering procedure (Wiener 1958) to our problem. We start with the usual crosscorrelation method to find the firstorder Wiener kernel (Lee and Schetzen 1965). At that point, a simple modification of the error function definition provides us with an equation that is simple to use for the SRM formalism. Note that the classic Wiener–Hopf method was already used in a similar context by Brillinger and Segundo (Brillinger 1988; Brillinger and Segundo 1979).
General case.
The classic formulation of the Wiener–Hopf optimal filtering problem is the following. Let us suppose that we record the output signal V_{t}^{data} from some device driven with an input signal I_{t}. We want to find the best linear filter F under the assumption that F is a finite impulse response filter so that (A1) is as close as possible from V_{t}^{data}. M_{−} and M_{+} are the largest negative and positive lags for which F is different from 0. We define the error for each bin by (A2) and the total error on some sample of data of length L by (A3) Then, the optimal linear filter F for the considered sample of data is the solution of the following system of equations (A4) Some algebra yields the Wiener–Hopf equation (A5) where is the empirical crosscorrelation between two vectors f and g at lag i. It is then very easy to find the optimal filter F by solving this linear system. For an application to our case, we need to take into account the fact that the optimal filter F is timedependent in the SRM formalism (the socalled kernel κ). More precisely, it depends on the time elapsed since the last emitted spike. We therefore replace F_{k} with F_{t−t̂,k} where t̂ is the last spike emitted when considering time t. Equation A1 becomes (A6) Again, we define the error for each bin by (A7) However, instead of defining the mean square error over a sample of data as in Eq. A3, we now sum over the set of firing times t_{i}, i = 1, … , T (A8) where s is simply t − t̂. Then, the optimal time dependent linear filter F for the whole spike train is the solution of the following system of equations (A9) After some algebra, we find a set of equations that are the equivalent of the Wiener–Hopf equation (Eq. A5) (A10) The quantity X[·, ·] is defined (for two vectors f and g) by (A11) Note that this latter quantity has 2 indices and is a generalization of the empirical crosscorrelation between 2 vectors. If the input I is a stationary random process, X[I, I]_{k} ∝ Corr [I, I]_{k} and strictly equivalent if L = T.
Two inputs with different response filters.
We suppose now that the device that we record from is driven by 2 different inputs and respond to each of these inputs with a different filter (A12) For the sake of simplicity, we assume that M_{−} and M_{+} are the same for both filters F^{1} and F^{2} and that both filters are causal (M_{−} = 0). Again, we want to find F^{1} and F^{2} so that V_{t}^{model} is as close as possible from V_{t}^{data}. Let us define F and I_{t} as (A13) (A14)
Equation A12 can now be rewritten as a vector product (A15) We define the error for each bin by (A16) and the total error on some sample of data of length L by (A17)
The optimal linear filters F^{1} and F^{2} for the considered sample of data are then the solution of the following system of equations (A18) Some algebra yields an equation that is equivalent to the classic Wiener–Hopf equation (Eq. A5) (A19) In its matrix form, this equation gives (A20) where T(f) is the hermitian Toeplitz matrix of the vector f. Note that in this case, the linear system includes terms of the crosscorrelations between inputs I^{1} and I^{2}.
For an application to our case, we need again to take into account the fact that the optimal filters F^{1} and F^{2} are timedependent in the SRM formalism (the socalled kernels ε). More precisely, they depend on the time elapsed since the last emitted spike. We follow the same way as above and replace both F_{k}^{1} and F_{k}^{2} by respectively F_{t−t̂,k}^{1} and F_{t−t̂,k}^{2}, where t̂ is the time of the last emitted spike when considering time t. Equation A12 is replaced by (A21) Using the same stratagem as before (with M_{−} = 0), we define (A22) (A23)
Thus V_{t}^{model} can be restated as a scalar product of the 2 vectors F_{t−t̂} and I_{t} (A24) The error for each bin is defined in Eq. A16. Just as in the case of a single input, we sum over the set of firing times t̂_{i}, i = 1,… , T for the total error function (A25) where s is short for t − t̂. The optimal filters F_{s}^{1} and F_{s}^{2} are then the solution of the following system of equations (A26) Some algebra yields a set of equations that are the equivalent of the Wiener–Hopf equation for this specific scenario (A27)
In its matrix form, this equation gives (A28)
See Eq. A11 for a definition of X[·, ·].
GRANTS
This work was supported by Swiss National Science Foundation (SNF) Grant NFN2100065268.
Acknowledgments
We thank M. Richardson for many helpful comments and discussions.
Footnotes
↵* All authors contributed equally to this work.
↵2 The Supplementary Material for this article (a figure and 4 tables) is available online at http://jn.physiology.org/cgi/content/full/00190.2004/DC1.

Address for reprint requests and other correspondence: R. Jolivet, Laboratory of Computational Neuroscience, EPFL, 1015 Lausanne, Switzerland (Email: renaud.jolivet{at}epfl.ch).
 Copyright © 2004 by the American Physiological Society