5 Nonlinear Integrate-and-Fire Models

5.2 Exponential Integrate-and-Fire Model

In the exponential integrate-and-fire model (156), the differential equation for the membrane potential is given by

τddtu=-(u-urest)+ΔTexp(u-ϑrhΔT)+RI;\tau{{\text{d}}\over{\text{d}}t}u=-(u-u_{\rm rest})+\Delta_{T}\,\exp\left({u-% \vartheta_{rh}\over\Delta_{T}}\right)+R\,I\,; (5.6)

The first term on the right-hand-side of Eq. (5.6) is identical to Eq. (5.3) and describes the leak of a passive membrane. The second term is an exponential nonlinearity with ‘sharpness’ parameter ΔT\Delta_{T} and ‘threshold’ ϑrh\vartheta_{rh}.

The moment when the membrane potential reaches the numerical threshold θreset\theta_{\rm reset} defines the firing time t(f)t^{(f)}. After firing, the membrane potential is reset to uru_{r} and integration restarts at time t(f)+Δabst^{(f)}+\Delta^{\rm abs} where Δabs\Delta^{\rm abs} is an absolute refractory time, typically chosen in the range 0<Δabs<5ms0<\Delta^{\rm abs}<5ms. If the numerical threshold is chosen sufficiently high, θresetϑ+ΔT\theta_{\rm reset}\gg\vartheta+\Delta_{T}, its exact value does not play any role. The reason is that the upswing of the action potential for uϑ+ΔTu\gg\vartheta+\Delta_{T} is so rapid, that it goes to infinity in an incredibly short time (515). The threshold θreset\theta_{\rm reset} is introduced mainly for numerical convenience. For a formal mathematical analysis of the model, the threshold can be pushed to infinity.

Example: Rheobase threshold and interpretation of parameters

The exponential integrate-and-fire model is a special case of the general nonlinear model defined in Eq. (5.2) with a function

f(u)=-(u-urest)+ΔTexp(u-ϑrhΔT).f(u)=-(u-u_{\rm rest})+\Delta_{T}\,\exp\left({u-\vartheta_{rh}\over\Delta_{T}}% \right)\,. (5.7)

In the absence of external input (I=0I=0), the differential equation of the exponential integrate-and-fire model (5.6) has two fixed points, defined by the zero-crossings f(u)=0f(u)=0; cf. Fig. 5.1A. We suppose that parameters are chosen such that ϑrhurest+ΔT\vartheta_{rh}\gg u_{\rm rest}+\Delta_{T}. Then the stable fixed point is at uurestu\approx u_{\rm rest} because the exponential term becomes negligibly small for uϑrh-ΔTu\ll\vartheta_{rh}-\Delta_{T}. The unstable fixed point which acts as a threshold for pulse input lies to the right-hand side of ϑrh\vartheta_{rh}.

If the external input increases slowly in a quasi-constant fashion, the two fixed points move closer together until they finally merge at the bifurcation point; cf. Fig. 5.1B. The voltage at the bifurcation point can be determined from the condition df/du=0df/du=0 to lie at u=ϑrhu=\vartheta_{rh}. Thus ϑrh\vartheta_{rh} is the threshold found with constant (rheobase) current, which justifies its name.

Fig. 5.3: Exponential and leaky integrate-and-fire model. The function f(u)f(u) is plotted for different choices of the ’sharpness’ of the threshold (ΔT=1, 0.5, 0.25, 0.05\Delta_{T}=1,\,0.5,\,0.25,\,0.05 mV) In the limit ΔT0\Delta_{T}\to 0 the exponential integrate-and-fire model becomes equivalent to a leaky integrate-and-fire model (dashed line). The inset shows a zoom onto the threshold region (dotted box).

Example: Relation to the Leaky Integrate-and-Fire Model

In the exponential integrate-and-fire model, the voltage threshold ϑ\vartheta for pulse input is different from the rheobase threshold ϑrh\vartheta_{rh} for constant input (Fig. 5.1). However, in the limit ΔT0\Delta_{T}\to 0, the sharpness of the exponential term increases and ϑ\vartheta approaches ϑrh\vartheta_{rh} (Fig. 5.3). In the limit, ΔT0\Delta_{T}\to 0, we can approximate the nonlinear function by the linear term

f(u)=-(u-urest)foru<ϑrhf(u)=-(u-u_{\rm rest})\quad{\rm for}\quad u<\vartheta_{rh} (5.8)

and the model fires whenever uu reaches ϑrh=ϑ\vartheta_{rh}=\vartheta. Thus, in the limit ΔT0\Delta_{T}\to 0, we return to the leaky integrate-and-fire model.

5.2.1 Extracting the Nonlinearity from Data

Why should we choose an exponential nonlinearity rather than any other nonlinear dependence in the function f(u)f(u) of the general nonlinear integrate-and-fire model? Can we use experimental data to determine the ‘correct’ shape of f(u)f(u) in Eq. (5.2)?

We can rewrite the differential equation (5.2) of the nonlinear integrate-and-fire model by moving the function f(u)f(u) to the left-hand side and all other terms to the right-hand-side of the equation. After rescaling with the time constant τ\tau, the nonlinearity f~(u)=f(u)/τ\tilde{f}(u)=f(u)/\tau is

f~(u(t))=1CI(t)-ddtu(t);\tilde{f}(u(t))={1\over C}\,I(t)-{{\text{d}}\over{\text{d}}t}u(t)\,; (5.9)

where C=τ/RC=\tau/R can be interpreted as the capacity of the membrane.

In order to determine the function f~(u)\tilde{f}(u), an experimentalist injects a time-dependent current I(t)I(t) into the soma of a neuron while measuring with a second electrode the voltage u(t)u(t). From the voltage time course, one finds the voltage derivative du/dtdu/dt.

Fig. 5.4: Extracting nonlinear integrate-and-fire models from data. The function f(u)f(u) characterizing the nonlinearity of an integrate-and-fire model according to Eq. (5.2) is derived from experimental data using random current injection into neurons. A Cortical pyramidal cells. Experimental data points (symbols) and fit by an exponential integrate-and-fire model. B As in A, but for an inhibitory interneuron. Data courtesy of Laurent Badel and Sandrine Lefort (35).

A measurement at time tt yields a value u(t)u(t) (which we use as value along the xx-axis of a plot) and a value [(I(t)/C)-(du/dt)][(I(t)/C)-(du/dt)] (which we plot along the yy axis). With a thousand or more time points per second, the plot fills up rapidly. For each voltage uu there are many data points with different values along the y-axis. The best choice of the parameter CC is the one that minimizes the width of this distribution. At the end, we average across all points at a given voltage uu to find the empirical function (35)

f~(u(t))=1CI(t)-ddtu(t)\tilde{f}(u(t))=\left\langle{1\over C}\,I(t)-{{\text{d}}\over{\text{d}}t}u(t)\right\rangle (5.10)

where the angular brackets indicate averaging. This function is plotted in Fig. 5.4. We find that the empirical function extracted from experiments is well approximated by a combination of a linear and exponential term

f~(u)=-u-urestτ+ΔTτexp(u-ϑrhΔT)\tilde{f}(u)=-{{u-u_{\rm rest}}\over\tau}+{\Delta_{T}\over\tau}\,\exp\left({u-% \vartheta_{rh}\over\Delta_{T}}\right) (5.11)

which provides an empirical justification of the choice of nonlinearity in the exponential integrate-and-fire model.

We note that the slope of the curve at the resting potential is related to the membrane time constant τ\tau while the threshold parameter ϑrh\vartheta_{rh} is the voltage at which the function f~\tilde{f} goes through its minimum.

Example: Refractory Exponential Integrate-and-Fire Model

The above procedure for determining the nonlinearity can be repeated for a set of data points restricted to a few milliseconds after an action potential (Fig. 5.6). After a spike, the threshold ϑrh\vartheta_{rh} is slightly higher which is one of the signs of refractoriness. Moreover, the location of the zero crossing urestu_{\rm rest} and the slope of the function f~\tilde{f} at urestu_{\rm rest} are different, which is to be expected since after a spike the sodium channel is inactivated while several other ion channels are open. All parameters return to the ‘normal’ values within a few tens of milliseconds. An exponential integrate-and-fire model where the parameters depend on the time since the last spike has been called the ‘refractory exponential integrate-and-fire model’ (35). The refractory exponential integrate-and-fire model predicts the voltage time course of a real neuron for novel time-dependent stimuli to a high degree of accuracy, if the input statistics is similar to the one used for parameter extraction (Fig. 5.5).

Fig. 5.5: Predicting the membrane voltage with an exponential integrate-and-fire model. A. Comparison of membrane voltage in experiments (thick line) with the predictions of the exponential integrate-and-fire model (thin line). The fit is excellent, except during a short period after a spike. B. Same as in A, but in a model with refractoriness. Modified from Badel et al. (34).
Fig. 5.6: Refractory effects in the exponential integrate-and-fire model. Top: Because of refractoriness immediately after a spike, the exponential integrate-and-fire model has a higher firing threshold and increased slope in the linear section. Data points and fit as in Fig. 5.4, but data points restricted to intervals 5-10ms (far left), 10-20ms (left), 20-30ms (right), or 30-50ms (far right) after a spike. As the time since the last spike increases, refractoriness decays and the parameters of the exponential integrate-and-fire model approach their standard values (dashed lines). Bottom: Sample voltage traces during and after a spike. Taken from Badel et al. (34).

5.2.2 From Hodgkin-Huxley to Exponential Integrate-and-Fire

In Section 4.2 of Chapter 4 we have already seen that the four-dimensional system of equations of Hodgkin and Huxley can be reduced to two equations. Here we show how to take a further step so as to arrive at a single nonlinear differential equation combined with a reset (246).

After appropriate rescaling of all variables, the system of two equations that summarizes a Hodgkin-Huxley model reduced to two dimensions can be written as

dudt\displaystyle{{\text{d}}u\over{\text{d}}t} =\displaystyle= F(u,w)+I\displaystyle F(u,w)+I (5.12)
dwdt\displaystyle{{\text{d}}w\over{\text{d}}t} =\displaystyle= ϵG(u,w)\displaystyle\epsilon\,G(u,w) (5.13)

which is just a copy of Eq. (4.34) in Section 4.6; note that time is measured in units of the membrane time constant τm\tau_{m} and that the resistance has been absorbed into the definition of the currrent II. The function F(u,w)F(u,w) is given by Eq. (4.3). The exact shape of the function G(u,w)G(u,w) has been derived in Section 4.2, but plays no role in the following. We recall that the fixed points are defined by the condition du/dt=dw/dt=0du/dt=dw/dt=0. For the case without stimulation I=0I=0, we denote the variables at the stable fixed point as urestu_{\rm rest} (resting potential) and wrestw_{\rm rest} (resting value of the second variable).

In the following we assume that there is a separation of time scales (ϵ1\epsilon\ll 1) so that the evolution of the variable ww is much slower than that of the voltage. As discussed in Section 4.6, this implies that all flow arrows in the two-dimensional phase plane are horizontal except those in the neighborhood of the uu-nullcline. In particular, after a stimulation with short current pulses, the trajectories move horizontally back to the resting state (no spike elicited) or horizontally leftward (upswing of an action potential) until they hit one of the branches of the uu-nullcline; cf. Fig. 4.22. In other words, the second variable stays at its resting value w=wrestw=w_{\rm rest} and can therefore be eliminated - unless we want to describe the exact shape of the action potential. As long as we are only interested in the initiation phase of the action potential we can assume a fixed value w=wrestw=w_{\rm rest}.

For constant ww, Eq. 5.12 becomes

dudt=F(u,wrest)+I=f(u)+I{{\text{d}}u\over{\text{d}}t}=F(u,w_{\rm rest})+I=f(u)+I (5.14)

which has the form of a nonlinear integrate-and-fire neuron. The resulting function f(u)f(u) is plotted in Fig. 5.7A. It has three zero-crossings: the first one (left) at urestu_{\rm rest}, corresponding to a stable fixed point; a second one (middle) which acts as a threshold ϑ\vartheta; and a third one to the right, which is again a stable fixed point and limits the upswing of the action potential. The value of the reset threshold θreset>ϑ\theta_{\rm reset}>\vartheta should be reached during the upswing of the spike and must be therefore be chosen between the second and third fixed point. While in the two-dimensional model the variable ww is necessary to describe the downswing of the action potential on a smooth trajectory back to rest, we replace the downswing in the nonlinear integrate-and-fire model by an artificial reset of the voltage variable to a value uru_{r} whenever uu hits θreset\theta_{\rm reset}.

If we focus on the region u<θresetu<\theta_{\rm reset}, the function f(u)=F(u,wrest)f(u)=F(u,w_{\rm rest}) is very well approximated by the nonlinearity of the exponential integrate-and-fire model (Fig. 5.7B).

Example: Exponential Activation of Sodium Channels

In the previous section, we followed a series of formal mathematical steps, from the two-dimensional version of the Hodgkin-Huxley model to a one-dimensional differential equation which looked like a combination of linear and exponential terms, i.e., an exponential integrate-and-fire model. For a more biophysical derivation and interpretation of the exponential integrate-and-fire model, it is, however, illustrative to start directly with the voltage equation of the Hodgkin-Huxley model, (2.4) - (2.5), and replace the variables hh and nn by their values at rest, hresth_{\rm rest} and nrestn_{\rm rest}, respectively. Furthermore, we assume that mm approaches instantaneously its equilibrium value m0(u)m_{0}(u). This yields

C   d   u   d   t=-gNa[m0(u)]3hrest(u-ENa)-gK(nrest)4(u-EK)-gL(u-EL)+I,C{{\text{d}}u\over{\text{d}}t}=-g_{\rm Na}[m_{0}(u)]^{3}\,h_{\rm rest}\,(u-E_{% \rm Na})-g_{\rm K}\,\left(n_{\rm rest}\right)^{4}\,(u-E_{\rm K})-g_{L}\,(u-E_{% L})+I\,, (5.15)

Potassium and leak currents can now be summed up to a new effective leak term geff(u-Eeffg^{\rm eff}\,(u-E^{\rm eff}. In the voltage range close to the resting potential the driving force (u-ENa)(u-E_{\rm Na}) of the sodium current can be well approximated by (urest-ENa)(u_{\rm rest}-E_{\rm Na}). Then the only remaining nonlinearity on the right-hand-side of Eq. (5.15) arises from the m0(u)m_{0}(u). For voltages around rest, m0(u)m_{0}(u) has, however, an exponential shape. In summary, the right-hand side of Eq. (5.15) can be approximated by a linear and an exponential term - and this gives rise to the exponential integrate-and-fire model (156).

Fig. 5.7: Approximating Hodgkin-Huxley by an exponential integrate-and-fire model. A. The value of the F(u,wrest)F(u,w_{\rm rest}) for fixed value of the second variable w=wrestw=w_{\rm rest} is plotted as a function of the voltage variable uu. The choice of a reset threshold θreset\theta_{\rm reset} is indicated. B. Solid line as in A, restricted to u<θresetu<\theta_{\rm reset}. The dashed line shows the approximation by an exponential integrate-and-fire model.