JN Email Content Delivery
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


J Neurophysiol 92: 959-976, 2004; doi:10.1152/jn.00190.2004
0022-3077/04 $5.00
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Supplemental Data
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Web of Science (41)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Jolivet, R.
Right arrow Articles by Gerstner, W.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Jolivet, R.
Right arrow Articles by Gerstner, W.

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
 
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
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
 
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
Detailed conductance-based neuron models, sometimes termed Hodgkin–Huxley-type models (Hodgkin and Huxley 1952Go), can reproduce electrophysiological measurements to a high degree of accuracy (for a review, see Bower and Beeman 1995Go; Koch and Segev 1989Go). 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 1966Go; Hill 1936Go; Stein 1965Go; Tuckwell 1988Go) are highly popular and have been used to discuss aspects of neural coding, memory, or network dynamics (for a review, see Gerstner and Kistler 2002Go; Maass and Bishop 1998Go).

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 1948Go)] in the limit of very low firing rates yields a canonical type I model (Ermentrout 1996Go; Ermentrout and Kopell 1986Go; Hoppensteadt and Izhikevich 1997Go) that is equivalent to a quadratic IF model (Hansel and Mato 2001Go; Latham et al. 2000Go). 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. 2003Go) or the Spike Response Model (SRM) (Gerstner and Kistler 2002Go; Kistler et al. 1997Go) 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)Go 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 2002Go) 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. 2003Go) 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. 2001Go). Quadratic IF neurons can approximate the frequency–current curve of a detailed conductance-based model (Hansel and Mato 2003Go). 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 2001Go; Kistler et al. 1997Go). 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 1990Go; Destexhe 1997Go; Ermentrout 1996Go; Ermentrout and Kopell 1986Go; Kistler et al. 1997Go; Latham et al. 2000Go). 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 1968Go; Destexhe and Paré 1999Go).

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 1988Go; Brillinger and Segundo 1979Go; Wiener 1958Go). 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é 1999Go; Destexhe et al. 2001Go; Harsch and Robinson 2000Go; Reyes 2003Go).


 METHODS
 
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
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 {vartheta} 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 {tau}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 {delta}refr and restart the integration at time + {delta}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 1993Go)

(3)
As before, the dynamics is stopped if u reaches the threshold {vartheta} and reinitialized at u = ureset. A specific instance of a NLIF model is the quadratic model (Feng 2001Go; Hansel and Mato 2001Go; Latham et al. 2000Go)

(4)
with parameters a0 > 0 and ueq < uc < {vartheta}. For Iext = 0 and initial conditions u < uc, the voltage decays to the resting potential ueq. For u > uc, the voltage increases up to {vartheta} 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 {Theta}-neuron, a canonical type I neuron model (Ermentrout 1996Go; Hoppensteadt and Izhikevich 1997Go).

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 > + {delta}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 + {delta}refr. If we replace the exponentials multiplied by a step function by "response kernels" {eta} and {kappa}, 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 1995Go; Kistler et al. 1997Go). The kernel {eta} models, in the general case, the spike itself and the afterhyperpolarization that follows the spike. In the specific case of Eq. 5, {eta} describes the hard reset of u(t) to ureset at time t = + {delta}refr. The kernel {kappa} 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 {tau}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 1998Go; Wehmeier et al. 1989Go)

(7)
with C = {tau}m/R. Starting the integration at u( + {delta}refr) = ureset, we find for t > + {delta}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 1962Go; Geisler and Goldberg 1966Go; Holden 1976Go; Stein 1967Go; Weiss 1966Go)

(9)
During an absolute refractory period {delta}refr, we may, for example, set {vartheta} to a large and positive value to avoid firing and let it relax back to its equilibrium value for t > + {delta}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 {kappa} is replaced by a sum and a response kernel {varepsilon} 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 {varepsilon} 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) = {Sigma}j {Sigma}f {delta}(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)Go. 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. 1999Go) 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. 1999Go). Finally, there is the linear leak current

(16)
which defines a passive membrane time constant {tau}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)Go 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. 2002Go). 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.


View this table:
[in this window]
[in a new window]
 
TABLE 1. Fast-spiking neuron model

 
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 {sigma} 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 1968Go; Destexhe and Paré 1999Go).

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 1993Go)

(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 2001Go)

(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 {tau}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é 1999Go). Numerical values are summarized in Table 2.


View this table:
[in this window]
[in a new window]
 
TABLE 2. Parameters of excitatory and inhibitory synapses

 
Coincidence factor {Gamma}

The coincidence factor {Gamma} between two spike trains (Kistler et al. 1997Go) 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 {Delta} between the 2 spike trains, and <Ncoinc> = 2{nu}{Delta}N1 is the expected number of coincidences generated by a homogeneous Poisson process with the same rate {nu} 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{nu}{Delta} normalizes {Gamma} 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 {Gamma} = 0. We compute the coincidence factor {Gamma} 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, {Gamma} 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 1990Go; Kepler et al. 1992Go; Rinzel 1985Go). 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 {delta}refr.

The MCIF model is defined by a voltage threshold {vartheta}, a refractory time {delta}refr, and the reset values from which the integration is restarted. We fix the threshold at {vartheta} = –40 mV; the exact value is not critical, and we could take values of –20 to –45 mV without significantly changing the results. With {vartheta} = –40 mV, a refractory time of {delta}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 {nu} 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 {sigma} 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 {Gamma} as a function of the amplitude of the random current. In the regime 25 < {sigma}I < 35 µA/cm2, {Gamma} 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 {sigma}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 {Delta} = ±2 ms) by the MCIF model. We see from Fig. 1F that the coincidence rate {Gamma} 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 {Gamma} 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{infty}(u), or as slow compared to u, in which case we take x as constant (see also Abbott and Kepler 1990Go; Kepler et al. 1992Go; Rinzel 1985Go). In our case, m(t) is the only fast variable (in the subthreshold range [–100, –40] mV, we find 0.08 ≤ {tau}m ≤ 0.25 ms, 4.28 ≤ {tau}h ≤ 14.45 ms, 44.66 ≤ {tau}n1 ≤ 144.10 ms and 0.44 ≤ {tau}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 {delta}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 {equiv} 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. 1217 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: {tau} = –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/{tau} and then combined with a threshold at, say, {vartheta} = –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 1996Go; Ermentrout and Kopell 1986Go; Hansel and Mato 2001Go; Latham et al. 2000Go).

In the following, we do not make any further approximations. Instead we work directly with the nonlinear function F(u). The refractory period {delta}refr = 4 has already been introduced above. To finish the definition of the model, we have to specify the threshold {vartheta} and the reset value. We take {vartheta} = –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. 1217 to the SRM defined in Eqs. 6 and 9, we need to determine the kernels {eta}(t) and {kappa}(t, s), and furthermore the (time-dependent) threshold {vartheta}(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 {vartheta} (e.g., {vartheta} = –50 mV). The kernel {eta}(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 {eta}(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 {delta}refr = 2 ms. For t ≥ + {delta}refr, similar to the approach of Destexhe (1997)Go, 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 {tau}x is a fixed time constant and xreset is the initial condition at t = + {delta}refr, and xeq = x{infty} (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, {tau}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.



View larger version (10K):
[in this window]
[in a new window]
 
FIG. 3. Dynamics during and after an isolated action potential that has been elicited at time t = 0. A: time course of the variables m (solid line) and h (dotted line) in the full model are, after an absolute refractory period of 2 ms, approximated by exponentials with time constants {tau}m = 8.0 ms (long-dashed line; note that this is a different time constant than the membrane time constant) and {tau}h = 8.0 ms (dash-dotted line). B: similarly, the dynamics of n1 (solid line) and n2 (dotted line) are approximated by exponentials with time constants {tau}n2 = 0.75 ms (dash-dotted line) and {tau}n1 = 100 ms (long-dashed line). In A and B, stars denote the starting point of the exponential approximation after the 2-ms refractory period. C: Spike Response Model with kernels 31 and 32 can be interpreted as an IF model with a time constant {tau} that depends on the time {delta}t since the last spike. Integration restarts after an absolute refractory period of 2 ms with a time constant of 0.1 ms. The time constant (solid line) relaxes first rapidly and then more slowly toward its equilibrium value of {tau} {approx} 4 ms (dashed line).

 

View this table:
[in this window]
[in a new window]
 
TABLE 3. Parameter values for the semianalytically derived SRM

 
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 {tau}n2/2. We insert the time-dependent conductance into Eq. 12 and find for t ≥ + {delta}refr

(28)
where the sum runs over the 4 ion channels INa, IK1, IK2, and I1. With the definitions {tau}(t) = C/{Sigma}j gj(t ) and Iion(t) = {Sigma}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 + {delta}refr, a value of Iion({delta}refr) = {Sigma}j gj({delta}refr)Ej. For t = {delta}refr, the effective time constant {tau}({delta}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({delta}refr). For t ≥ + {delta}refr, we therefore integrate Eq. 29 with the initial condition

(30)
This yields the SRM in the form of Eq. 6 for t > {delta}refr with the kernels given by

(31)

(32)

The kernels {kappa} and {eta} that we have constructed so far are limited to t > + {delta}refr. During the absolute refractory period < t < + {delta}refr the neuron is not responsive to external input. We therefore set {kappa} to zero. The action potential itself that occurs in the interval < t < + {delta}refr is, for the sake of simplicity, approximated by a triangular voltage pulse. Finally, we introduce a dynamic threshold

(33)
with constants {vartheta}0, {vartheta}1, {vartheta}refr, and {tau}{vartheta}. During the absolute refractory period {delta}refr, the threshold is set to a value {vartheta}refr = 100 mV that is sufficiently high to prevent the neuron from firing. After refractoriness, the threshold starts at {vartheta}0 + {vartheta}1 and relaxes with a time constant of {tau}{vartheta} to an asymptotic value of {vartheta}0. The initial value {vartheta}0 + {vartheta}1 was arbitrarily set at 0 mV. {vartheta}0 = –50 mV and {tau}{vartheta} = 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 {eta}, describing the shape of a spike; the kernel {kappa}, describing the input integration process and finally the threshold {vartheta} 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 {eta} from the course of the membrane voltage. Second, we use a linear approximation for the integration of the input to extract {kappa} 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 {vartheta} 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 {eta} 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 {eta} (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).



View larger version (23K):
[in this window]
[in a new window]
 
FIG. 4. Numerical reduction of the full conductance-based model to a Spike Response Model (SRM). A: sample of the spike train used for the reduction (top line) and the corresponding time derivative (bottom line). Dashed line indicates the threshold used to detect and align spikes (80 mV/ms). Horizontal bar is 50 ms. Vertical bars are 100 mV (top) and 300 mV/ms (bottom). B: about 300 spikes aligned with the criterion discussed in A. C: kernel {eta} resulting from averaging the membrane voltage time course of the action potentials shown in B. D: kernel {kappa}({delta}t, s) as a function of time s extracted by the method described in the text for different time delays {delta}t = t: {delta}t = 0 ms (squares), {delta}t = 6 ms (crosses), and {delta}t = 10 ms (circles). Solid lines are exponential fits with {kappa}({delta}t, s) = {kappa}0 exp[–s/{tau}({delta}t)]. E: time constant {tau}({delta}t) of the kernel {kappa}. Diamonds are points resulting of fit of kernels as in D. Time constant {tau}({delta}t) of the kernel {kappa} can be approximated by the function {tau}({delta}t) = ({tau}{infty}/2){1 + tanh [{alpha}({delta}t{beta})]} with {tau}{infty} = 3.72 ms, {alpha} = 0.17, and {beta} = 6.64 ms (dashed line). F: segment of the spike train of the full model used for reduction (solid line) is compared to that of the reduced SRM (long-dashed line). Subthreshold behavior of the membrane voltage is well reproduced. Timing of the first spike is off by 1.6 ms and that of the second spike by 0.6 ms. Both spikes are therefore counted as coincident.

 
2) After the extraction of {eta}, we move to the extraction of {kappa}. 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 {sigma}1 are constants and {xi}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 {eta}-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 {sigma}I{xi}t and a family of filters parameterized by t. The problem of finding the best linear filter ({kappa}t–,s) of some device given the output vector (Vt,sdata) and the input vector ({sigma}I{xi}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. 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 {kappa}t–,s is a solution of the following linear system

(39)
where smax is the maximum time lag for which {kappa}{delta}t,s is nonzero, {delta}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 {kappa}{delta}t,s for each {delta}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 {tau} of this function versus the delay {delta}t and fit the dependency on {delta}t by the function {tau}({delta}t) = {alpha}1(1 + tanh ({alpha}2({delta}t{alpha}3)]) with free parameters {alpha}1, {alpha}2, and {alpha}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 {delta}refr is always kept at 2 ms and {vartheta}refr at a value of 100 mV to prevent firing. {vartheta}0, {vartheta}1, and {tau}{theta} are free parameters. To find the best fit for these 3 parameters, we simulate the full spike train with the SRM and compute 1 – {Gamma}. For a definition of {Gamma}, see Eq. 22. Then, we use the downhill simplex method (Nelder and Mead 1965Go) algorithm to find a set of parameters that minimize 1 – {Gamma}.

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 {sigma}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 {varepsilon}-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 {varepsilon}exc and {varepsilon}inh both at the same time (see APPENDIX). The filters {varepsilon}exc and {varepsilon}inh are fitted by double exponentials (see Fig. 9). The kernel {eta} 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 {nu}exc = 0.3 Hz and {nu}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.



View larger version (37K):
[in this window]
[in a new window]
 
FIG. 9. Conductance injection. A: raster plot of the excitatory neurons (left scale) and histogram of spike arrival times Sexc(t) with 0.2 ms time resolution (right scale). B: same as in A for inhibitory neurons. Black box corresponds to the temporal segment plotted in C. C: membrane potential as a function of time. Response of the full conductance-based model (solid line) to the spike input shown in A and B is compared to the prediction of the SRM (dashed line). All the spikes are predicted correctly within ±2 ms by the SRM and the predicted membrane voltage is most of the time almost indistinguishable from that of of the full conductance-based model. For this specific scenario, the coincidence factor evaluated over 1 s of data is {Gamma} = 0.72. D: {varepsilon}exc-kernels for stochastic spike arrival scenario extracted with the numerical method detailed in the APPENDIX. {varepsilon}exc-kernels are plotted for various t delays (t = 4 ms, solid line; t = 6 ms, dashed line; and t = 60 ms, dotted line) plotted on top of the {eta}-kernel (solid line). Arrows indicate the timing of activation of a single excitatory synapse. Inset shows the same 3 kernels aligned relatively to the onset of activation. Immediately after action potential emission, the kernel {varepsilon}exc is larger in amplitude than that at rest because of increased driving force of the synaptic current during the hyperpolarizing spike afterpotential. Kernel {varepsilon}exc is also shorter because of a shorter effective membrane time constant just after emission of a spike (see Fig. 4C). These effects vanish roughly 5–10 ms after the spike initiation. E: same as in D except that inhibitory {varepsilon}-kernels are plotted. Note that amplitude of the kernels plotted on top of {eta} is 5 times larger than normal for legibility reasons. Horizontal thick solid line indicates the reversal potential of inhibitory synapses (Einh = –80 mV). Close to the spike, membrane voltage is usually below Einh and activation of an inhibitory synapse induces a depolarization of the neuron. This effect reverses around t = 5 ms and leads to biphasic kernels.

 

 RESULTS
 
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
We have compared a detailed conductance-based model of a fast-spiking neuron (Erisir et al. 1999Go), 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 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 2-ms 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 fast-spiking neuron (Erisir et al. 1999Go) by a NLIF with a single variable u, given by the differential equation

(43)
The NLIF neuron fires when u hits a threshold of {vartheta} = –45 mV. Integration then restarts after a refractory period of 4 ms at ureset = –85 mV. The function F(u) is shown in Fig. 5F and parameterized by 3 constants, hav, n1,av, and n2,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 short-dashed 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 long-dashed line. Thus although it is always possible to approximate the dynamics of the full model by a one-dimensional 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 fast-spiking neuron model for 3 different input scenarios.



View larger version (26K):
[in this window]
[in a new window]
 
FIG. 5. NLIF 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 NLIF model (dashed line) is compared to that of the full model (solid line). Action potentials in the NLIF 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 4 ms. B: periodic firing. Gain function (frequency vs. current) of the full fast-spiking neuron model (solid line) is compared to that of the NLIF model with parameters adjusted to either the constant subthreshold current regime (long-dashed line with circles) or to a periodic firing regime (long-dashed line with squares). C: pulse input. Response of the full model (solid line) is compared to that of the NLIF model (dashed line) for input Iext = 20 µA/cm2 of 2-ms duration applied at t = 10 ms. Action potential is replaced by a triangular pulse spanning the refractory period of {delta}refr = 4 ms. Note the difference in the resting potential between NLIF and the full model, attributed to the fact that we used the set of parameters optimal for firing at 40 Hz. D: spike train of the NLIF model (dashed line) is compared to that of the full fast-spiking neuron model (solid line; {sigma}I = 45 µA/cm2, µI = 0 µA/cm2). The two traces are most of the time indistinguishable. E: mean firing rate {nu} of the NLIF (solid line with squares) is compared to that of the full model (dashed line with circles). F: function F(u) of the NLIF neuron of Eq. 24 with parameters adapted to the resting state (solid line) is compared to that adapted to a state with subthreshold constant current application (short-dashed line) or periodic firing at 40 Hz (long-dashed line). Zero crossings F(u) = 0 define the resting potential ueq and the critical voltage uc for spike initiation by current pulses. Note that the function F adapted to periodic firing has a different resting potential than that of the other 2 functions. G: contour plot of the coincidence factor {Gamma} for the NLIF model when compared to the full fast-spiking neuron model as a function of both mean current drive µI and amplitude of the random fluctuations {sigma}I. Each of the data points (squares and triangles) is the mean computed over 5 spike trains of 10 s each. Surface is interpolated with cubic splines (MATLAB, MathWorks, Natick, MA) and the circle highlights the stimulation paradigm used for parameter optimization. Points indicated by squares are reported in DATA SUPPLEMENTS, Table 3.

 
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 2-ms 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 2-ms pulse is above 8.5 µA/cm2, which is close to the pulse initation threshold of 8.8 µA/cm2 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 Iext ≥ 15 µA/cm2, whereas the full model triggers spikes for Iext ≥ 8.8 µA/cm2. 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 1968Go; Destexhe and Paré 1999Go). 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 hav, n1,av, and n2,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 1965Go) using a fluctuating input current of zero mean and variance {sigma}I = 25 µA/cm2 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 time-dependent 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 spike-by-spike basis. In particular, we maximize the coincidence rate {Gamma} 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 {delta}refr = 2 ms; the performance does not change significantly if we take {delta}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 {sigma}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 [{sigma}I > 50 (µA/cm2)] of the input current. For a more detailed comparison we plot the coincidence rate {Gamma} as a function of both the mean and the variance of the stimulating current (Fig. 5G). The reproduction of spike trains on a spike-by-spike basis is good (coincidence rate {Gamma} > 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 {sigma}I = 25 µA/cm2, 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 {sigma}I < 20 µA/cm2 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., {vartheta}, {tau}m, a0, uc, 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 conductance-based model of a fast-spiking neuron model (Erisir et al. 1999Go) may also be approximated by a single differential equation

(44)
where the ion currents Iion are exponential functions of the time t that has elapsed since the last spike at (see Eqs. 28 and 29). We emphasize, that, in contrast to a standard leaky IF model the "time constant" {tau} is a function of t. Integration of Eq. 44 yields the equation of the SRM

(45)
The kernels {eta} and {kappa}, are given by exponential functions (see METHODS for details). If u(t) reaches a dynamic threshold {vartheta}(t) from below, then the neuron fires and is set to a new value = t, thus effectively resetting u(t), the time constant, and the ion currents. The parameters {vartheta}0, {vartheta}1, and {tau}{theta} 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 2-ms duration, the model neuron emits at most a single spike so that only the asymptotic threshold {vartheta}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 {vartheta} = –50 mV the SRM produces action potentials whenever the amplitude of the current pulse is above 12.4 µA/cm2. Delayed spikes caused by current pulses of 2 ms with amplitude Iext between 8.8 and 12.4 µA/cm2 are not reproduced. If we adjusted the threshold {vartheta}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 Iext ≤ 7 µA/cm2 is reproduced to a fair degree of accuracy by the SRM.



View larger version (26K):
[in this window]
[in a new window]
 
FIG. 6. Predictions of the SRM compared to that of the full conductance-based model. A: spike train of the full model (solid line) is compared with that of the SRM (dashed line) for constant input Iext = 10 µA/cm2. Even though the interspike interval is correctly reproduced, the trajectory of the membrane potential is significantly different. B: gain function (frequency vs. current) of the full model (solid line) and the SRM with constant (dotted line with triangles) and dynamic threshold (long-dashed line with triangles). C: pulse input. Response of the SRM (dashed lines) to short current pulses (2 ms duration) is compared to that of the full model (solid lines). Action potentials caused by strong pulses (left: I = 20 µA/cm2; middle: I = 12.5 µA/cm2) are reproduced nicely by the SRM, whereas the delayed pulse that the full model exhibits for weaker input (right: I = 8.8 µA/cm2) is missed. D: spike train of the SRM model (dashed line) is compared to that of the full fast-spiking neuron model (solid line; {sigma}I = 45 µA/cm2, µI = 0 µA/cm2). The two traces are most of the time indistinguishable. E: mean firing rate {nu} of the SRM (solid line with squares) as a function of the amplitude {sigma}I of current fluctuations compared to that of the full model (dashed line with circles). F: normalized histogram of the error in prediction of the subthreshold membrane voltage bin after bin for one specific spike train (light gray area) with a fitted Gaussian distribution (solid line; {chi}2 = 2.5 x 10–5; mean = 1.0 mV; SD = 2.2 mV).

 
We now turn to constant current input. For constant stimuli, a SRM with fixed threshold {vartheta}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 {vartheta}0 = –50 mV with a time constant {tau}{vartheta} = 6 ms and starting at a value of {vartheta}0 + {vartheta}1 = 0 mV at time + {delta}refr. With the dynamic threshold, we get a fair approximation of the gain function of the full fast-spiking neuron model. The approximation for currents that are just suprathreshold is bad, but for Iext ≥ 5 µA/cm2 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 time-dependent 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 {sigma}I = 25 µA/cm2. For the numerical optimization of the SRM defined in Eq. 45 we need to determine the spike shape {eta}, the input-response filter {kappa}, as well as the parameters {vartheta}0, {vartheta}1, and {tau}{vartheta} of the dynamic threshold (see METHODS). The kernel {eta} is the average spike shape and (see Fig. 4, B and C). From Fig. 4D, we see that the kernel {kappa} can be approximated by an exponential {kappa}({delta}t, s) = exp[–s/{tau}({delta}t)] with a time constant {tau}({delta}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 {vartheta}0 = –53 mV, {vartheta}1 = 12.7 mV, and {tau}{vartheta} = 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 {Gamma} = 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 {sigma}I and µI, as quantified by the coincidence rate {Gamma} (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 conductance-based neuron model.

The amplitude of the fluctuations of the stimulus (here, {sigma}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 2003Go; Bryant and Segundo 1976Go; Mainen and Sejnowski 1995Go; Tanabe and Pakdaman 2001Go). Second, when the constant part of the driving current dominates the fluctuations (µI >> {sigma}I), whenever a spike is missed or added by the SRM, errors propagate further in time and the coincidence factor {Gamma} 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 one-dimensional 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 {vartheta}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).



View larger version (43K):
[in this window]
[in a new window]
 
FIG. 7. A: coincidence factor {Gamma} for an SRM with constant threshold. B: coincidence factor {Gamma} for the numerically optimized LIF. For details, see Fig. 5G.

 
The SRM with constant threshold is very close to the standard LIF neuron. We simply use {kappa}(t,s) = exp(–s/{tau}) for t > + {delta}refr with {tau} = 4 ms and {delta}refr = 2 ms, and {kappa}[(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 {vartheta}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 {eta} (i.e., we neglect the interaction between action potentials). Moreover, the kernel {kappa} is replaced by limt–->{infty} {kappa}t–,s. The membrane voltage is then

(47)
We refer to this as minimal model 1. Specifically, {kappa}{infty},s = {kappa}0 exp(–s/{tau}) with {tau} {approx} 4 ms (i.e., a simple low-pass filter). As before, spikes are triggered whenever the variable u reaches a constant threshold {vartheta} 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, {Gamma} 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 {Gamma} 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 low-pass filter we take {kappa} 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 high-conductance 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 {vartheta}0, which has been optimized on a spike train generated by random input current with zero mean and variance {sigma}I = 25 µA/cm2. Minimal model 2 performs poorly (see DATA SUPPLEMENTS, Table 2). In particular, even at the point for which the parameter {vartheta}0 has been optimized, {Gamma} 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 single-variable 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 A1/{tau}adapt. In the absence of spikes A relaxes with time constant {tau}adapt to an equilibrium value A0. The sum on the right-hand 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 A0 + A1v. Equations such as Eq. 49 yield standard phenomological models of adaptation (Benda and Herz 2003Go; Rauch et al. 2003Go).

To include adaptation into the SRM framework, we make the time constant {tau}{infty} of the exponential response kernel {kappa} (Fig. 4E) depend on A. In the case of the Erisir model of fast-spiking interneurons, the variable n1 of slow K+ channels accumulates over consecutive spikes, so that the total conductance is increased and the effective membrane time constant {tau}{infty} is shortened. We computed the response kernel {kappa}({delta}t, s) (see METHODS) for different stimulation parameters and plotted the parameter 1/{tau}{infty} as a function of the output frequency (Fig. 8A). A linear fit 1/{tau}{infty}(v) = A0 + A1v gives us the parameters A0 and A1 and the identification of the adaptation variable A in Eq. 49 as A = 1/{tau}{infty}. The free parameter {tau}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 {tau}{infty} 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 {Gamma}, the nonadapting SRM yields {Gamma} = 0.63 when the threshold is optimized for intermediate stimulation (as in the rest of the article; Fig. 8C), whereas the adapting SRM yields {Gamma} = 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 2-compartment neuron model with Ca2+ channels and Ca2+-gated K+ channels responsible for slow adaptation (Wang 1998Go). Injection of random current in the soma allows us to identify the time constant {tau}{infty}; a test with changing input confirms that the extended SRM can follow the adaptation of the full model.



View larger version (44K):
[in this window]
[in a new window]
 
FIG. 8. Extending the framework to adaptation. A: 1/{tau}{infty} is plotted vs. the output frequency of the full Erisir model (squares) for different parameters (µI, {sigma}I) of random current stimulation. Dependency is fitted by a linear relation (solid line). The two circles correspond to the stimulus used in B. B: Erisir model is stimulated with random input current whose parameters are instantaneously switched after 5 s. Driving current (top) around the transition time and response of the full Erisir model (bottom, total duration of segment is 700 ms). C: spike train of the full model (solid line) is compared to that predicted by the SRM (dashed line) optimized for initial stimulation paradigm. From left to right: 20 ms samples corresponding to the gray boxes plotted in B (bottom), respectively, before, shortly after, and more than 300 ms after the transition. D: same as in C, except that {tau}{infty} adapts to the output frequency with an arbitrarily set time constant {tau}adapt = 450 ms. In C and D, the value of {tau}{infty} (ms) is written for each sample. E–H: same as A–D when SRM is compared to the full Wang model (Wang 1998Go).

 
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é 1999Go). 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 2004Go). Nevertheless, the numerical method exploited in previous sections can be extended to the case of random conductance injection. Instead of optimizing a kernel {kappa}(t, s) that describes the linear response of the membrane potential to an input current, we now optimize two kernels {varepsilon}exc(t, t tj(f)) and {varepsilon}inh (t, ttj(f)) that describe the response of the membrane potential to spike arrival at an excitatory or inhibitory synapse. Here tj(f) denotes the time of arrival of spike number f at synapse number j. In other words {varepsilon}exc as a function of ttj(f) describes the excitatory postsynaptic potential (EPSP), whereas {varepsilon}inh as a function of ttj(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 fast-spiking neuron (Erisir et al. 1999Go) while stimulated by spike arrival at 8,000 excitatory model synapses (rate vexc = 0.3 Hz each) and 2,000 inhibitory model synapses (rate {nu}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 {varepsilon}exc and {varepsilon}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 {Gamma} = 0.82.

We then keep the kernels {varepsilon}exc and {varepsilon}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 ({Gamma} = 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. 2003Go) and a reduction of the membrane time constant (Destexhe and Paré 1999Go; Destexhe et al. 2001Go). 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 {kappa} 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 {varepsilon}-kernels, which explains the rapid decrease of {Gamma} for low and high values of {nu}inh. In particular, for reduced input rates (e.g., {nu}inh = 1 Hz), the input conductance of the full model is reduced and hence, the membrane time constant is increased. Because the {varepsilon}-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).



View larger version (16K):
[in this window]
[in a new window]
 
FIG. 10. Conductance injection. A: mean firing rate {nu}post of the SRM (solid line with squares) compared to that of the full model (dashed line with circles). Excitatory discharge frequency is kept constant at {nu}exc = 0.3 Hz, whereas we vary the inhibitory discharge frequency between 0.5 and 5.5 Hz. Each plotted symbol gives the mean computed over 5 spike trains of 1 s each with 1 SD (error bars). Arrow highlights the stimulation paradigm used for extraction of {varepsilon}-kernels and optimization of the threshold. B: coincidence factor {Gamma}. All the rest as in A. C: spike train of the SRM model (dashed line) is compared to that of the full fast-spiking neuron model (solid line; {nu}exc = 0.3 Hz and {nu}inh = 2.5 Hz). D: normalized histogram of the error in prediction of the subthreshold membrane voltage bin after bin for the same spike train as in C (light gray area). Fitted Gaussian distribution with {chi}2 = 0.001; mean = –1.3 mV, and SD = 3.7 mV (solid line).

 

 DISCUSSION
 
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
In this paper, we mainly focused on deterministic single-variable IF models (i.e., a spike is elicited if a single variable u crosses a threshold {vartheta} 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. 2003Go; FitzHugh 1961Go; Nagumo et al. 1962Go) 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 step-by-step procedure for numerical fitting. An extension toward neurons with hyperpolarization-activated 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. 2003Go; Brillinger 1988Go; Brillinger and Segundo 1979Go; Gerstner 2001Go; Tuckwell 1988Go). 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 1976Go; Mainen and Sejnowski 1995Go).

Reduction of neuronal complexity to a single-variable 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 2003Go). 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 random-current 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 {Gamma} defined in Eq. 22. An alternative approach could be to derive a quality measure from information theory as proposed recently (Arcas et al. 2003Go). 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 single-variable 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 {Gamma} significantly below one. The performance of the MCIF model could be improved by replacing the fixed reset values mreset, hreset, n1,reset, n2,reset, and ureset by values that depend on the gating variables m, h, n1, and n2 at time (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 single-variable IF models.

The one-dimensional NLIF model that we used in this paper is very similar to the exponential IF model proposed recently (Fourcaud-Trocmé et al. 2003Go). A direct comparison of the NLIF model with the SRM shows that both yield a similar performance. The time-dependent 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 ({Gamma} ≥ 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.

Single-variable 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 time-dependent 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 time-dependent 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 time-dependent input current are available.


 APPENDIX
 
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
Optimal filtering problem

In this appendix, we adapt the classic Wiener–Hopf optimal filtering procedure (Wiener 1958Go) to our problem. We start with the usual cross-correlation method to find the first-order Wiener kernel (Lee and Schetzen 1965Go). 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 1988Go; Brillinger and Segundo 1979Go).

General case.

The classic formulation of the Wiener–Hopf optimal filtering problem is the following. Let us suppose that we record the output signal Vtdata from some device driven with an input signal It. 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 Vtdata. 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 cross-correlation 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 time-dependent in the SRM formalism (the so-called kernel {kappa}). More precisely, it depends on the time elapsed since the last emitted spike. We therefore replace Fk with Ft,k where 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 ti, i = 1, ... , T

(A8)
where s is simply 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 cross-correlation between 2 vectors. If the input I is a stationary random process, X[I, I]k {propto} 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 F1 and F2 and that both filters are causal (M = 0). Again, we want to find F1 and F2 so that Vtmodel is as close as possible from Vtdata. Let us define F and It 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 F1 and F2 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 cross-correlations between inputs I1 and I2.

For an application to our case, we need again to take into account the fact that the optimal filters F1 and F2 are time-dependent in the SRM formalism (the so-called kernels {varepsilon}). More precisely, they depend on the time elapsed since the last emitted spike. We follow the same way as above and replace both Fk1 and Fk2 by respectively Ft–,k1 and Ft–,k2, where 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 Vtmodel can be restated as a scalar product of the 2 vectors Ft– and It

(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 i, i = 1,... , T for the total error function

(A25)
where s is short for t . The optimal filters Fs1 and Fs2 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
 
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
This work was supported by Swiss National Science Foundation (SNF) Grant N FN-2100-065268.


 ACKNOWLEDGMENTS
 
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
We thank M. Richardson for many helpful comments and discussions.


 FOOTNOTES
 
* All authors contributed equally to this work. Back

Address for reprint requests and other correspondence: R. Jolivet, Laboratory of Computational Neuroscience, EPFL, 1015 Lausanne, Switzerland (E-mail: renaud.jolivet{at}epfl.ch).

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. Back


 REFERENCES
 
 TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS
 DISCUSSION
 APPENDIX
 GRANTS
 ACKNOWLEDGMENTS
 REFERENCES
 
Abbott L and Kepler T. Model neurons: from Hodgkin–Huxley to Hopfield. In: Statistical Mechanics of Neural Networks, edited by Garrido L. Berlin: Springer-Verlag, 1990.

Abbott L and van Vreeswijk C. Asynchronous states in a network of pulse-coupled oscillators. Phys Rev E 48: 1483–1490, 1993.[CrossRef]

Arcas B, Fairhall A, and Bialek W. Computation in a single neuron: Hodkin and Huxley revisited. Neural Comp 15: 1715–1749, 2003.[CrossRef][Web of Science][Medline]

Benda J and Herz A. A universal model for spike-frequency adaptation. Neural Comp 15: 2523–2564, 2003.[CrossRef][Web of Science][Medline]

Bower J and Beeman D. The Book of Genesis. New York: Springer-Verlag, 1995.

Brette R and Guignon E. Reliability of spike timing is a general property of spiking model neurons. Neural Comp 15: 279–308, 2003.[CrossRef][Web of Science][Medline]

Brillinger D. Maximum likelihood analysis of spike trains of interacting nerve cells. Biol Cybern 59: 189–200, 1988.[CrossRef][Web of Science][Medline]

Brillinger D and Segundo J. Empirical examination of the threshold model of neuronal firing. Biol Cybern 35: 213–220, 1979.[CrossRef][Web of Science][Medline]

Bryant H and Segundo J. Spike inititation by transmembrane current: a white noise analysis. J Physiol 260: 279–314, 1976.[Abstract/Free Full Text]

Calvin W and Stevens C. Synaptic noise and other sources of randomness in motoneuron interspike intervals. J Neurophysiol 31: 574–587, 1968.[Free Full Text]

Dayan P and Abbott L. Theoretical Neuroscience. Cambridge, MA: MIT Press, 2001.

Destexhe A. Conductance-based integrate-and-fire models. Neural Comp 9: 503–514, 1997.[CrossRef][Web of Science][Medline]

Destexhe A and Paré D. Impact of network activity on the integrative properties of neocortical pyramidal neurons in vivo. J Neurophysiol 81: 1531–1547, 1999.[Abstract/Free Full Text]

Destexhe A, Rudolph M, Fellous J, and Sejnowski T. Fluctuating synaptic conductances recreate in vivo-like activity in neocortical neurons. Neuroscience 107: 13–24, 2001.[CrossRef][Web of Science][Medline]

Erisir A, Lau D, Rudy B, and Leonard C. Specific K+ channels are required to sustain high frequency firing in fast-spiking neocortical interneurons. J Neurophysiol 82: 2476–2489, 1999.[Abstract/Free Full Text]

Ermentrout B. Type I membranes, phase resetting curves, and synchrony. Neural Comp 8: 979–1001, 1996.[Web of Science][Medline]

Ermentrout B and Kopell N. Parabolic bursting in an excitable system coupled with a slow oscillation. SIAM J Appl Math 46: 233–253, 1986.

Feng J. Is the integrate-and-fire model good enough—a review. Neural Comp 13: 955–975, 2001.

FitzHugh R. Impulses and physiological states in models of nerve membrane. Biophys J 1: 445–466, 1961.

Fourcaud-Trocmé N, Hansel D, van Vreeswijk C, and Brunel N. How spike generation mechanisms determine the neuronal response to fluctuating input. J Neurosci 23: 11628–11640, 2003.[Abstract/Free Full Text]

Fuortes M and Mantegazzini F. Interpretation of the repetitive firing of nerve cells. J Gen Physiol 45: 1163–1179, 1962.[Abstract/Free Full Text]

Geisler C and Goldberg J. A stochastic model of repetitive activity of neurons. Biophys J 6: 53–69, 1966.

Gerstner W. Time structure of the activity in neural network models. Phys Rev E 51: 738–758, 1995.[CrossRef]

Gerstner W. A framework for spiking neuron models—the spike response model. In: Handbook of Biological Physics, edited by Moss F and Gielen S. Amsterdam: Elsevier, 2001, vol. 4, chap. 12, p. 469–516.

Gerstner W and Kistler W. Spiking Neurons Models: Single Neurons, Populations, Plasticity. Cambridge UK: Cambridge Univ. Press, 2002.

Hansel D and Mato G. Existence and stability of persistent states in large neuronal networks. Phys Rev Lett 86: 4175–4178, 2001.[CrossRef][Web of Science][Medline]

Hansel D and Mato G. Asynchronous states and the emergence of synchrony in large networks of interacting excitatory and inhibitory. Neural Comp 15: 1–56, 2003.[CrossRef][Web of Science][Medline]

Harsch A and Robinson H. Postsynaptic variability of firing in rat cortical neurons: the roles of input synchronization and synaptic NMDA receptor conductance. J Neurosci 20: 6181–6192, 2000.[Abstract/Free Full Text]

Hill A. Excitation and accomodation in nerve. Proc R Soc Lond B Biol Sci 119: 305–355, 1936.[Free Full Text]

Hodgkin A. The local electric changes associated with repetitive action in a non-medullated axon. J Physiol 107: 165–181, 1948.

Hodgkin A and Huxley A. A quantitative description of ion currents and its applications to conduction and excitation in nerve membranes. J Physiol 117: 500–544, 1952.[Free Full Text]

Holden A. Models of the stochastic activity of neurons. In: Lecture Notes in Biomathematics, vol. 12. New York: Springer-Verlag, 1976.

Hoppensteadt F and Izhikevich E. Weakly Connected Neural Networks. New York: Springer-Verlag, 1997.

Keat J, Reinagel P, Reid R, and Meister M. Predicting every spike: a model for the responses of visual neurons. Neuron 30: 803–817, 2001.[CrossRef][Web of Science][Medline]

Kepler T, Abbott L, and Marder E. Reduction of conductance-based neuron models. Biol Cybern 66: 381–387, 1992.[CrossRef][Web of Science][Medline]

Kistler W, Gerstner W, and van Hemmen J. Reduction of Hodgkin–Huxley equations to a single-variable threshold model. Neural Comp 9: 1015–1045, 1997.[Web of Science]

Koch C and Segev I. Methods in Neuronal Modeling. Cambridge, MA: MIT Press, 1989.

Latham P, Richmond B, Nelson P, and Nirenberg S. Intrinsic dynamics in neuronal networks. I. Theory. J Neurophysiol 83: 808–827, 2000.[Abstract/Free Full Text]

Lee Y and Schetzen M. Measurement of the wiener kernels of a non-linear system by cross-correlation. Int J Control 2: 237–254, 1965.

Maass W and Bishop C. Pulsed Neural Networks. Cambridge, MA: MIT Press, 1998.

Mainen Z and Sejnowski T. Reliability of spike timing in neocortical neurons. Science 268: 1503–1506, 1995.[Abstract/Free Full Text]

Monier C, Chavane F, Baudot P, Graham L, and Frégnac Y. Orientation and direction selectivity of synaptic inputs in visual cortical neurons: a diversity of combinations produces spike tuning. Neuron 37: 663–680, 2003.[CrossRef][Web of Science][Medline]

Nagumo J, Arimoto S, and Yoshizawa S. An active pulse transmission line simulating nerve axon. Proc IRE 50: 2061–2070, 1962.[CrossRef]

Nelder J and Mead R. A simplex method for function minimization. Comput J 7: 308–313, 1965.

Rauch A, Camera G, Lüscher H, Senn W, and Fusi S. Neocortical pyramidal cells respond as integrate-and-fire neurons to in-vivo–like input currents. J Neurophysiol 90: 1598–1612, 2003.[Abstract/Free Full Text]

Reinagel P and Reid R. Precise firing events are conserved across neurons. J Neurosci 22: 6837–6841, 2002.[Abstract/Free Full Text]

Reyes A. Synchrony-dependent propagation of firing rate in iteratively constructed networks in vitro. Nat Neurosci 6: 593–599, 2003.[CrossRef][Web of Science][Medline]

Richardson MJE. The effects of synaptic conductance on the voltage distribution and firing rate of spiking neurons. Phys Rev E 69: 044405, 2004.

Rinzel J. Excitation dynamics: insights from simplified membrane models. Theor Trends Neurosci Fed Proc 44: 2944–2946, 1985.

Robinson H and Kawai N. Injection of digitally synthesized synaptic conductance transients to measure the integrative properties of neurons. J Neurosci Methods 49: 157–165, 1993.[CrossRef][Web of Science][Medline]

Stein R. A theoretical analysis of neuronal variability. Biophys J 5: 173–194, 1965.

Stein R. Some models of neuronal variability. Biophys J 7: 37–68, 1967.

Stevens C and Zador A. Novel integrate-and-fire like model of repetitive firing in cortical neurons. In: Proc. 5th Joint Symp Neural Comput 1998 [can be downloaded from: http://www.sloan.salk.edu/~zador/publications.html].

Tanabe S and Pakdaman K. Noise-enhanced neuronal reliability. Phys Rev E 64: 041904, 2001.[CrossRef]

Tuckwell H. Introduction to Theoretic Neurobiology. Cambridge, MA: Cambridge Univ. Press, 1988.

Wang X. Calcium coding and adaptative temporal computation in cortical pyramidal neurons. J Neurophysiol 79: 1549–1566, 1998.[Abstract/Free Full Text]

Wang Y, Gupta A, Toledo-Rodriguez M, Wu CZ, and Markram H. Anatomical, physiological, molecular and circuit properties of nest basket cells in the developing somatosensory cortex. Cereb Cortex 12: 395–410, 2002.[Abstract/Free Full Text]

Wehmeier U, Dong D, Koch C, and van Essen D. Modeling the mammalian visual system. In: Methods in Neuronal Modeling. Cambridge, MA: MIT Press, 1989, p. 335–359.

Weiss T. A model of the peripheral auditory system. Kybernetik 3: 153–175, 1966.[CrossRef][Web of Science][Medline]

Wiener N. Nonlinear Problems in Random Theory. Cambridge, MA: MIT Press, 1958.




This article has been cited by other articles:


Home page
J. Neurophysiol.Home page
L. Badel, S. Lefort, R. Brette, C. C. H. Petersen, W. Gerstner, and M. J. E. Richardson
Dynamic I-V Curves Are Reliable Predictors of Naturalistic Pyramidal-Neuron Voltage Traces
J Neurophysiol, February 1, 2008; 99(2): 656 - 666.
[Abstract] [Full Text] [PDF]


Home page
J. Neurophysiol.Home page
B. Doiron, A.-M. M. Oswald, and L. Maler
Interval Coding. II. Dendrite-Dependent Mechanisms
J Neurophysiol, April 1, 2007; 97(4): 2744 - 2757.
[Abstract] [Full Text] [PDF]


Home page
J. Neurosci.Home page
J. G. Mancilla, T. J. Lewis, D. J. Pinto, J. Rinzel, and B. W. Connors
Synchronization of Electrically Coupled Pairs of Inhibitory Interneurons in Neocortex
J. Neurosci., February 21, 2007; 27(8): 2058 - 2073.
[Abstract] [Full Text] [PDF]


Home page
J. Neurophysiol.Home page
G. La Camera, A. Rauch, D. Thurbon, H.-R. Luscher, W. Senn, and S. Fusi
Multiple Time Scales of Temporal Response in Pyramidal and Fast Spiking Cortical Neurons
J Neurophysiol, December 1, 2006; 96(6): 3448 - 3464.
[Abstract] [Full Text] [PDF]


Home page
J. Neurophysiol.Home page
Q. J. M. Huys, M. B. Ahrens, and L. Paninski
Efficient Estimation of Detailed Single-Neuron Models
J Neurophysiol, August 1, 2006; 96(2): 872 - 890.
[Abstract] [Full Text] [PDF]


Home page
J. Neurosci.Home page
J. W. Pillow, L. Paninski, V. J. Uzzell, E. P. Simoncelli, and E. J. Chichilnisky
Prediction and Decoding of Retinal Ganglion Cell Responses with a Probabilistic Spiking Model
J. Neurosci., November 23, 2005; 25(47): 11003 - 11013.
[Abstract] [Full Text] [PDF]


Home page
J. Neurophysiol.Home page
R. Brette and W. Gerstner
Adaptive Exponential Integrate-and-Fire Model as an Effective Description of Neuronal Activity
J Neurophysiol, November 1, 2005; 94(5): 3637 - 3642.
[Abstract] [Full Text] [PDF]


Home page
J. Neurosci.Home page
S. J. Slee, M. H. Higgs, A. L. Fairhall, and W. J. Spain
Two-Dimensional Time Coding in the Auditory Brainstem
J. Neurosci., October 26, 2005; 25(43): 9978 - 9988.
[Abstract] [Full Text] [PDF]


Home page
Proc. Natl. Acad. Sci. USAHome page
T. Toyoizumi, J.-P. Pfister, K. Aihara, and W. Gerstner
Generalized Bienenstock-Cooper-Munro rule for spiking neurons that maximizes information transmission
PNAS, April 5, 2005; 102(14): 5239 - 5244.
[Abstract] [Full Text] [PDF]


Home page
J. Neurophysiol.Home page
M. Giugliano, P. Darbon, M. Arsiero, H.-R. Luscher, and J. Streit
Single-Neuron Discharge Properties and Network Activity in Dissociated Cultures of Neocortex
J Neurophysiol, August 1, 2004; 92(2): 977 - 996.
[Abstract] [Full Text] [PDF]


This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Supplemental Data
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Web of Science (41)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Jolivet, R.
Right arrow Articles by Gerstner, W.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Jolivet, R.
Right arrow Articles by Gerstner, W.


HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
Visit Other APS Journals Online
Copyright © 2004 by the The American Physiological Society.