J Neurophysiol 92: 959-976, 2004;
doi:10.1152/jn.00190.2004
0022-3077/04 $5.00
Generalized Integrate-and-Fire Models of Neuronal Activity Approximate Spike Trains of a Detailed Model to a High Degree of Accuracy
Renaud Jolivet 1,*,
Timothy J. Lewis2,* and
Wulfram Gerstner1,*
1Laboratory of Computational Neuroscience, Swiss Federal Institute of Technology, École Polytechnique Fédérale de Lausanne 1015 Lausanne, Switzerland; and 2Center for Neural Science and Courant Institute of Mathematical Sciences, New York University, New York, New York 10003
Submitted 27 February 2004;
accepted in final form 18 March 2004
 |
ABSTRACT
|
|---|
We demonstrate that single-variable integrate-and-fire models can quantitatively capture the dynamics of a physiologically detailed model for fast-spiking cortical neurons. Through a systematic set of approximations, we reduce the conductance-based model to 2 variants of integrate-and-fire models. In the first variant (nonlinear integrate-and-fire 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 conductance-based model. To quantitatively test the predictive power of the SRM and of the nonlinear integrate-and-fire model, we compare spike trains in the simple models to those in the full conductance-based 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 integrate-and-fire model based on spike trains in the full conductance-based model. This technique can be used to tune simple models to reproduce spike trains of real neurons.
|
INTRODUCTION
|
|---|
Detailed conductance-based neuron models, sometimes termed Hodgkin–Huxley-type 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 integrate-and-fire (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–Huxley-type models by an essentially one-dimensional fire-and-reset 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 (Fourcaud-Trocmé et al. 2003) or the Spike Response Model (SRM) (Gerstner and Kistler 2002; Kistler et al. 1997) perform for a realistic, time-dependent 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 time-dependent 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 conductance-based 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 2-fold.
First, we attempt to illustrate the relation of phenomenological spiking neuron models to conductance-based 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 step-by-step analytical derivation of 2 formal spiking neuron models starting from a detailed conductance-based model of a fast-spiking cortical interneuron. In both cases, we compare the behavior of the reduced model to that of the detailed conductance-based 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 conductance-based model of an interneuron. The optimal parameters characterizing the reduced models are extracted from a sample spike train generated with the conductance-based 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 conductance-based 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
|
|---|
Integrate-and-fire 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 
is called the firing time of the neuron. We focus in this article on models that are fully described by a single variable u. Well-known 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 integrate-and-fire 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 Iext with the external driving current, we may interpret 
m as the membrane time constant, R as the input resistance, and ueq 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 ureset. 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 = , we interrupt the dynamics (Eq. 2) during an absolute refractory time 
refr and restart the integration at time + refr with the new initial condition ureset. 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 is recorded. In a general nonlinear integrate-and-fire 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 = ureset. A specific instance of a NLIF model is the quadratic model (Feng 2001; Hansel and Mato 2001; Latham et al. 2000)
 | (4) |
with parameters a0 > 0 and ueq < uc < . For Iext = 0 and initial conditions u < uc, the voltage decays to the resting potential ueq. For u > uc, the voltage increases up to where an action potential is triggered; uc 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 time-dependent input Iext(t). Let us denote the last firing time by . For 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; ureset is the initial condition of the integration at time + 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 ureset at time 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 – (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 time-dependent time constant (i.e., a SRM; Stevens and Zador 1998; Wehmeier et al. 1989)
 | (7) |
with C = m/R. Starting the integration at u( + refr) = ureset, we find for t > + refr
 | (8) |
which is a special case of Eq. 6. A time-dependent 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 > + refr.
If input is provided by the activation of a conductance-based 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 wj the strength of synapse j. For the sake of simplicity, we will assume throughout this article that all synapses have the same strength wj = 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 – tj(f)). Equation 10 can thus be restated
 | (11) |
This latter formulation will be used for the numerical fitting of the SRM to the full conductance-based model in case of random spike arrival.
Full conductance-based model
As our reference model, we take the conductance-based model of a fast-spiking cortical interneuron proposed by Erisir and coworkers (1999). We have chosen this specific model so that, because fast-spiking 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 fast-spiking neuron is comparatively simple, we can hope to illustrate the steps necessary for a reduction to spiking models in a transparent fashion.
The fast-spiking 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 INa is described by the equation
 | (13) |
with two gating variables m and h. A slow potassium current is modeled by
 | (14) |
with gating variable n1. Because of its slow time constant, the variable n1 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 gK2 that is much larger than gK1. 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/gI. Each gating variable follows the equation
 | (17) |
where x stands for m,h,n1, or n2. 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 fast-spiking neocortical interneurons. We adopt their set of parameters, except for the leak current where we have chosen gl = 0.25 mS/cm2 (instead of their value of 1.25 mS/cm2) 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 first-order method is sufficient.
Current injection
In the current injection scenario, the neuron is driven with an uncorrelated Gaussian-distributed random current. We vary both the mean µ and the standard deviation 
of the Gaussian current. This kind of highly variable time-dependent 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 Isyn is given by (Robinson and Kawai 1993)
 | (18) |
where gexc (ginh) is the total excitatory (inhibitory) conductance and Eexc (Einh) is the corresponding reversal potential. The synaptic conductances gexc and ginh consist of the summed input from Nexc excitatory synapses and Ninh 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 Aexc for each presynaptic spike arriving at the synapse at time tj,k, and decays with a time constant exc. The result is a low-pass–filtered version of the presynaptic spike train {tj,k}. The dynamics of inhibitory synapses are defined similarly.
Presynaptic spike trains are described by random homogeneous Poisson processes. At each time step, 
independent random variables are generated and distributed among the 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 N1 is the number of spikes in the reference spike train S1, N2 is the number of spikes in the spike train S2 that is compared to the reference spike train, Ncoinc is the number of coincidences with precision 
between the 2 spike trains, and 
Ncoinc
= 2
N1 is the expected number of coincidences generated by a homogeneous Poisson process with the same rate as the spike train S2. In this article, the reference spike train S1 is always generated by the full conductance-based model of a fast-spiking interneuron, whereas S2 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 S1 generated by the full conductance-based model and the train S2 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 fast-spiking 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 fast-spiking 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 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 ureset = –85 mV.
The first important approximation involves the reset of the gating variables m, h, n1, and n2, 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 hreset 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 Iext = 5 µA/cm2 that leads to a mean firing rate of about 40 Hz. The reset values are mreset = 0.0; hreset = 0.16; n1,reset = 0.874; n2,reset = 0.2; and ureset = –85 mV. This set of parameters leads to a near-perfect fit of the time course of the membrane potential during periodic firing at 40 Hz and approximates the gain function of the full fast-spiking neuron model to a high degree of accuracy (Fig. 1, A and B).
View larger version (22K):
[in this window]
[in a new window]
|
FIG. 1. Multicurrent integrate-and-fire (MCIF) model compared to the full conductance-based model. A: periodic firing. Neuron is driven by a constant input of Iext = 5 µA/cm2, which leads to regular firing at about 40 Hz. Trajectory of the membrane potential of the MCIF model (dashed line) is compared to that of the full model. (solid line). Action potentials in the MCIF model are replaced by triangular pulses starting at the threshold value of –45 mV, peaking at +65 mV, and ending at the reset value of –85 mV. Spike duration is 1.7 ms. B: gain function (frequency as a function of input current Iext) of the MCIF model (dashed line with circles) is compared to that of the full model (solid line). C: pulse input. A current pulse of 2-ms duration is delivered at t = 10 ms. MCIF model (dashed line) reproduces perfectly the firing times of the full model (solid line) to both strong (Iext = 20 µA/cm2, left) and weak (Iext = 8.8 µA/cm2, delayed action potential, right) superthreshold current pulses. Parameters of the reset have been adapted to a case where the neuron fires at about 40 Hz. Spike afterpotential is reproduced to a fair degree of accuracy, but not perfectly. D: for random input, the mean firing rate of the full model (solid line) is compared to that of the MCIF model for different amplitudes of the random current (dashed line with circles). E: spike train of the full model (solid line) is compared to that of the MCIF model (dashed line), whereas both neuron models receive exactly the same realization of a random current. For this input scenario, about 95% of the spike times are correct within ±2 ms. F: coincidence rate as a function of the amplitude of the random current. In the regime 25 < I < 35 µA/cm2, is about 0.95.
|
|
For stimulation with isolated 2-ms current pulses, the MCIF model reproduces both the threshold behavior and the firing times of the full fast-spiking neuron model. The hyperpolarizing spike-afterpotential, 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/cm2), 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, n1, and n2. To reduce it further to a single-variable 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 m0[u(t)] (see Fig. 2A). The treatment of the other gating variables deserves some extra discussion.
View larger version (23K):
[in this window]
[in a new window]
|
FIG. 2. Reduction from MCIF model to nonlinear IF (NLIF) model. All 4 variables, m, n2, h, and n1, are plotted in the periodic-firing regime when the neuron is stimulated by a constant current of 5 µA/cm2, which is switched on at t = 0. A: variable m (solid line) is compared to m0(u) (dashed line). Two trajectories coincide nearly perfectly. B: outside the refractory period, the variable n2 (solid line) is approximated by a constant value n2,eq. Vertical dotted lines delimit the 4-ms refractory period during and after spikes, during which the dynamics is not modeled. C: variable h (solid line) is compared to hav = 0.45 (dashed line). D: variable n1 (solid line) is compared to n1,av = 0.8 (dashed line).
|
|
We start with the variable n2. A thorough inspection of the time course of n2(t) shows that n2 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 n2 falls within the refractory period. Between spikes we can therefore replace n2 by its equilibrium value at rest n2 
n2,eq = n0(ueq).
The gating variables h and n1 vary only slowly so that, for a given input scenario, the variables may be replaced by their average value hav and n1,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., hav = 0.87 and n1,av = 0.00057). If the neuron is stimulated by a constant bias current that shifts the membrane potential just below threshold, we should take hav = 0.54 and n1,av = 0.019. Finally, in a periodic firing regime (constant stimulation with I = 5 µA/cm2) reasonable averages are hav = 0.45 and n1,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 = m0(u) and constant values for h, n1, and n2, the dynamics of the full fast-spiking 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 hav and n1,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 ueq and the critical voltage uc 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 = ueq. 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 ureset = –85 mV.
Analytical reduction to a spike response model
For a fitting of the full fast-spiking neuron model defined in Eqs. 12–17 to the SRM defined in Eqs. 6 and 9, we need to determine the kernels (t – ) and (t – , s), and furthermore the (time-dependent) threshold (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 as the time when the membrane potential crosses an (arbitrarily set) threshold (e.g., = –50 mV). The kernel (t – ) is defined by the time course of the membrane potential for 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 – ) = u() for 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, n1, and n2 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 + refr, similar to the approach of Destexhe (1997), we fit the approach to equilibrium by an exponential
 | (25) |
where x = m, h, n1, n2 stands for each of the 4 gating variables. The parameter x is a fixed time constant and xreset is the initial condition at t = + refr, and xeq = x (ueq). In other words, the voltage-dependent differential Eq. 17 for x is replaced by the linear voltage-independent equation
 | (26) |
Numerical values for xeq, x, and xreset 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 IK2 is
 | (27) |
where g2(t – ) is an exponential function with time constant n2/2. We insert the time-dependent conductance into Eq. 12 and find for t + refr
 | (28) |
where the sum runs over the 4 ion channels INa, IK1, IK2, and I1. With the definitions (t – ) = C/j gj(t – ) and Iion(t – ) = jgj(t – ) Ej, we arrive at
 | (29) |
which is a linear differential equation with time-dependent time constant (cf. Eq. 7). The effective time constant is shown in Fig. 3C as a function of t – .
According to Eqs. 27 and 28, the ion currents have, at time + refr, a value of Iion(refr) = j gj(refr)Ej. For 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 Iion(refr). For t + refr, we therefore integrate Eq. 29 with the initial condition
 | (30) |
This yields the SRM in the form of Eq. 6 for t > – refr with the kernels given by
 | (31) |
 | (32) |
The kernels and that we have constructed so far are limited to t > + refr. During the absolute refractory period < 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 < + 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., I0 = 5 µA/cm2 and I0 = 30 µA/cm2) 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 conductance-based 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 utdata and the discretized driving current Itext of a cell. We assume that both utdata and Itext 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 left-hand side, we subtract the effects of the action potentials and arrive at an expression for the subthreshold regime, which we denote as Vt–,tSRM
 | (37) |
Similarly, for the data utdata, we subtract the time course of the average action potential that yields a subthreshold behavior Vt–,tdata
 | (38) |
The right-hand side of Eq. 37 can be interpreted as a numerical convolution between It and a family of filters parameterized by t – . The problem of finding the best linear filter (t–,s) of some device given the output vector (Vt–,sdata) and the input vector (It) 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 – . 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 – . We find that t–,s is a solution of the following linear system
 | (39) |
where smax is the maximum time lag for which t,s is nonzero, 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 cross-correlation 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 time-dependent 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 Gaussian-distributed random current (cf. Fig. 4F). Parameters of the current are µI = 0 µA/cm2 and 1 = 25 µA/cm2. The value of the driving current is changed every 0.2 ms. Both the membrane voltage and the driving current are sampled at 5-kHz 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 Stexc and Stinh are the histograms of spike arrivals at excitatory and inhibitory synapses, respectively. More precisely, a value Stexc = 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, Sexc (Sinh) 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 conductance-based model of a fast-spiking 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 integrate-and-fire (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 th
TOP
ABSTRACT
INTRODUCTION
4 ms (dashed line).
)]} with 