| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
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 |
|---|
|
|
|
INTRODUCTION |
|---|
|
), can reproduce electrophysiological measurements to a high degree of accuracy (for a review, see Bower and Beeman 1995; Koch and Segev 1989). Unfortunately, because of their intrinsic complexity, these models are usually difficult to analyze and are computationally expensive in numerical implementations. For this reason, simple phenomenological spiking neuron models such as integrate-and-fire (IF) models (Geisler and Goldberg 1966; Hill 1936; Stein 1965; Tuckwell 1988) are highly popular and have been used to discuss aspects of neural coding, memory, or network dynamics (for a review, see Gerstner and Kistler 2002; Maass and Bishop 1998).
By replacing the rich dynamics of Hodgkin–Huxley-type models by an essentially one-dimensional fire-and-reset process, many details of the electrophysiology of neurons will be missed. In particular, standard leaky IF models do not correctly reproduce neuronal dynamics close to the firing threshold. A systematic reduction of the neuronal dynamics of type I models [i.e., neurons with a smooth frequency–current curve (Hodgkin 1948)] in the limit of very low firing rates yields a canonical type I model (Ermentrout 1996; Ermentrout and Kopell 1986; Hoppensteadt and Izhikevich 1997) that is equivalent to a quadratic IF model (Hansel and Mato 2001; Latham et al. 2000). Although quadratic IF models give, by construction, a correct description of neuronal dynamics close to the firing threshold, it is unclear how well quadratic and other generalized IF models such as the exponential IF model (Fourcaud-Trocmé et al. 2003) or the Spike Response Model (SRM) (Gerstner and Kistler 2002; Kistler et al. 1997) perform for a realistic, time-dependent input scenario where the neuron could spend a significant amount of time far away from the firing threshold. Keat and colleagues (2001) have shown that a phenomenological model of neuronal activity can predict every spike of lateral geniculate nucleus (LGN) neurons with a millisecond precision. However, these neurons produce very stereotyped spike trains with short periods of intense activity followed by long periods of silence (Reinagel and Reid 2002) so that their approach could be limited to LGN neurons only. Recently, it was shown that an IF model can predict the mean rate of pyramidal cells recorded in in vitro experiments (Rauch et al. 2003) over a broad range of different time-dependent inputs. Moreover, the experimental distributions of membrane potentials in the subthreshold regime are well reproduced by leaky IF models (Destexhe et al. 2001). Quadratic IF neurons can approximate the frequency–current curve of a detailed conductance-based model (Hansel and Mato 2003). Finally, it was shown that generalized IF models can approximate the dynamics of the classic Hodgkin–Huxley model of the squid giant axon with high accuracy (Feng 2001; Kistler et al. 1997). The aim of the present paper is 2-fold.
First, we attempt to illustrate the relation of phenomenological spiking neuron models to conductance-based models (Abbott and Kepler 1990; Destexhe 1997; Ermentrout 1996; Ermentrout and Kopell 1986; Kistler et al. 1997; Latham et al. 2000). To do so, we will use a step-by-step analytical derivation of 2 formal spiking neuron models starting from a detailed conductance-based model of a fast-spiking cortical interneuron. In both cases, we compare the behavior of the reduced model to that of the detailed conductance-based model for 3 different input scenarios: a strong isolated current pulse, constant superthreshold input, and random current. Although isolated current pulses and constant drive are standard experimental paradigms, a random current is thought to be more realistic in that it reflects an approximation to the random conductance background caused by input spikes of cortical neurons in vivo (see, e.g., Calvin and Stevens 1968; Destexhe and Paré 1999).
Second, we apply a numerical technique to map the reduced models to the detailed conductance-based model of an interneuron. The optimal parameters characterizing the reduced models are extracted from a sample spike train generated with the conductance-based neuron model by a procedure which generalizes previous approaches (Brillinger 1988; Brillinger and Segundo 1979; Wiener 1958). We compare quantitatively the predictions of the reduced models with that of the full conductance-based model in the random input scenario. We explore the regime of validity of the reduced models by varying the mean and the standard deviation of the input current in a biologically realistic range, and we test the robustness of the reduced models to further simplifications. A generalization of our approach to adapting neurons is indicated. Finally, we extend our numerical technique to a conductance injection scenario, which has become increasingly popular in recent years (Destexhe and Paré 1999; Destexhe et al. 2001; Harsch and Robinson 2000; Reyes 2003).
|
|
METHODS |
|---|
|
In formal spiking neuron models, action potentials are generated by a threshold process. The neuron fires whenever the variable u reaches a threshold 
from below
 | (1) |
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) |
m as the membrane time constant, R as the input resistance, and ueq as the equilibrium potential of the leakage conductance. Integration of Eq. 2 yields the membrane potential as a function of time. Firing is defined by the threshold condition (Eq. 1). After firing, the membrane potential is reset to a value ureset. The LIF neuron may also incorporate an absolute refractory period, in which case we proceed as follows. If u reaches the threshold at time t = , we interrupt the dynamics (Eq. 2) during an absolute refractory time 
refr and restart the integration at time + refr with the new initial condition ureset. We emphasize that firing is a formal event defined by the threshold condition (Eq. 1). The time course of the action potential is disregarded in standard IF models and only the firing time is recorded. In a general nonlinear integrate-and-fire model (NLIF model), Eq. 2 is replaced by (Abbott and van Vreeswijk 1993)
 | (3) |
and reinitialized at u = ureset. A specific instance of a NLIF model is the quadratic model (Feng 2001; Hansel and Mato 2001; Latham et al. 2000)
 | (4) |
. For Iext = 0 and initial conditions u < uc, the voltage decays to the resting potential ueq. For u > uc, the voltage increases up to where an action potential is triggered; uc can therefore be interpreted as the critical voltage for spike initiation by a short current pulse. The quadratic IF model is closely related to the 
-neuron, a canonical type I neuron model (Ermentrout 1996; Hoppensteadt and Izhikevich 1997). Spike response model
In contrast to the NLIF model, the LIF model defined in Eq. 2 can be analytically integrated for arbitrary time-dependent input Iext(t). Let us denote the last firing time by . For t > + refr, the membrane potential is
 | (5) |
0 and zero otherwise; ureset is the initial condition of the integration at time + refr. If we replace the exponentials multiplied by a step function by "response kernels" 
and 
, we may rewrite Eq. 5 in the form
 | (6) |
; Kistler et al. 1997). The kernel models, in the general case, the spike itself and the afterhyperpolarization that follows the spike. In the specific case of Eq. 5, describes the hard reset of u(t) to ureset at time t = + refr. The kernel describes the response of the membrane potential to an input current pulse. In case of the LIF neuron, both kernels are characterized by an exponential decay with time constant m (compare Eq. 6 and 5).
Although both the NLIF model (Eq. 3) and the SRM (Eq. 6) contain the LIF model (Eq. 2) as a special case, the direction of the generalizations is somewhat different. In the NLIF model, parameters are made voltage dependent, whereas in the SRM they depend on t – (i.e., the time since the last spike). To illustrate the relation, compare the NLIF models in Eq. 3 and 4 to an IF model with a time-dependent time constant (i.e., a SRM; Stevens and Zador 1998; Wehmeier et al. 1989)
 | (7) |
m/R. Starting the integration at u( + refr) = ureset, we find for t > + refr
 | (8) |
As a further generalization, we replace in the SRM the fixed threshold by a dynamic one (Fuortes and Mantegazzini 1962; Geisler and Goldberg 1966; Holden 1976; Stein 1967; Weiss 1966)
 | (9) |
refr, we may, for example, set to a large and positive value to avoid firing and let it relax back to its equilibrium value for t > + refr.
If input is provided by the activation of a conductance-based synapse (e.g., attributed to presynaptic spike arrival), rather than current injection, a different formulation of the state of the neuron is provided. The convolution with the response kernel is replaced by a sum and a response kernel 
that takes into account the filtering ascribed to the synapses
 | (10) |
can be interpreted as the postsynaptic potential generated by each spike and it also depends on the time elapsed since the last emitted spike (see above). This latter equation can be transformed back into a convolution product. Let us define the sequence of presynaptic spikes by S(t) = 
j f (t – tj(f)). Equation 10 can thus be restated
 | (11) |
Full conductance-based model
As our reference model, we take the conductance-based model of a fast-spiking cortical interneuron proposed by Erisir and coworkers (1999). We have chosen this specific model so that, because fast-spiking neurons show little adaptation, we avoid most of the complications caused by slow ionic processes that would not be well captured by the class of idealized spiking neuron models reviewed above. In fact, because the model of a fast-spiking neuron is comparatively simple, we can hope to illustrate the steps necessary for a reduction to spiking models in a transparent fashion.
The fast-spiking model neuron (Erisir et al. 1999) follows the Hodgkin–Huxley formalism and consists of a single homogeneous compartment with a nonspecific leak current and 3 active ionic membrane currents
 | (12) |
 | (13) |
 | (14) |
 | (15) |
). Finally, there is the linear leak current
 | (16) |
m = C/gI. Each gating variable follows the equation
 | (17) |
so as to reproduce the behavior of fast-spiking neocortical interneurons. We adopt their set of parameters, except for the leak current where we have chosen gl = 0.25 mS/cm2 (instead of their value of 1.25 mS/cm2) so as to set the passive membrane time constant to about 4 ms. This larger value of the passive membrane time constant seems to be realistic for a broad class of interneurons (Wang et al. 2002). With this set of parameters, the Erisir model exhibits the typical behavior of a type I neuron model, that is, a continuous frequency–current curve that allows firing at very low frequencies and delayed action potentials for pulse input that is just superthreshold.
|
Current injection
In the current injection scenario, the neuron is driven with an uncorrelated Gaussian-distributed random current. We vary both the mean µ and the standard deviation 
of the Gaussian current. This kind of highly variable time-dependent input is thought to be more realistic than other current injection stimulation protocols (Calvin and Stevens 1968; Destexhe and Paré 1999).
Conductance injection (stochastic presynaptic spike arrival)
An even more realistic input scenario is to consider stochastic spike arrival at excitatory and inhibitory synapses. Such a scenario is equivalent to driving the neuron with a highly variable random conductance multiplied by a driving force that depends on the instantaneous membrane voltage of the cell. More precisely, the total synaptic input current Isyn is given by (Robinson and Kawai 1993)
 | (18) |
 | (19) |
 | (20) |
with (Dayan and Abbott 2001)
 | (21) |
increases by an amount Aexc for each presynaptic spike arriving at the synapse at time tj,k, and decays with a time constant exc. The result is a low-pass–filtered version of the presynaptic spike train {tj,k}. The dynamics of inhibitory synapses are defined similarly.
Presynaptic spike trains are described by random homogeneous Poisson processes. At each time step, 
independent random variables are generated and distributed among the N > synapses to generate slightly correlated spike trains (Destexhe and Paré 1999). Numerical values are summarized in Table 2.
|
The coincidence factor between two spike trains (Kistler et al. 1997) is defined by
 | (22) |
between the 2 spike trains, and 
Ncoinc
= 2
N1 is the expected number of coincidences generated by a homogeneous Poisson process with the same rate as the spike train S2. In this article, the reference spike train S1 is always generated by the full conductance-based model of a fast-spiking interneuron, whereas S2 is a spike train generated by one of the reduced models (IF models). The factor N = 1 – 2 normalizes to a maximum value of one, which is reached if and only if the spike train of the reduced model reproduces exactly that of the full model. A homogeneous Poisson process with the same number of spikes as the reduced model model would yield = 0. We compute the coincidence factor by comparing the 2 complete spike trains: the spike train S1 generated by the full conductance-based model and the train S2 predicted by the reduced model.Therefore, in this article, gives a measure of the ability of the reduced model to predict the spike train of the full model. Analytical reduction to an IF model
To find a NLIF model that approximates the dynamics of the full fast-spiking neuron model, we proceed in 2 steps. In a first step, we keep all variables, but introduce a threshold for spike initiation. We call this the multicurrent IF (MCIF) model. In a second step, we separate gating variables into fast and slow ones (Abbott and Kepler 1990; Kepler et al. 1992; Rinzel 1985). The fast variables are turned into instantaneous ones, whereas the slow variables are replaced by constants. The result is the desired NLIF model, which depends on only a single variable.
step 1.
For the first step, we make use of the observation that the shape of an action potential of the fast-spiking neuron model is always roughly the same, independently of the way the spike is initiated. Instead of calculating the shape of an action potential again and again, we can therefore simply stop the costly numerical integration of the nonlinear differential equations as soon as a spike is triggered and restart the integration after the downstroke of the spike about 1.5–2 ms later. We call such a scheme a multicurrent IF model. The interval between the spike trigger time and the restart of the integration corresponds to an absolute refractory period refr.
The MCIF model is defined by a voltage threshold , a refractory time refr, and the reset values from which the integration is restarted. We fix the threshold at = –40 mV; the exact value is not critical, and we could take values of –20 to –45 mV without significantly changing the results. With = –40 mV, a refractory time of refr = 1.7 ms is suitable and an appropriate reset value for the voltage variable is ureset = –85 mV.
The first important approximation involves the reset of the gating variables m, h, n1, and n2, because their time courses are not as stereotyped as that of the membrane potential. To illustrate this point, let us focus on the variable h. At 1.7 ms after spike initiation, the variable h is at a value of 0.16 during periodic firing at about 40 Hz; however, it is at a value of 0.26 when only a single action potential is triggered (not shown). Thus, the optimal value of hreset depends on the choice of input scenario. If the biological neuron under consideration has a mean firing rate of 40 Hz, then we should choose a reset value appropriate for the particular regime of periodic firing. If, on the other hand, the biological neuron is near rest most of the time and emits only an occasional spike, we should choose different reset values. In the following, we adjust the reset values based on a scenario with constant drive current Iext = 5 µA/cm2 that leads to a mean firing rate of about 40 Hz. The reset values are mreset = 0.0; hreset = 0.16; n1,reset = 0.874; n2,reset = 0.2; and ureset = –85 mV. This set of parameters leads to a near-perfect fit of the time course of the membrane potential during periodic firing at 40 Hz and approximates the gain function of the full fast-spiking neuron model to a high degree of accuracy (Fig. 1, A and B).
|
We now turn to random input. The amplitude I of the fluctuations determines the mean firing rate of the neuron model. We see from Fig. 1D that the mean rate of the MCIF model follows closely that of the full model. For a more detailed comparison of the MCIF with the full model, we stimulate both models with exactly the same random current (i.e., identical initiation of the random number generators). In Fig. 1E, we see that the voltage time course of the 2 models is indistinguishable most of the time. Occasionally, the MCIF model misses a spike or adds an extra spike. For this specific input scenario (where the input fluctuations have strength 25 µA/cm2), about 95% of the spike times or the full model are reproduced correctly (with a resolution of = ±2 ms) by the MCIF model. We see from Fig. 1F that the coincidence rate varies as a function of the fluctuation amplitude, but stays above 0.75 in the whole range that we considered. As expected, the value of is highest when the neuron fires at an average rate between 35 and 55 Hz (i.e., in the regime for which parameters have been optimized).
step 2.
The above MCIF neuron is not yet the desired NLIF model because it still depends on all 5 variables, u, m, h, n1, and n2. To reduce it further to a single-variable model that depends only on the membrane potential u, we need to consider a gating variable x either as fast compared to u, in which case we replace x by x
(u), or as slow compared to u, in which case we take x as constant (see also Abbott and Kepler 1990; Kepler et al. 1992; Rinzel 1985). In our case, m(t) is the only fast variable (in the subthreshold range [–100, –40] mV, we find 0.08 
m 0.25 ms, 4.28 h 14.45 ms, 44.66 n1 144.10 ms and 0.44 n2 4.19 ms). We eliminate m by replacing m(t) by its equilibrium value m0[u(t)] (see Fig. 2A). The treatment of the other gating variables deserves some extra discussion.
|
refr = 4 ms, most of the excursion trajectory of n2 falls within the refractory period. Between spikes we can therefore replace n2 by its equilibrium value at rest n2 
n2,eq = n0(ueq). The gating variables h and n1 vary only slowly so that, for a given input scenario, the variables may be replaced by their average value hav and n1,av. The average, however, does depend on the input scenario. To illustrate this point, we consider 3 different situations. First, for a neuron at rest that receives only occasionally an isolated current pulse, the average values are close to the resting values. In other words, we may use the parameters at rest (i.e., hav = 0.87 and n1,av = 0.00057). If the neuron is stimulated by a constant bias current that shifts the membrane potential just below threshold, we should take hav = 0.54 and n1,av = 0.019. Finally, in a periodic firing regime (constant stimulation with I = 5 µA/cm2) reasonable averages are hav = 0.45 and n1,av = 0.8 (see Fig. 2, C and D). The last pair of values will be used in the following as our standard set of parameters unless stated otherwise.
With m = m0(u) and constant values for h, n1, and n2, the dynamics of the full fast-spiking neuron model defined in Eqs. 12–17 reduces to
 | (23) |
 | (24) |
= –1/(dF/du) where the derivative is to be evaluated at u = ueq. In principle the function F could be further approximated by a linear function with slope –1/ and then combined with a threshold at, say, = –45 mV. This would yield a standard linear LIF model. Alternatively, F could be approximated by a quadratic function that would lead us to the quadratic IF neuron. The canonical type I model would correspond to a quadratic IF neuron where parameters are optimized in the limit of a very low firing frequency (Ermentrout 1996; Ermentrout and Kopell 1986; Hansel and Mato 2001; Latham et al. 2000).
In the following, we do not make any further approximations. Instead we work directly with the nonlinear function F(u). The refractory period refr = 4 has already been introduced above. To finish the definition of the model, we have to specify the threshold and the reset value. We take = –45 mV (the exact value is not critical). After firing, integration is restarted at ureset = –85 mV.
Analytical reduction to a spike response model
For a fitting of the full fast-spiking neuron model defined in Eqs. 12–17 to the SRM defined in Eqs. 6 and 9, we need to determine the kernels (t – ) and (t – , s), and furthermore the (time-dependent) threshold (t – ) must be adjusted. As a first step, we stimulate the full model by a short suprathreshold current pulse to determine the time course of the action potential and afterhyperpolarization. Let us define as the time when the membrane potential crosses an (arbitrarily set) threshold (e.g., = –50 mV). The kernel (t – ) is defined by the time course of the membrane potential for t > (i.e. during and after the action potential) in the absence of external input. If we were interested in a purely phenomenological model, we could simply record the numerical time course u(t) and define (t – ) = u() for t > (see Numerical optimization method). Instead of such a purely numerical method, it is instructive to take a semianalytical approach and study the 4 gating variables m, h, n1, and n2 immediately after action potential generation. About 2 ms after initiation of a spike, all 4 variables have passed their maximum or minimal values and are on their way back to equilibrium (cf. Fig. 3, A and B). We set refr = 2 ms. For t + refr, similar to the approach of Destexhe (1997), we fit the approach to equilibrium by an exponential
 | (25) |
x is a fixed time constant and xreset is the initial condition at t = + refr, and xeq = x (ueq). In other words, the voltage-dependent differential Eq. 17 for x is replaced by the linear voltage-independent equation
 | (26) |
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.
|
|
 | (27) |
) is an exponential function with time constant n2/2. We insert the time-dependent conductance into Eq. 12 and find for t + refr
 | (28) |
(t – ) = C/j gj(t – ) and Iion(t – ) = jgj(t – ) Ej, we arrive at
 | (29) |
.
According to Eqs. 27 and 28, the ion currents have, at time + refr, a value of Iion(refr) = j gj(refr)Ej. For t – = refr, the effective time constant (refr) of the membrane voltage is extremely short (Fig. 3C) so that we may assume that it is at its steady state value given Iion(refr). For t + refr, we therefore integrate Eq. 29 with the initial condition
 | (30) |
– refr with the kernels given by
 | (31) |
 | (32) |
The kernels and that we have constructed so far are limited to t > + refr. During the absolute refractory period < t < + refr the neuron is not responsive to external input. We therefore set to zero. The action potential itself that occurs in the interval < t < + refr is, for the sake of simplicity, approximated by a triangular voltage pulse. Finally, we introduce a dynamic threshold
 | (33) |
0, 1, refr, and . During the absolute refractory period refr, the threshold is set to a value refr = 100 mV that is sufficiently high to prevent the neuron from firing. After refractoriness, the threshold starts at 0 + 1 and relaxes with a time constant of to an asymptotic value of 0. The initial value 0 + 1 was arbitrarily set at 0 mV. 0 = –50 mV and = 6 ms were then chosen so that for 2 different inputs (i.e., I0 = 5 µA/cm2 and I0 = 30 µA/cm2) the mean firing rate of the SRM was approximately correct. Numerical optimization method
In the SRM framework, 3 quantities are needed to define the model. These are the kernel , describing the shape of a spike; the kernel , describing the input integration process and finally the threshold for spike initiation. The previous section described a semianalytical method to reduce the full conductance-based model to a SRM, obtaining approximations of the kernels. In this section, we describe a numerical method that finds kernels that optimize the fit of the SRM to a spike train. The method proceeds in 3 steps. First, we use the fact that spikes have always roughly the same shape to extract from the course of the membrane voltage. Second, we use a linear approximation for the integration of the input to extract as the best linear filter between the driving current and the fluctuations of the membrane potential in the subthreshold regime. Third, we fit the parameters of the threshold to find the best reproduction of the specific spike train that we use. The resulting fit is finally evaluated on a set of new spike trains that have not been used during optimization. We now discuss the 3 steps in more detail.
Throughout this subsection, we suppose that we have at disposal a simultaneous recording of the discretized membrane voltage utdata and the discretized driving current Itext of a cell. We assume that both utdata and Itext are recorded at the same sampling frequency.
1) To extract the discretized version of the kernel from the data, we align all the spikes. The spike triggered average of the membrane potential (i.e., the "mean shape" of an action potential) yields, after smoothing, the kernel (see Fig. 4, B and C). Detection and alignment are realized relative to an arbitrary threshold condition on the first time derivative of the membrane voltage (see Fig. 4A).
|
, we move to the extraction of . We start by discretizing the SRM equations. The equivalent of Eq. 6 in discrete time is
 | (34) |
 | (35) |
1 are constants and 
1 are independent random variables drawn from a normal distribution with mean 0 and unit variance. For the sake of simplicity, we set µI = 0 for the rest of this discussion. However, the method can easily be generalized. Equation 34 can therefore be rewritten
 | (36) |
-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) |
,tdata
 | (38) |
The right-hand side of Eq. 37 can be interpreted as a numerical convolution between It and a family of filters parameterized by t – . The problem of finding the best linear filter (t–,s) of some device given the output vector (Vt–,sdata) and the input vector (It) is known in the literature as the Wiener–Hopf optimal filtering problem. However, the classic formulation of this problem is not well suited for the present question because of the explicit dependency of the filter on t – . In the APPENDIX, we show that the Wiener–Hopf approach can nevertheless be adapted to the present situation. The formulation we propose circumvents a couple of technical issues of the classic formulation. In particular, there is no need to choose a specific time window in which to perform the extraction of the filter. Instead, all segments are aligned for a given t – . We find that t–,s is a solution of the following linear system
 | (39) |
t,s is nonzero, t = t – , and X[·, ·] is defined by
 | (40) |
t,s for each t and s by solving Eq. 39. To smooth the results, we fit them with a suitable function. In the present article, we use a single exponential decay in the variable s (see Fig. 4D). It is then possible to plot the time constant of this function versus the delay t and fit the dependency on t by the function (t) = 
1(1 + tanh (2(t – 3)]) with free parameters 1, 2, and 3. The SRM can then be turned back into a differential equation with a time-dependent time constant (Eq. 7; see also DISCUSSION), which is very efficient for simulations.
3) Finally, we choose a specific threshold condition and optimize it in terms of spike train reproduction. We take a dynamic threshold defined by Eq. 33, where refr is always kept at 2 ms and refr at a value of 100 mV to prevent firing. 0, 1, and 
are free parameters. To find the best fit for these 3 parameters, we simulate the full spike train with the SRM and compute 1 – . For a definition of , see Eq. 22. Then, we use the downhill simplex method (Nelder and Mead 1965) algorithm to find a set of parameters that minimize 1 – .
For this article, we have used for the parameter optimization a single spike train of 10 s length generated with Gaussian-distributed random current (cf. Fig. 4F). Parameters of the current are µI = 0 µA/cm2 and 1 = 25 µA/cm2. The value of the driving current is changed every 0.2 ms. Both the membrane voltage and the driving current are sampled at 5-kHz sampling frequency. For this set of parameters, the mean rate is 34.6 Hz.
In case of conductance injection protocol (stochastic spike arrival), steps 1 and 3 of the method are left unchanged, whereas step 2 is slightly modified. The state of the neuron in discrete time is given by (see Eq. 10)
 | (41) |
-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) |
exc and inh both at the same time (see APPENDIX). The filters exc and inh are fitted by double exponentials (see Fig. 9). The kernel was maintained from current injection simulations. The rest of parameter optimization was done on a single spike train of 1 s length generated with discharge frequencies exc = 0.3 Hz and inh = 2.5 Hz. Both the membrane voltage and the driving current are sampled at 5 kHz sampling frequency. In this specific spike train, the mean rate is 29.0 Hz. Finally, Sexc (Sinh) is the poststimulus time histogram PSTH computed over all excitatory (inhibitory) synapses with a given time binning (here, each bin equals 0.2 ms). See also Eq. 11.
|
|
|
RESULTS |
|---|
|
), 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. 1999) by a NLIF with a single variable u, given by the differential equation
 | (43) |
= –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.
|
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 1968; Destexhe and Paré 1999). To choose a function F that yields a good performance of the NLIF model for random input, we take advantage of the fact that F is characterized by 3 constants 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 1965) using a fluctuating input current of zero mean and variance 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 that measures the similarity of the 2 spike trains with a temporal resolution of 2 ms (see METHODS). The parameters found by the numerical optimization are given in DATA SUPPLEMENTS, Table 1.A sample spike train is shown in Fig. 5D. Most of the spikes coincide, whereas some spikes are missed and other spikes are added. Action potentials of the NLIF model are indicated by triangular pulses that span a refractory period of refr = 2 ms; the performance does not change significantly if we take refr in the range of 1 to 4 ms.
We now turn to a systematic exploration of the performances of the numerically optimized NLIF for fluctuating input with different mean µI and variance I. In Fig. 5E we plot the mean firing rate as a function of the fluctuation amplitude for both the full model and the NLIF model. The NLIF model gives a fair reproduction of the mean rates of the full model except for large variance [I > 50 (µA/cm2)] of the input current. For a more detailed comparison we plot the coincidence rate as a function of both the mean and the variance of the stimulating current (Fig. 5G). The reproduction of spike trains on a spike-by-spike basis is good (coincidence rate > 0.7) over a broad range of stimulation parameters. Thus, even though the function F(u) has been optimized using data for fixed µI = 0 and I = 25 µA/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 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., , 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. 1999) may also be approximated by a single differential equation
 | (44) |
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" is a function of t – . Integration of Eq. 44 yields the equation of the SRM
 | (45) |
and , are given by exponential functions (see METHODS for details). If u(t) reaches a dynamic threshold (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 0, 1, and of the dynamic threshold (see METHODS, Eq. 33) are free parameters that we now adjust so as to optimize the performance for a given stimulation paradigm.
For pulse input of 2-ms duration, the model neuron emits at most a single spike so that only the asymptotic threshold 0 matters. For a sufficiently strong suprathreshold current pulse, spike initiation and hyperpolarizing spike afterpotential of the full model are reproduced by the SRM to a high degree of accuracy (cf. Fig. 6C). With = –50 mV the SRM produces action potentials whenever the amplitude of the current pulse is above 12.4 µA/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 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.
|
0 = –50 mV fails to approximate the gain function of the full model to a satisfactory degree (cf. Fig. 6B). We therefore use a dynamic threshold (Eq. 33; see METHODS) which approaches a value of 0 = –50 mV with a time constant = 6 ms and starting at a value of 0 + 1 = 0 mV at time + 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 I = 25 µA/cm2. For the numerical optimization of the SRM defined in Eq. 45 we need to determine the spike shape , the input-response filter , as well as the parameters 0, 1, and of the dynamic threshold (see METHODS). The kernel is the average spike shape and (see Fig. 4, B and C). From Fig. 4D, we see that the kernel can be approximated by an exponential (t, s) = exp[–s/(t)] with a time constant (t) shown in Fig. 4E. The threshold parameters are optimized so that most spikes occur with the correct timing (Fig. 4F). The optimal set of threshold parameters is 0 = –53 mV, 1 = 12.7 mV, and = 54 ms (see also DATA SUPPLEMENTS, Table 3). We keep these parameters fixed in the following. For the specific spike train used for parameter optimization, we find a coincidence factor = 0.83.
We now turn to a systematic exploration of the performances of the numerically optimized SRM. A typical spike train is shown in Fig. 6D. To quantify the performance for different inputs, we first keep the variance of the input fixed and change the mean input µI. Figure 6E shows that the mean firing rates as a function of µI of the full model is comparable to that of the SRM. Moreover, the SRM reproduces the spike timing of the full model to a high degree of accuracy over a broad range of different I and µI, as quantified by the coincidence rate (cf. Fig. 6G). In DATA SUPPLEMENTS, Fig. 1, it is furthermore shown that both the interval distribution and the coefficient of variation of the full model are correctly predicted by the simple model indicating that missed or added spikes do not modify significantly the spike pattern. Finally, in the subthreshold regime, the membrane potential of the reduced model approximates that of the full model within ±2 mV (cf. Fig. 6F). To summarize, the SRM predicts the subthreshold fluctuations, the correct number of spikes, most of them (>70%) with the correct timing and with a firing pattern very similar to the one produced by the full conductance-based neuron model.
The amplitude of the fluctuations of the stimulus (here, I) is a crucial factor for the quality of the predictions within the SRM framework. This can be easily understood. First, if the variance of the input signal is large, the amplitude of fluctuations of the membrane voltage is large, which facilitates the emission of spikes at a correct timing with the threshold condition (Eq. 33) (Brette and Guignon 2003; Bryant and Segundo 1976; Mainen and Sejnowski 1995; Tanabe and Pakdaman 2001). Second, when the constant part of the driving current dominates the fluctuations (µI >> I), whenever a spike is missed or added by the SRM, errors propagate further in time and the coincidence factor decreases dramatically.
Comparison with simple threshold models
In the previous subsections, we have seen that both the NLIF and the SRM show a good performance for random input current, if the fluctuation amplitude is large enough. We may therefore wonder whether this is a universal property that would hold for any 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) |
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).
|
(t – ,s) = exp(–s/) for t > + refr with = 4 ms and refr = 2 ms, and [(t – , s)] = 0 otherwise. This is equivalent to Eq. 2 of the LIF model. After firing, the membrane potential is reset to a value of –85 mV. The optimal parameter 0 for the LIF neuron is –47.33 ± 0.08. Figure 7B shows the quality of predictions for the LIF neuron and, as expected, the results are close to the results of the SRM with constant threshold.
Let us now simplify even further. We use neither reset nor spike afterpotential (i.e., we neglect the interaction between action potentials). Moreover, the kernel is replaced by limt–
t–,s. The membrane voltage is then
 | (47) |
,s = 0 exp(–s/) with 
4 ms (i.e., a simple low-pass filter). As before, spikes are triggered whenever the variable u reaches a constant threshold from below. However, in contrast to the LIF model, there is no reset after firing. Although the minimal model 1 still performs well in the neighborhood of the optimization point, decreases quickly as the stimulation parameters change (see DATA SUPPLEMENTS, Table 2). It even goes below 0 for high frequencies and thus performs worse than a homogeneous Poisson neuron model (see definition of in Eq. 22). This is attributed to the fact that for larger stimulation, the firing frequencies of minimal model 1 are too high and spikes occur systematically too early.
We now consider a minimal model 2, which is even simpler: instead of a low-pass filter we take as instantaneous (i.e., it takes a fixed positive value 1 during one time step and is zero thereafter). This is an approximation to a neuron in the high-conductance state where the voltage follows the current quasi instantaneously. In this case, the variable u is given by
 | (48) |
0, which has been optimized on a spike train generated by random input current with zero mean and variance I = 25 µA/cm2. Minimal model 2 performs poorly (see DATA SUPPLEMENTS, Table 2). In particular, even at the point for which the parameter 0 has been optimized, attains a value of only 0.16 ± 0.01 (mean ± SD). The results with the minimal models indicate that the SRM and the NLIF model add features beyond simple threshold models and these features (i.e., spike–spike interaction for the SRM and nonlinear voltage dependency and voltage reset for the NLIF) are necessary if we want a model that works over a broad range of different inputs. Extending the framework to slow processes
Real neurons often exhibit a rich repertoire of ion channels, in particular, slow ion channels that contribute to frequency adaptation, a widespread phenomenon in biological neurons. To take into account such slow processes, we need to leave the framework of single-variable models and introduce an adaptation variable A. We use a simple relaxation dynamics
 | (49) |
adapt. In the absence of spikes A relaxes with time constant 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 2003; Rauch et al. 2003).
To include adaptation into the SRM framework, we make the time constant of the exponential response kernel (Fig. 4E) depend on A. In the case of the Erisir model of 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 is shortened. We computed the response kernel (t, s) (see METHODS) for different stimulation parameters and plotted the parameter 1/ as a function of the output frequency (Fig. 8A). A linear fit 1/(v) = A0 + A1v gives us the parameters A0 and A1 and the identification of the adaptation variable A in Eq. 49 as A = 1/. The free parameter adapt is set to a value of 450 ms. Figure 8, B–D shows the comparison between nonadapting SRM (Fig. 8C) and adapting SRM (Fig. 8D) for current injection with Gaussian white noise. After some time, parameters of the input are switched, leading to an increase in output frequency (Fig. 8B). The nonadapting SRM maintains as tuned for the first part of the stimulation paradigm and thus performs poorly when stimulation parameters are changed (Fig. 8C, from left to right). The performance of the adapting SRM, however, are much better once adaptation has taken place (Fig. 8D). In terms of the coincidence factor , the nonadapting SRM yields = 0.63 when the threshold is optimized for intermediate stimulation (as in the rest of the article; Fig. 8C), whereas the adapting SRM yields = 0.76 (Fig. 8D). The same methodology can also be applied to strongly adapting neuron models, as illustrated in Figure 8, E–H for the case of a 2-compartment neuron model with Ca2+ channels and Ca2+-gated K+ channels responsible for slow adaptation (Wang 1998). Injection of random current in the soma allows us to identify the time constant ; a test with changing input confirms that the extended SRM can follow the adaptation of the full model.
|
A more realistic input scenario than random current injection would be a model of stochastic conductance changes (Destexhe and Paré 1999). We preferred in the previous sections to work with the more common random current input because this allows us to characterize the input in a transparent fashion. Furthermore, for any random input scenario with stationary white noise characteristics, a random conductance scenario can be replaced, to a high degree of accuracy, by a random current scenario (Richardson 2004). Nevertheless, the numerical method exploited in previous sections can be extended to the case of random conductance injection. Instead of optimizing a kernel (t – , s) that describes the linear response of the membrane potential to an input current, we now optimize two kernels exc(t – , t – tj(f)) and inh (t – , t – tj(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 exc as a function of t – tj(f) describes the excitatory postsynaptic potential (EPSP), whereas inh as a function of t – tj(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. 1999) while stimulated by spike arrival at 8,000 excitatory model synapses (rate vexc = 0.3 Hz each) and 2,000 inhibitory model synapses (rate inh = 2.5 Hz each). Best fits for the dynamic threshold parameters are indicated in DATA SUPPLEMENTS, Table 4.
As a result of the optimization procedure we get the shape of the kernels exc and inh shown in Fig. 9, D and E. Both EPSP and IPSP can be well approximated by double exponentials. For this specific spike train, the coincidence factor is = 0.82.
We then keep the kernels exc and inh fixed and turn to a systematic exploration of the performances. Figure 10 shows the quality of predictions for different values of the discharge frequency of inhibitory synapses when holding the discharge frequency of excitatory synapses constant. The mean firing rates do not match as well as for current injection (Fig. 10A) except at the point at which parameters have been optimized. For this set of input parameters, the coincidence rate is good ( = 0.7; see Fig. 10B). The subthreshold behavior of the membrane voltage is also nicely predicted (Fig. 10, C and D). However, conductance injection is known to produce a suppression of subthreshold membrane voltage fluctuations (Monier et al. 2003) and a reduction of the membrane time constant (Destexhe and Paré 1999; Destexhe et al. 2001). We already mentioned that large membrane voltage fluctuations are a key point for correct prediction of the timing of the spikes. The fact that the effective membrane time constant depends on the stimulation paradigm implies that the EPSPs and IPSPs that we compute are optimal only in a limited range of input characteristics. This is an important difference to the current injection scenario where changes of the kernel with the input characteristics are less important. It is clear then that, in the case of conductance injection, our approach is valid only over a limited range of input characteristics with mean conductances close to those used to compute the -kernels, which explains the rapid decrease of for low and high values of inh. In particular, for reduced input rates (e.g., inh = 1 Hz), the input conductance of the full model is reduced and hence, the membrane time constant is increased. Because the -kernels of the SRM do not change their time constant, the mean membrane potential of the SRM and thus the predicted mean firing rate are too low (see Fig. 10A).
|
|
|
DISCUSSION |
|---|
|
from below. If the evolution of u depends on the instantaneous voltage, we arrive at the NLIF model; if it depends on the time since the last spike, we arrive at the SRM. Models with 2 (Arcas et al. 2003; FitzHugh 1961; Nagumo et al. 1962) dimensions or more offer a wide spectrum of behavior but are often more difficult to optimize numerically. We proposed a straightforward extension of the SRM to the case of adapting neurons allowing a 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. 2003; Brillinger 1988; Brillinger and Segundo 1979; Gerstner 2001; Tuckwell 1988). This is important in regimes where the neuron shows a low degree of spike train reproducibility under repeated stimulation with the same signal (Bryant and Segundo 1976; Mainen and Sejnowski 1995).
Reduction of neuronal complexity to a 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 2003). Similarly, the SRM is based on a linearization about a reference state and there is no a priori reason why it should work well far from that reference point. Nevertheless, we found in this paper that the combination of linear summation with a dynamic threshold works for 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 defined in Eq. 22. An alternative approach could be to derive a quality measure from information theory as proposed recently (Arcas et al. 2003). Although information theory provides a systematic theoretical framework, we think that our coincidence measure offers several practical advantages: in particular, it is easy to evaluate and allows a straightforward interpretation.
On our random input scenario, the MCIF model clearly yields the best performance. Although it is straightforward to implement and rapid to simulate, it is difficult to analyze mathematically. Strictly speaking, it does not fall in the class of 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 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. 2003). 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 ( 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 |
|---|
|
In this appendix, we adapt the classic Wiener–Hopf optimal filtering procedure (Wiener 1958) to our problem. We start with the usual cross-correlation method to find the first-order Wiener kernel (Lee and Schetzen 1965). At that point, a simple modification of the error function definition provides us with an equation that is simple to use for the SRM formalism. Note that the classic Wiener–Hopf method was already used in a similar context by Brillinger and Segundo (Brillinger 1988; Brillinger and Segundo 1979).
General case.
The classic formulation of the Wiener–Hopf optimal filtering problem is the following. Let us suppose that we record the output signal 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) |
 | (A2) |
 | (A3) |
 | (A4) |
 | (A5) |
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 ). 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) |
|
| (A7) |
 | (A8) |
. Then, the optimal time dependent linear filter F for the whole spike train is the solution of the following system of equations
 | (A9) |
 | (A10) |
 | (A11) |
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) |
 | (A13) |
 | (A14) |
Equation A12 can now be rewritten as a vector product
 | (A15) |
|
| (A16) |
|
| (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) |
 | (A19) |
 | (A20) |
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 ). 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) |
 | (A22) |
|
| (A23) |
Thus Vtmodel can be restated as a scalar product of the 2 vectors Ft– and It
 | (A24) |
i, i = 1,... , T for the total error function
|
| (A25) |
. The optimal filters Fs1 and Fs2 are then the solution of the following system of equations
 | (A26) |
 | (A27) |
In its matrix form, this equation gives
 | (A28) |
See Eq. A11 for a definition of X[·, ·].
|
|
GRANTS |
|---|
|
|
|
ACKNOWLEDGMENTS |
|---|
|
|
|
FOOTNOTES |
|---|
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.
|
|
REFERENCES |
|---|
|
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.
Calvin W and Stevens C. Synaptic noise and other sources of randomness in motoneuron interspike intervals. J Neurophysiol 31: 574–587, 1968.
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.
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.
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.
Fuortes M and Mantegazzini F. Interpretation of the repetitive firing of nerve cells. J Gen Physiol 45: 1163–1179, 1962.
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.
Hill A. Excitation and accomodation in nerve. Proc R Soc Lond B Biol Sci 119: 305–355, 1936.
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.
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.
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.
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.
Reinagel P and Reid R. Precise firing events are conserved across neurons. J Neurosci 22: 6837–6841, 2002.
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.
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.
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:
|
|
 |
|
  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]
|
|
|
|
|
|
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]
|
|
|
|
 |
|
 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]
|
|
|
|
|
|
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]
|
|
|
|
|
|
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]
|
|
|
|
|
|
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]
|
|
|
|
|
|
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]
|
|
|
|
|
|
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]
|
|
|
|
 |
|
 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]
|
|
|
|
|
|
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]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Visit Other APS Journals Online |
TOP
ABSTRACT
INTRODUCTION
)]} with 
2 = 2.5 x 10–5; mean = 1.0 mV; SD = 2.2 mV).
