Abstract

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.

1.  Introduction

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.

The model is event based and assumes that communication between neurons is completely dependent on action potentials generated in the axo-somatic region. As a result, the time of occurrence of the action potential contains all the information a single neuron communicates to the neurons it innervates. In the model used here, an event is generated when the neuron passes the firing threshold; axonal propagation and presynaptic delays are accounted for by delivering the event to the postsynaptic cell with an appropriate time delay, which is typically on the order of a few milliseconds. The postsynaptic part of the synapse responds by generating a synaptic current in the postsynaptic cell. Physiologically, the time course of synaptic currents is determined by the association rate of the neurotransmitter to the receptor, the dissociation rate of the neurotransmitter from the receptor, the removal rate of neurotransmitter from the synaptic cleft, and the driving synaptic reversal potential (Destexhe, Mainen, & Sejnowski, 1994). For current-based integrate-and-fire models, to which the Carnevale and Hines model belongs, it is assumed that the synaptic current follows receptor opening, while the induced change in driving force due to a changing membrane potential is ignored. Excitatory synapses and inhibitory synapses are included as membrane potential independent currents Eν and Iμ, respectively. The subscripts ν and μ indicate the different excitatory and inhibitory synaptic subtypes, which can each have their own kinetics. Because the currents follow the opening and closing of the receptors, their time course is completely specified by the average time course of the open states of the receptor-neurotransmitter complex. To describe these, we use the variables eν and Iμ for the excitatory and inhibitory synapses, respectively. On the arrival of an event at an excitatory synapse, the receptors in the excitatory synapse are instantaneously opened, reflecting fast transient AMPA binding dynamics. In the model, this is reflected by adding the weight of the synapse we ⩾ 0 to Eν. On the arrival of an event at an inhibitory synapse, however, we first get a fast increase in the amount of receptor-neurotransmitter complex in the closed state, followed by a transition of these receptor-neurotransmitter complexes to the open state. This change in amount of receptor-neurotransmitter complex is modeled by adding the weight of the particular connection wi ⩽ 0 to an auxiliary variable jμ describing the closed state of the receptor-neurotransmitter complex. The sign of the weight indicates the inhibitory nature of the resulting current. After this event handling, the variables eν, jμ, and iμ evolve according to the following differential equations, which link the slowest relevant timescales of the kinetic models underlying receptor opening and closing of the AMPA and GABA receptor to a time development model for the synaptic currents:
formula
2.1
with the time constants τeν for the excitatory decay time and τjμ, τiμ for the inhibitory rise and decay time. In addition, there is a parameter ajμ (where the appearance of τjμ−1 might have been anticipated), which acts as normalization constant and is chosen such that an event-induced change in jμ by an amount wi results in a maximal change of wi in iμ. For technical reasons, we will assume that τjμ < τiμ, which basically restricts the mechanism underlying the double exponential current to mechanisms in which the slow time constant acts on i and not on j. Although this can be at odds with the actual biophysics, it poses no real constraint on this phenomenological model, since its spike output is sensitive only to the shape of the inhibitory current, which is insensitive to an interchange of τjμ and τiμ provided the normalization constants are adjusted accordingly.
In the model, the membrane potential is represented by the variable m, and the resting potential and spike threshold are identified with m = 0 and m = 1, respectively. When the membrane potential deviates from the resting potential, then, in the absence of synaptic currents, it decays back to resting potential and will do so with the membrane time constant τm. Furthermore, the excitatory E = ∑νEν = ∑νaeνeν and inhibitory I = ∑μIμ = ∑μaiμiμ synaptic currents act on the membrane potential. The constants aeν and aiμ are normalization constants and are chosen such that an isolated instantaneous change in eν by an amount we > 0 induces a maximum depolarization of we and, similarly, such that an isolated instantaneous change in jν by an amount wi < 0 induces a maximum hyperpolarization of wi (Carnevale & Hines, 2006; Hines & Carnevale, 2004). Alternatively, they can be chosen to normalize the charge transfer (Hansel et al., 1998). Putting everything together, the differential equation describing the membrane potential becomes
formula
2.2
This differential equation, together with the accompanying differential equations for the synapses, has an exact solution for the subthreshold dynamics. These differential equations are linear in the currents, and therefore the solution is completely analogous to the case with a single excitatory and a single inhibitory current (Carnevale & Hines, 2006; Hines & Carnevale, 2004). Using the reciprocals of the time constants kx = 1/τx, we can write the exact solution as follows:
formula
2.3
formula
2.4
formula
2.5
formula
2.6
In this expression, t0 refers to the last time preceding t at which m, eν, jμ, and iμ were evaluated; the values of those variables at t0 are denoted by m0, eν,0, jμ,0, iμ,0. The constants aeν, ajμ, and aiμ used in the differential equations are absorbed into new constants b together with part of the kx, that is, . The challenge now is to derive from the subthreshold behavior the point in time where the membrane potential crosses spiking threshold, the topic of the next section.

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.

To answer the question when Newton iteration underestimates threshold passage time, we need to look at the different components contributing to the membrane potential derivative m′ = dm/dt. We have three kinds of currents: a leak current, which acts in the direction of the resting membrane potential; excitatory synaptic currents, which depolarize the membrane; and inhibitory currents, which hyperpolarize the membrane. If the Newton iteration estimate of threshold passage time te,
formula
3.1
is required to underestimate threshold passage time, then that puts limitations on the possible time course of the currents contributing to m′. If the membrane potential does not cross threshold θ between t0 and te, then it satisfies the inequality
formula
3.2
Using simple arguments, we can derive conditions on the time constants for which the sum of synaptic and leak currents is decreasing. From the inequality above, we can immediately see that if during t0 < t < te the total current shows no growth (i.e., m′(t) ⩽ m0), then no threshold passage will take place in this time interval. Therefore, under the stronger assumption that the total current is decreasing the Newton iteration will underestimate threshold passage time.

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′ < m0 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) ⩽ m0 up to reaching the threshold, and we know that no threshold passage takes place before the estimated time te.

Figure 1:

Illustration of the physical analysis. The figures show the derivative of the membrane potential, m′, against time. (a) The entire time course is shown. (b) Only the initial part is shown. If m0, m0 > 0, then m′ will stay below m0 up to the Newton estimated threshold passage time indicated by the vertical dark gray bar. The top black line represents the case where initially all synaptic currents are excitatory and the inhibitory synaptic current is zero for all times; it is shown up to the point where threshold is reached. The lower-lying black lines are at larger initial inhibitory (and hence excitatory) synaptic currents (mi = aii = −0.20, − 0.15, − 0.10, − 0.05) but without the growth component (j = 0). The accompanying gray lines show the effect of adding the growth component (j = −0.025, − 0.050, − 0.075, − 0.100). Parameters used: τm = 30, τe = 3, τi = 9, τj = 2, m0 = 0.5, m0 = 0.25. Simulations performed in NEURON.

Figure 1:

Illustration of the physical analysis. The figures show the derivative of the membrane potential, m′, against time. (a) The entire time course is shown. (b) Only the initial part is shown. If m0, m0 > 0, then m′ will stay below m0 up to the Newton estimated threshold passage time indicated by the vertical dark gray bar. The top black line represents the case where initially all synaptic currents are excitatory and the inhibitory synaptic current is zero for all times; it is shown up to the point where threshold is reached. The lower-lying black lines are at larger initial inhibitory (and hence excitatory) synaptic currents (mi = aii = −0.20, − 0.15, − 0.10, − 0.05) but without the growth component (j = 0). The accompanying gray lines show the effect of adding the growth component (j = −0.025, − 0.050, − 0.075, − 0.100). Parameters used: τm = 30, τe = 3, τi = 9, τj = 2, m0 = 0.5, m0 = 0.25. Simulations performed in NEURON.

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 m0 > 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′ < m0 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′ ⩽ m0 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:

Corollary 1.
For t, μ, λ ∈ ℜ the functions
formula
3.3
are defined for μ ≠ λ and by including the limits limμ→λfλ(μ, t) into fλ can be extended to a function continuous in μ. For fixed t, the extended function fλ is a monotonically decreasing continuous function of μ.

Proof.

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

The next step is to prove that the derivative of the function fλ with respect to μ is negative almost everywhere for t ≠ 0. This derivative is given by
formula
3.4
On the right-hand side of this equation, the first factor is clearly positive, and therefore we need to prove that the second factor is negative almost everywhere to show that the function is monotonically decreasing. To show that the second factor is negative almost everywhere, we prove that it has a nonpositive maximum at λ = μ. If we examine the derivative of the second factor,
formula
3.5
we find that it is 0 for μ = λ, positive for μ < λ, and negative for μ>λ, showing that there is indeed a maximum (t(μ − λ) + 1 − e(μ−λ)t) = 0 in the second factor at λ = μ. As a result, we find that dfλ(μ, t)/dμ < 0 for μ ≠ λ, and the only step left is to prove that dfλ(μ, t)/dμ is continuous.
The expression for dfλ(μ, t)/dμ is indeterminate at μ = λ, but the preconditions for l'Hôpital's rule are satisfied, indicating that this derivative is continuous at this point. The values of the dfλ(μ, t)/dμ at μ = λ are most easily calculated by inserting the Taylor expansion for t(μ − λ) + 1 − e(μ−λ)t into equation 3.4:
formula
3.6
The indeterminacy at μ = λ is now canceled, and we find
formula
3.7
This expression is negative, and thus the derivative dfλ(μ, t)/dμ is negative everywhere.

Lemma 2.
(generalized Carnevale-Hines lemma). If μ2 > μ1 and μ1, μ2 ≠ λ, then
formula
3.8

Proof.

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 + m0(tt0). Instead of showing that this is true because m′(t) ⩽ m0 as in the physical analysis, we now show it on the basis of the exact solution itself.

As before, the terms describing growth of the inhibition, (i.e., those terms related to jμ), can be discarded, but now we use the Carnevale-Hines lemma to achieve this. From the lemma, we find that for τjμ < τiμ, the factors,
formula
3.9
found in the exact solution for m(t) (see equation 2.6) are all positive. The factors multiplying these expression, biμ and bjμ, are also positive, but the factors jμ,0 are negative. The terms therefore in which these appear are negative, and we obtain an upper bound for m by dropping these terms:
formula
3.10
In the next step, we use the assumption that τeν ≤ τiμ for all combinations of μ and ν to be able to apply the Carnevale-Hines lemma. Using an arbitrary excitatory decay time τeν′ we replace the dependence on τiμ in the iμ,0-related terms with a τeν′ dependence. If we fix τeν′ by choosing the largest excitatory decay time for it, we can use the lemma a second time and replace the τeν dependence in the eν,0-related terms with a τeν′ dependence as well:
formula
3.11
Now we can use m′ = −kmm + ∑νaeνeν + ∑μaiμiμ to obtain
formula
3.12
We assumed that m′ ⩽ 0, and because the other factors multiplying it are all positive, we can remove the associated term from the inequality so that we obtain
formula
3.13
This expression is now fully equivalent to the one found in the original treatment (Carnevale & Hines, 2006; Hines & Carnevale, 2004), where it is shown that the derivative with respect to time of the term multiplying m0 is negative, and because it is 1 at t = t0, it will be smaller than 1 at larger times while it is also positive. Because m0 < 1, this shows that m(t) stays between −∞ and 1 and hence below threshold.

3.4.  Movement Toward Threshold Slows Down.

If m′>0 the step from equation 3.12 to equation 3.13 is not allowed, but using our corollary, we can replace the last term in equation 3.12 with the alpha function belonging to the slowest decay time out of τm and τe, ν′:
formula
3.14
By the same argument as used before, the first two terms combined are smaller than or equal to m0; also, is smaller than 1, so taking this together with m0, (tt0) ⩾ 0, we obtain
formula
3.15
When m0 < 0, the second term on the right-hand side of equation 3.14 is negative, and we can drop it from the inequality, which is the correct thing to do if τe, ν′ < τm. If τe, ν′ > τm, we can reorganize the first two terms in such a way that we can replace it with eke, ν′(tt0):
formula
3.16
from which we see that we have a free choice whether to keep the factor or ekm(tt0) from the first two terms. It will be convenient to make a choice that fits the last term:
formula
3.17
This expression shows that in the area where the linear extrapolation is below zero, we have no guarantee that m(t) is under the linear extrapolation, but the actual upper bound is below zero, and no threshold passage occurs. When the linear extrapolation crossed zero, we again find m(t) ⩽ m0 + m′(tt0), and we see that Newton iteration is underestimating the threshold crossing time. This seemingly strange behavior of the upper bound does not indicate a magical zero crossing of m(t) at the point where m0 + m′(tt0) changes sign, but reflects our lack of knowledge about the system. If the system was dominated by excitatory currents during the whole period, then we know from our physical analysis that the linear extrapolation was above the real curve all the time. If, however, in this hyperpolarized situation, the system was dominated by a leak current that becomes unmasked due to fast decay of inhibitory currents, then the real curve went above the linear extrapolation. But because the whole movement is toward resting potential and will never cross the resting potential, the linear extrapolation has to cross the membrane potential before passing through zero.

3.5.  Summary.

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.

Figure 2:

Illustration of the argument based on the Carnevale-Hines lemma. The figures show the membrane potential, m, against time. (a) The whole time course is shown. (b) Only the initial part of the time course is shown. If m0, m0 > 0, then m will stay below m0 + m0t up to the estimated threshold passage time indicated by the point where the tangent at t0 (dark gray line) crosses threshold θ (dashed dark gray line). The top black line represents the case where initially all synaptic currents are excitatory and the inhibitory synaptic current is zero for all times. It is shown up to the point where threshold is reached. The lower-lying black lines are at larger initial inhibitory synaptic currents (mi = aii = −0.20, − 0.15, − 0.10, − 0.05) but without the growth component (j = 0). The accompanying gray lines show the effect of adding the growth component (j = −0.025, − 0.050, − 0.075, − 0.100). Parameters used: τm = 30, τe = 3, τi = 9, τj = 2, m0 = 0.5, m0 = 0.25. Simulations performed in NEURON.

Figure 2:

Illustration of the argument based on the Carnevale-Hines lemma. The figures show the membrane potential, m, against time. (a) The whole time course is shown. (b) Only the initial part of the time course is shown. If m0, m0 > 0, then m will stay below m0 + m0t up to the estimated threshold passage time indicated by the point where the tangent at t0 (dark gray line) crosses threshold θ (dashed dark gray line). The top black line represents the case where initially all synaptic currents are excitatory and the inhibitory synaptic current is zero for all times. It is shown up to the point where threshold is reached. The lower-lying black lines are at larger initial inhibitory synaptic currents (mi = aii = −0.20, − 0.15, − 0.10, − 0.05) but without the growth component (j = 0). The accompanying gray lines show the effect of adding the growth component (j = −0.025, − 0.050, − 0.075, − 0.100). Parameters used: τm = 30, τe = 3, τi = 9, τj = 2, m0 = 0.5, m0 = 0.25. Simulations performed in NEURON.

4.  Discussion

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.

Within the Carnevale-Hines integration scheme based on Newton iteration, slowly growing excitatory currents (like NMDA currents) cannot be included. We can, however, under some conditions, replace the Newton iteration by a higher-order extrapolation scheme. Let us take, for example, a synaptic model for NMDA currents using the same differential equations as those used for the inhibitory current (admittedly ignoring the magnesium block, which is outside the scope of our treatment here). We now take i, j to represent the double exponential excitatory currents and auxiliary variable, respectively, and leave out the inhibitory currents. To do this, we only need to change the sign of j, so we assume j ⩾ 0. If we now further assume that τi > τm, then the exact solution has an upper bound:
formula
4.1
This quadratic upper bound for the exact solution was found by using appropriate α-functions as upper bounds for the double exponentials, that is, by repetitively applying the generalized Carnevale-Hines lemma. We again obtain an estimate for threshold crossing if we determine where this upper bound crosses threshold. This new estimate does not yet solve the problem of unmasking the slowly decaying part of the existing excitatory current by rapid decay of the inhibitory current, but we expect that unmasking can be dealt with in a similar way, leading to extra quadratic terms in the upper bound. The reason for this is that only existing decaying currents can be unmasked; those currents lead to a double exponential contribution to the development of the membrane potential, and putting an upper bound on a combination of two double exponentials involves two times in succession the replacement of a double exponential by an α-function leading to a (tt0)2 term. If this reasoning is correct, then in even more general physiological conditions——an estimate for threshold crossing can be found by taking the positive root of a quadratic polynomial.

Acknowledgments

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.

References

Abbott
,
L.
(
1999
).
Lapicque's introduction of the integrate-and-fire model neuron (1907)
.
Brain Res. Bull.
,
50
(
5/6
),
303
304
.
Bosman
,
L. W. J.
,
Heinen
,
K.
,
Spijker
,
S.
, &
Brussaard
,
A. B.
(
2005
).
Mice lacking the major adult GABA(A) receptor subtype have normal number of synapses, but retain juvenile IPSC kinetics until adulthood
.
J. Neurophysiology
,
94
(
1
),
338
346
.
Bosman
,
L. W. J.
,
Rosahl
,
T. W.
, &
Brussaard
,
A. B.
(
2002
).
Neonatal development of the rat visual cortex: Synaptic function of GABA(A) receptor alpha subunits
.
J. Physiology-London
,
545
(
1
),
169
181
.
Brette
,
R.
(
2006
).
Exact simulation of integrate-and-fire models with synaptic conductances
.
Neural Comput.
,
18
(
8
),
2004
2027
.
Brette
,
R.
,
Rudolph
,
M.
,
Carnevale
,
T.
,
Hines
,
M.
,
Beeman
,
D.
, &
Bower
,
J.
, et al
(
2007
).
Simulation of networks of spiking neurons: A review of tools and strategies
.
J. Comput. Neurosci.
,
23
(
3
),
349
398
.
Carnevale
,
N. T.
,
Hines
,
M. L.
(
2006
).
The NEURON Book
,
Cambridge
:
Cambridge University Press
.
Destexhe
,
A.
,
Mainen
,
Z. F.
, &
Sejnowski
,
T. J.
(
1994
).
Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism
.
J. Comput. Neurosci.
,
1
,
195
230
.
Hansel
,
D.
,
Mato
,
G.
,
Meunier
,
C.
, &
Neltner
,
L.
(
1998
).
On numerical simulations of integrate-and-fire neural networks
.
Neural Comput.
,
10
(
2
),
467
483
.
Hausser
,
M.
, &
Roth
,
A.
(
1997
).
Estimating the time course of the excitatory synaptic conductance in neocortical pyramidal cells using a novel voltage jump method
.
J. Neurosci.
,
17
(
20
),
7606
7625
.
Hines
,
M. L.
, &
Carnevale
,
N. T.
(
2002
).
Efficient discrete event simulation of spiking neurons in neuron
.
Poster for the Society for Neuroscience Meeting. Available online at http://www.neuron.yale.edu/neuron/papers/discrete/discrete_event_poster.pdf
.
Hines
,
M. L.
, &
Carnevale
,
N. T.
(
2004
).
Discrete event simulation in the NEURON environment
.
Neurocomputing
,
58–60
,
1117
1122
.
Karube
,
F.
,
Kubota
,
Y.
, &
Kawaguchi
,
Y.
(
2004
).
Axon branching and synaptic bouton phenotypes in GABAergic nonpyramidal cell subtypes
.
J. Neurosci.
,
24
(
12
),
2853
2865
.
Kleppe
,
I. C.
, &
Robinson
,
H. P. C.
(
1999
).
Determining the activation time course of synaptic AMPA receptors from openings of colocalized NMDA receptors
.
Biophys. J.
,
77
(
3
),
1418
1427
.
Lapicque
,
L.
(
1907
).
Recherches quantitatives sur l'excitation électrique des nerfs traitée comme une polarization
.
J. Physiol. Pathol. Gen.
,
9
,
620
635
.
Morrison
,
A.
,
Straube
,
S.
,
Plesser
,
H. E.
, &
Diesmann
,
M.
(
2007
).
Exact subthreshold integration with continuous spike times in discrete-time neural network simulations
.
Neural Comput.
,
19
(
1
),
47
79
.
Rudolph
,
M.
, &
Destexhe
,
A.
(
2006a
).
Analytical integrate-and-fire neuron models with conductance-based dynamics for event-driven simulation strategies
.
Neural Comput.
,
18
(
9
),
2146
2210
.
Rudolph
,
M.
, &
Destexhe
,
A.
(
2006b
).
Event-based simulation strategy for conductance-based synaptic interactions and plasticity
.
Neurocomputing
,
69
(
10–12
),
1130
1133
.
Tateno
,
T.
, &
Robinson
,
H. P. C.
(
2006
).
Rate coding and spike-time variability in cortical neurons with two types of threshold dynamics
.
J. Neurophysiology
,
95
(
4
),
2650
2663
.
Trevelyan
,
A. J.
, &
Jack
,
J. J. B.
(
2002
).
Detailed passive cable models of layer 2/3 pyramidal cells in rat visual cortex at different temperatures
.
J. Physiol.
,
539.2
,
623
636
.
Tsodyks
,
M.
,
Mitkov
,
I.
, &
Sompolinsky
,
H.
(
1993
).
Pattern of synchrony in inhomogeneous networks of oscillators with pulse interactions
.
Physical Review Letters
,
71
(
8
),
1280
1283
.
van Elburg
,
R. A. J.
,
Bosman
,
L. W.
,
Ridder
,
M. C.
,
Brussaard
,
A. B.
, &
van Ooyen
,
A.
(
2009
).
GABAA receptor plasticity provides homeostasis of neuronal activity in an integrate-and-fire model for neocortical microcircuitry
.
Unpublished manuscript
.
Zhang
,
Z. W.
(
2004
).
Maturation of layer V pyramidal neurons in the rat prefrontal cortex: Intrinsic properties and synaptic function
.
J. Neurophysiology
,
91
(
3
),
1171
1182
.