An event-based integration scheme for an integrate-and-fire neuron model with exponentially decaying excitatory synaptic currents and double exponential inhibitory synaptic currents has been introduced by Carnevale and Hines. However, the integration scheme imposes nonphysiological constraints on the time constants of the synaptic currents, which hamper its general applicability. This letter addresses this problem in two ways. First, we provide physical arguments demonstrating why these constraints on the time constants can be relaxed. Second, we give a formal proof showing which constraints can be abolished. As part of our formal proof, we introduce the generalized Carnevale-Hines lemma, a new tool for comparing double exponentials as they naturally occur in many cascaded decay systems, including receptor-neurotransmitter dissociation followed by channel closing. Through repeated application of the generalized lemma, we lift most of the original constraints on the time constants. Thus, we show that the Carnevale-Hines integration scheme for the integrate-and-fire model can be employed for simulating a much wider range of neuron and synapse types than was previously thought.
The most salient feature of neurons is their ability to summate synaptic inputs arriving from other neurons and to respond with the generation of an action potential or spike when the membrane potential reaches a certain threshold value. After its generation, a spike will travel down the neuron's axon and serve as input for other neurons. In its most basic form, spike generation is captured by the so-called integrate-and-fire model. This model was first conceived by Lapicque (1907) a hundred years ago. Lapicque modeled the subthreshold behavior of the membrane potential as a capacitance in parallel with a resistor, mimicking basic electrical properties of the cell membrane. At that time, the spike-generating mechanism was not known, and therefore only a phenomenological description of the process could be given. On the basis of electrophysiological experiments, Lapicque assumed that when the membrane potential reached a threshold value, the cell would generate (fire) a spike, after which the membrane potential would reset to resting level (Abbott, 1999; Lapicque, 1907). Integrate-and-fire models are still widely used today in simulations and for the analytical study of neural network dynamics.
The event-based integration scheme that we analyze here was introduced by Carnevale and Hines (Hines & Carnevale, 2002, 2004) for the widely used NEURON simulation environment (Carnevale & Hines, 2006). Event-based integration methods are also used in many other neural simulators (e.g., NEST, XPP and Mvaspike); (Brette et al., 2007). In event-based models, the synaptic coupling between neurons is mediated by events, which are triggered when the membrane potential crosses the spiking threshold. Events are subsequently communicated to the postsynaptic cells, where they initiate a change in the synapse, leading to excitation or inhibition of the postsynaptic cell. Besides the Carnevale-Hines scheme, many other integration schemes have been applied in integrate-and-fire models. These schemes range from purely numerical integration schemes, such as Euler and Runge-Kutta, to numerically exact calculations based on root-finding algorithms for determining when the membrane potential crosses the spiking threshold (Brette, 2006; Hansel, Mato, Meunier, & Neltner. 1998; Morrison, Straube, Plesser, & Diesmann, 2007; Rudolph, & Destexhe, 2006a, 2006b; Tsodyks, Mitkov, & Sompolinsky, 1993). The Carnevale-Hines scheme takes a middle ground between these two extremes. It uses explicit knowledge of the exact solution to determine whether threshold crossing will occur but avoids the expensive explicit calculation of the threshold passage time. Instead, the scheme employs the computationally cheaper Newton iteration to obtain a spike time estimate.
The Carnevale-Hines scheme is developed to solve an integrate-and-fire model, which includes excitatory synapses as exponentially decaying currents and inhibitory synapses as currents obeying a double exponential function. With these currents, it is in principle possible to describe the main class of excitatory synapses (characterized by AMPA receptors) and the main class of inhibitory synapses (characterized by GABA receptors). As mentioned before, the Carnevale-Hines scheme uses Newton iteration to estimate the threshold crossing times (i.e., the events). In their proof of the correctness of the Newton iteration estimate, Carnevale and Hines used the following constraints on the time constants of the synapses and the membrane time constants: (Hines, & Carnevale, 2002, 2004). These constraints on the parameters lead to the loss of many physiological relevant realizations of the conceptual model, as discussed below.
In cortical areas of mammals, the excitatory AMPA currents have a fast rise time—between 0.1 and 0.8 ms—followed by a fast decay of 1 to 3 ms (Hausser & Roth, 1997; Kleppe & Robinson, 1999). The inhibitory GABAergic currents have a rise time between 1 and 2 ms (Bosman, Rosahl, & Brussaard, 2002) and a decay time varying from about 5 to 30 ms (Bosman, Heinen, Spijker, & Brussaard, 2005). Also, the membrane time constants of different cortical neurons vary over a wide range, from close to 5 ms to well over 40 ms (Karube, Kubota, & Kawaguchi, 2004; Tateno & Robinson, 2006; Trevelyan & Jack, 2002; Zhang, 2004). These physiological data show that the decay times for AMPA synapses are similar to or larger than the GABA rise times and, furthermore, that the GABA decay times can be larger than the membrane time constant. However, the original treatment of the Carnevale-Hines integration scheme requires that the excitatory decay time should be smaller than the inhibitory rise time and that the inhibitory decay time should be smaller than the membrane time constant. Consequently, the full physiological range of parameters as found in cortical neurons is not accessible to the standard implementation of the model. These limitations became apparent to us when studying the effects of GABA receptor maturation on microcircuit processing (van Elburg, Bosman, Ridder, Brussaard, & van Ooyen, 2009) and subsequently led us to reexamine the Carnevale-Hines integration scheme.
The aim of our analysis here is to remove the unphysiological constraints on the time constants while keeping the strength of the event-based integration scheme. Our analysis proceeds in four stages. In the first stage, we introduce the basics of the model. In the second stage, we provide an analysis of the physics of the problem. In the third stage, we present the generalized Carnevale-Hines lemma and its proof. In the fourth stage, the actual analysis of the integration scheme is carried out, in which we use the generalized Carnevale-Hines lemma to abolish the unphysiological constraints on the decay times.
2. Model and Integration Scheme
In this section, we expand on the short introduction of the event-based integrate-and-fire model and subsequently discuss the exact solution of the model's subthreshold dynamics. In our treatment of the model and in the exact solution, we essentially follow the original treatment (Hines & Carnevale, 2002, 2004), except that our notation is slightly different to account for the possibility of having different synaptic subtypes.
2.1. Event-Based Synaptic Current Driven Integrate-and-Fire Model.
2.2. Carnevale-Hines Integration Scheme.
The Carnevale-Hines integration scheme has been developed to solve the problem of finding threshold passage times from the exact solution of the model presented in the previous subsection. The scheme cycles through the following steps: in the first step, an event arrives at the neuron; in the second step, the actual membrane potential and synaptic currents at the arrival time of the event are calculated from the exact solution; in the third step, a check takes place on threshold crossing; in the fourth step, the actual event is handled by updating synaptic currents and calculating a new threshold crossing estimate based on these currents; in the fifth step, a self-event, which will arrive at the estimated threshold crossing time, is generated. Let us examine steps 3 and 4 in more detail. In the third step, a test takes place to establish whether the membrane potential m is close to the threshold θ = 1 (i.e., m > θ − ϵ, with ϵ an arbitrary number satisfying 0 < ϵ ≪ θ). If sufficiently close to threshold, then it is assumed that the membrane potential reached threshold, the membrane potential m is reset to zero, and an event is sent to the outgoing synaptic connections of this cell. In the fourth step, the event is handled. If the event was a synaptic activation, then the weight of the synapse is, in line with the model description above, added to either an eν or a jμ. If the event was a self-event related to an estimated threshold crossing time, then the currents obtained from the exact solution are kept unaltered. After this update of the synaptic currents, an estimate is made about when to evaluate the subthreshold solution again. This estimate is based on the time derivative of the membrane potential m′, which is equal to the total of the synaptic and leak currents. If the current is not depolarizing m′ ⩽ 0, then the evaluation of the membrane potential is postponed until the arrival of a new synaptic event. If the current is depolarizing m′>0, a future spike time te is estimated using Newton iteration te = (1 − m)/m′ + t0. This estimated time is used in step 5 to generate a self-event that is sent by the neuron to itself and indicates the latest time to which evaluation of the membrane potential for threshold detection can be postponed. Importantly, for this approach to work, it is necessary for the estimated threshold crossing time to be before the actual threshold crossing. At the estimated threshold crossing time, then, detection of threshold crossing is a simple comparison of the exact solution with the threshold, by returning to the beginning of the event handling cycle (realized through issuing a self-event at the time of the Newton estimate calculation) and detecting whether threshold was reached.
The essential assumption of this mechanism is that every Newton iteration step is underestimating the threshold passage time, which allows it to approach the threshold crossing with several iterations without the risk of overestimating it. The discussion in the next section and the subsequent mathematical proof focus on showing that provided the excitatory synaptic currents decay faster than the inhibitory currents (τeν < τiμ for all combinations of ν and μ), the Newton iteration step is underestimating threshold passage time.
3. Threshold Passage Time Is Underestimated by Newton Iteration
The formal proof we present in this section is largely analogous to the proof in the original analysis (Carnevale & Hines, 2006; Hines & Carnevale, 2002), except that we base it on a generalization of the underlying Carnevale-Hines lemma, for which we give a proof here. This generalization allows us to lift the unphysiological constraints in the original analysis. The actual proof consists of two parts. The first part shows that when at t0 the derivative m′ = dm/dt ⩽ 0, the membrane potential will stay below threshold at least until a new event arrives. The second part shows that when at t0 the derivative m′ satisfies m′>0, the Newton iteration formula te = (1 − m)/m′ + t0 underestimates threshold passage time. Before presenting the formal proof based on the exact solution for m(t), we analyze the physics, which provides better insight into the actual underlying mechanisms.
3.1. When Will Newton Iteration Underestimate Threshold Passage Time? A Physical Analysis.
Let us first analyze the situation in which no inhibitory currents are present (e.g., the top black lines in Figure 1). We know that the excitatory currents are decaying; their contribution to m′ will reduce over time, and the only uncertainty is in the leak currents. If m′ ⩽ 0, the excitatory synaptic currents cannot overcome the leak current at the present membrane potential, and therefore after decaying further, they are definitely unable to do so. On the other hand, if m′>0, the leak current is growing and therefore is reducing m′, so m′ < m′0 until threshold is reached or until m′ reverses sign and the membrane potential moves away from threshold never to reach it. From this, it follows that either the threshold will never be reached and every finite estimate is underestimating threshold passage time, or m′ satisfies the inequality m′(t) ⩽ m′0 up to reaching the threshold, and we know that no threshold passage takes place before the estimated time te.
Now let us examine the situation where, in addition, double exponential inhibitory currents are present. Again, examples are given in Figure 1, where gray lines indicate cases where the inhibitory currents are growing and black lines cases where there are only decaying inhibitory currents. The first observation is that a growth of the inhibitory current leads to reduction of m′, so adding the growth component associated with jμ (see equation 2.1) will only strengthen the arguments we can obtain when assuming it to be equal to zero. So assuming jμ = 0 and m′0 > 0 and m0 ⩾ 0, we know that the individual synaptic currents are decaying toward 0.
If, however, the excitatory contribution were to decay slowly while the inhibition decayed rapidly, then the actual resulting synaptic current might be a growing depolarizing current, and we might overestimate firing time. It is this possibility of unmasking a more slowly decaying depolarizing current through the fast decay of an inhibitory current that leads to the only real restriction on the time constants in the model: the decay times for excitatory synaptic currents should be faster than those for inhibitory synaptic currents. If the inhibitory currents decay more slowly than the excitatory currents, we know that the total synaptic current decreases faster than expected on the basis of the excitatory decay time constants alone and might even become hyperpolarizing. Furthermore, while m′>0, the leak current is growing and therefore also reducing m′.
Taken together, we find that m′ < m′0 until threshold is reached or until m′ reverses sign and the membrane potential moves away from threshold, never to reach it. So under the assumption that inhibitory synaptic currents decay more slowly than excitatory synaptic currents, the condition m′ ⩽ m′0 is fulfilled, and we know that the estimated threshold passage time te based on Newton iteration underestimates the threshold passage time.
If, contrary to our earlier assumption, m0 < 0, then we have two separate cases. In the first case, the excitatory synaptic current is larger than the inhibitory synaptic current, in which situation all the arguments above apply. In the second case, we have a dynamics dominated by the leak current. For this leak-dominated phase, the argument above does not apply, but because the leak current acts toward the resting potential and not toward the threshold, it is clear that no threshold passage will take place and any estimate will underestimate threshold passage time. Importantly, no comparison was needed between the time constants of the synaptic currents and the membrane time constants, which shows that the value of the membrane time constant does not affect the validity of the argument.
From these arguments, we also see why it is difficult to include NMDA-like currents, which are best modeled by an excitatory double exponential current. Their growing excitatory contribution accelerates the rate at which the membrane potential approaches the threshold after activation of the synapse, causing the Newton iteration to overestimate the threshold passage time.
3.2. Generalized Carnevale-Hines Lemma.
The exact expression of the membrane potential is built from a large number of double exponential functions. Establishing upper (lower) bounds on such a sum of double exponentials is strongly simplified by the Carnevale-Hines lemma, which allows us to replace one double exponential with another that for all times is larger (smaller). The lemma follows from the following corollary, which gives us the monotonic development of double exponentials when viewed as a function of one of the decay times:
We start by showing that the function fλ can be extended to a continuous function of μ by adding the point μ = λ. The nominator and denominator used in the definition of fλ are 0 at μ = λ, and their derivatives with respect to μ exist; furthermore, the derivative of the denominator is equal to 1 and therefore nonzero. This means that the preconditions of l'Hôpital's rule are satisfied. From l'Hôpital's rule, we know that the limit exists, and we find limμ→λfλ(μ, t) = te−μt (i.e., the alpha function). After extending fλ with this limit, continuity in μ follows from observing that l'Hôpital's rule is based on the fact that the left and right limits are equal, and hence no discontinuity occurs at μ = λ in the extended fλ.
Equality follows from choosing t = 0 for which we have fλ(μ, 0) = 0 for all values of μ. For t ≠ 0, the inequality is a direct consequence of the corollary.
3.3. Movement Away from Threshold Implies Threshold Will Never Be Reached.
We start this analysis from the exact solution of the membrane potential m(t) given in equation 2.6. The purpose is to show that the exact solution has an upper bound given by the line me(t) = m0 + m′0(t − t0). Instead of showing that this is true because m′(t) ⩽ m′0 as in the physical analysis, we now show it on the basis of the exact solution itself.
3.4. Movement Toward Threshold Slows Down.
Let us summarize the behavior of the membrane potential in this model as analyzed by both the integrated-current based argument and the Carnevale-Hines lemma. When excitatory synaptic currents dominate (see equations 3.15 and 3.17), as in the early stages of the examples given in Figure 2, the membrane potential increases but always stays below the linear extrapolation. When inhibitory synaptic currents dominate (see equation 3.13), the membrane potential moves away from threshold, and if the currents are sufficiently strong, as in the lower curves shown in Figure 2, the polarity of the membrane will be reversed. When in the last part of these curves leak currents dominate (see equations 3.13 and 3.17), the membrane potential moves toward resting potential. If the latter takes place at negative membrane potentials, the membrane potential can exceed linear extrapolation, because decay of inhibitory synaptic current can quickly unmask an upward leak current. Although in this situation, the linear extrapolation is exceeded, no overestimation of threshold passage time will occur because no threshold passage takes place in the leak-dominated phase. In conclusion, we established that either no threshold crossing will take place or that threshold crossing takes place after the Newton estimate for threshold crossing time.
We have shown that the nonphysiological constraints on the time constants in the event-based integrate-and-fire Carnevale-Hines integration scheme can be relaxed. The only restriction is that the slowest decaying excitatory synaptic current (νs) should decay faster than the fastest decaying inhibitory current (μf), that is, (. This makes the model applicable to a much wider range of neuron and synapse types than previously thought. Only minor modifications have to be made in the initialization phase of the Carnevale-Hines integration scheme as implemented in NEURON, none of which affects the speed of the algorithm. These modifications will be incorporated in NEURON's standard IntFire4 mechanism (M. L. Hines to R. A. J. van Elburg, September 2008). The code generating the figures in this letter is available from ModelDB under accession number 115357.
Although we have analyzed a current-based model, the argument about threshold passage time, which is based on our physical analysis, can, with slight modifications, be applied to conductance-based integrate-and-fire models as well. Conductance-based models have synaptic reversal potentials, but these do not weaken our argument, as can be seen as follows. On approaching the reversal potential, the driving force for the excitatory current is reduced. The excitatory currents will therefore be reduced further than would be expected from channel kinetics alone. The driving force of the inhibitory currents is enhanced on approaching the threshold potential. The inhibitory currents will therefore be less reduced than would be expected from channel kinetics alone. Thus, including reversal potentials will only strengthen the argument about threshold passage time. The only crucial ingredient still missing, then, in order to apply the Carnevale-Hines integration scheme to conductance-based integrate-and-fire models is the lack of an exact solution for the membrane potential or another numerically efficient way to calculate the membrane potential after a time step of arbitrary size.
R.vE. was supported by the Computational Life Sciences Program of the Netherlands Organization for Scientific Research (NWO, CLS2003, 635.100.005) and by the Dutch Companion Project funded by the Dutch agency for innovation and sustainable development SenterNovem (SenterNovem, IS053013). We thank Markus Diesmann for comments on an earlier version of this letter.