Abstract

The quadratic adaptive integrate-and-fire model (Izhikevich, 2003, 2007) is able to reproduce various firing patterns of cortical neurons and is widely used in large-scale simulations of neural networks. This model describes the dynamics of the membrane potential by a differential equation that is quadratic in the voltage, coupled to a second equation for adaptation. Integration is stopped during the rise phase of a spike at a voltage cutoff value Vc or when it blows up. Subsequently the membrane potential is reset, and the adaptation variable is increased by a fixed amount. We show in this note that in the absence of a cutoff value, not only the voltage but also the adaptation variable diverges in finite time during spike generation in the quadratic model. The divergence of the adaptation variable makes the system very sensitive to the cutoff: changing Vc can dramatically alter the spike patterns. Furthermore, from a computational viewpoint, the divergence of the adaptation variable implies that the time steps for numerical simulation need to be small and adaptive. However, divergence of the adaptation variable does not occur for the quartic model (Touboul, 2008) and the adaptive exponential integrate-and-fire model (Brette & Gerstner, 2005). Hence, these models are robust to changes in the cutoff value.

1.  Introduction

The question of finding a computationally simple and biologically realistic model of a neuron has been widely studied. The key problem is to find a neuron model that represents a compromise between simulation efficiency and flexibility to reproduce different firing patterns (Koch & Segev, 1998; Izhikevich, 2004; Rinzel & Ermentrout, 1989). Among the variety of computational neuron models, nonlinear spiking models with adaptation have recently been studied by several authors (Izhikevich, 2004; Brette & Gerstner, 2005; Touboul, 2008) and seem to stand out. They are relatively simple—mathematically tractable, efficiently implemented, and able to reproduce a large number of electrophysiological signatures. These models satisfy the equations
formula
1.1
where a and b are nonnegative parameters and F(v) is a regular strictly convex function satisfying the following assumption:

Assumption.

There exist ε > 0 and α > 0 for which F(v) ⩾ αv1+ε when v → ∞ (we will say that F grows faster than v1+ε when v → ∞).

A spike is emitted at time t* when the membrane potential v reaches a cutoff value Vc or when it blows up. At this time, the membrane potential is reset to a constant value c and the adaptation variable is updated to w(t*) + d, where w(t*) is the value of the adaptation variable at the time of the spike and d > 0 is the spike-triggered adaptation parameter.

In the absence of a cutoff value, the membrane potential blows up in finite time (see a proof of this fact in section 2), which allows replacement of the strict voltage threshold of classical (linear) integrate-and-fire neurons (Lapicque, 1907; Stein, 1967) by a more realistic smooth spike initiation (see Latham, Richmond, Nelson, & Nirenberg, 2000; Fourcaud-Trocme, Hansel, van Vreeswijk, & Brunel, 2003). Among these models, the quadratic adaptive model (Izhikevich, 2004) corresponds to the case where F(v) = v2 and has been recently used by Eugene Izhikevich and coworkers (Izhikevich & Edelman, 2008) in very large-scale simulations of neural networks. The adaptive exponential model (Brette & Gerstner, 2005) corresponds to the case where F(v) = evv. It is of interest because its parameters can be related to electrophysiological quantities and has been successfully fit to intracellular recordings of pyramidal cells (Clopath, Jolivet, Rauch, Lüscher, & Gerstner, 2007; Jolivet et al., 2008; Badel et al., 2008). The quartic model (Touboul, 2008) corresponds to the case where F(v) = v4 + 2av. It has the advantage of being able to reproduce all the behaviors featured by the other two and also self-sustained subthreshold oscillations, which are of particular interest in modeling certain nerve cells.

In these models, the reset mechanism makes critical the value of the adaptation variable at the time of the spike. Indeed, when a spike is emitted at time t*, the new initial condition of system 1.1 is (c, w(t*) + d). Therefore, this value governs the subsequent evolution of the membrane potential and, hence, the spike pattern produced. For instance Touboul and Brette (2008a, 2008b) show that the sequence of reset locations after each spike time shapes the spiking signature of the neuron.

In order to characterize the adaptation variable at the times of the spikes, we precisely study in section 2 the orbits of equation 1.1 in the phase plane (v, w). In the absence of a cutoff, we prove that the adaptation variable diverges at the time of the spike for the quadratic adaptive model and that it has a finite limit in the cases of the exponential and quartic models. The divergence of the adaptation variable in Izhikevich's model gives to the cutoff value Vc a very similar role as strict voltage thresholds for linear integrate-and-fire neurons. It has several theoretical and numerical consequences, which we discuss in section 3. We observe in particular that because of the divergence of the adaptation variable, the system is very sensitive to changes in the cutoff value Vc.

2.  Adaptation Variable at the Times of the Spikes

Heuristically, in the case where no adaptation is taken into account in equation 1.1 (e.g., by setting a = 0), the membrane potential variable blows up in finite time for some given initial conditions. The presence of the adaptation variable tends to slow the divergence of the membrane potential. Roughly speaking, the adaptation variable increases when v increases, which, because of the coupling, makes the derivative of the membrane potential smaller. However, we show in the sequel that the fact that the membrane potential blows up, is still true in the presence of adaptation. Therefore for some initial conditions the membrane potential blows up, and the adaptation variable either blows up or remains bounded depending on the nonlinear function F.

More precisely, in Touboul (2008), we have seen that there exists possibly one stable fixed point for system 1.1, which corresponds to a resting state of the neuron. In Touboul and Brette (2008b), we prove that all the orbits of the system that do not converge to this stable fixed point will be trapped after a finite time in a zone fully included in the half-plane {w < bv} called the spiking zone. Consider, for instance, the case where the subthreshold system has no fixed point.1 In that case, the vector field on the line w = bv points inward the half-plane {(v, w); wbv}, which implies that any trajectory having its initial condition in this zone remains trapped there. The other cases are slightly more complex and involve the description of the stable manifold of the saddle fixed point. (See Touboul & Brette, 2008b, for rigorous proofs of the existence of the spiking zone.)

Let us now consider the orbit of equation 1.1 whose initial condition at time t0 is inside the spiking zone. In this zone, since wbv, we have
formula
2.1
It is simple to prove that the solution of the equation
formula
2.2
blows up in finite time under assumption 1.2 Using Gronwall's (1919) theorem, we conclude that v(t) ⩾ u(t), and hence v blows up in finite time.
Let us now study the orbit of a solution (v(t), w(t)) of the differential system 1.1 having its initial condition in the spiking zone and characterize the behavior of w(t) as a function of v(t). In the spiking zone, we have seen that w(t) ⩽ bv(t), and therefore F(v) − w + IF(v) − bv + I, which tends to infinity when v tends to infinity. Since v(t) blows up, for any k > 0, there exists a time t1 ∈ [t0, t*) such that we will have F(v(t)) − w(t) + Ik > 0 for all t ∈ [t1, t*). We denote (v1v(t1), w1w(t1)). After time t1, because of this inequality, the trajectory in the phase plane can be written as the graph of a function W(v) that satisfies the equation
formula
2.3
(i.e., w(t) = W(v(t)) for t ∈ [t1, t*)). Since w(t) is increasing for t ∈ [t1, t*), we necessarily have
formula
2.4
Therefore Gronwall's theorem (Gronwall, 1919) ensures us that the solution of equation 2.3 will be lower-bounded for vv1 by the solution of the linear ordinary differential equation,
formula
2.5
which reads
formula
2.6
where . Because of assumption 1, the integrand involved in the expression of g, namely, , is integrable; hence, the function g(v) has a finite limit g(∞) when v → ∞. The exponential terms will hence converge when v → ∞. But the integral involved in the particular solution diverges in the quadratic case (or when F(v) grows slower than v2), since the integrand involved in formula 2.6 is equivalent when u → ∞ to
formula

Hence, the solution of the linear differential equation 2.5 tends to infinity when v → ∞ faster than a logarithmic function of v, and so does W(v), and hence w(t) blows up when v(t) blows up.

Let us now upper-bound the adaptation variable on the orbits of the system. Using the same notations, since w1w(t) ⩽ bv(t) for t ∈ [t1, t*), we have
formula
2.7
and hence
formula
2.8

In the case where F(u) = u2, this integral is bounded by a logarithmic function of v, and in the case where F(u) grows faster than u2+ε, this integral converges when v → ∞. Hence, using the fact that W is nondecreasing and upper-bounded, we conclude that W(v) has a finite limit when v → ∞.

Therefore, we have proved that in the case of the quadratic adaptive model, the adaptation variable blows up at the explosion time of the membrane potential variable v as a logarithmic function of v (see Figure 1h) whatever the parameters chosen. In the case of the quartic and exponential models, the adaptation variable has a finite limit at the explosion time of v. Therefore, for the quadratic adaptive model to be well defined, one has to introduce a finite voltage cutoff (or threshold) Vc, and the value of the adaptation variable at the time of the spike is simply given by W(Vc), which depends on the parameters of the system and the initial conditions. In the case of the quadratic model, it is an unbounded increasing function of Vc, whereas in the quartic and exponential models, this function converges to a finite value when Vc goes to infinity.

Figure 1:

Sensitivity of the spike patterns with respect to the cutoff value for the quadratic model for different set of parameters. Simulations done with a Euler scheme with time steps ranging from 10−4 to 10−3. (a, b) These have the same parameters {a = 0.02; b = 0.19; c = −60; d = 1.419, I = 10.25} but different cutoff values (30 and 45, respectively). A small increase of the cutoff results in a sharp transition from spiking to bursting linked with a period doubling bifurcation for the stationary sequence of reset adaptation values represented in e. (c, d) These are generated using the parameters {a = 0.1; b = 0.26; c = −60; d = 0;} with cutoff values 32.9 and 33, respectively. Changing the cutoff by an extremely small amount results in two very different global behaviors. (f) This represents the stationary sequence of reset adaptation values for {a = 0.02, b = 0.19, c = −57.7, d = 1.15, I = 10.377} and for cutoff values ranging from 20 to 100; an intricate bifurcation structure appears. (g) This shows the convergence of the firing rate to the intrinsic firing rate in the case of the quartic model, while the firing rate of the quadratic model regularly decreases to 0. (h) This is the divergence of the variable w as a function of v in semilogarithm axis.

Figure 1:

Sensitivity of the spike patterns with respect to the cutoff value for the quadratic model for different set of parameters. Simulations done with a Euler scheme with time steps ranging from 10−4 to 10−3. (a, b) These have the same parameters {a = 0.02; b = 0.19; c = −60; d = 1.419, I = 10.25} but different cutoff values (30 and 45, respectively). A small increase of the cutoff results in a sharp transition from spiking to bursting linked with a period doubling bifurcation for the stationary sequence of reset adaptation values represented in e. (c, d) These are generated using the parameters {a = 0.1; b = 0.26; c = −60; d = 0;} with cutoff values 32.9 and 33, respectively. Changing the cutoff by an extremely small amount results in two very different global behaviors. (f) This represents the stationary sequence of reset adaptation values for {a = 0.02, b = 0.19, c = −57.7, d = 1.15, I = 10.377} and for cutoff values ranging from 20 to 100; an intricate bifurcation structure appears. (g) This shows the convergence of the firing rate to the intrinsic firing rate in the case of the quartic model, while the firing rate of the quadratic model regularly decreases to 0. (h) This is the divergence of the variable w as a function of v in semilogarithm axis.

3.  Consequences

The divergence of the adaptation variable at the times of the spikes has a significant impact on the theoretical, qualitative, and computational analysis of the quadratic model.

First, models such as the quartic and the exponential adaptive models for which the adaptation variable converges at the times of the spikes can be defined without introducing a cutoff value. We call intrinsic model the system where the interspike behavior is governed by equation 1.1 and where spikes are emitted at the explosion time of the membrane potential variable. At these times, the membrane potential is reset, and the value of the adaptation variable at the explosion time is increased by the spike-triggered adaptation parameter d. The intrinsic system is mathematically well posed and tractable (see Touboul & Brette 2008b). It has the great interest of not requiring any voltage cutoff, as real neurons do (see Izhikevich, 2007). Numerical simulations of the model necessitate stopping the integration during the rise phase of the spike emission at a finite cutoff Vc. The behavior of the system with cutoff mimics the behavior of the intrinsic system provided the cutoff value is chosen large enough. The times of the spikes, the value of the adaptation variable at these times, and therefore the spike patterns fired are robust to the choice of the cutoff value because of the convergence property of the adaptation variable.

In the case of the Izhikevich quadratic adaptive integrate-and-fire model, since the adaptation variable blows up, the model necessarily requires the introduction of a finite voltage cutoff. Changing the cutoff value results in sensitively changing the value of the adaptation variable at the times of the spikes because of the divergence property of the adaptation variable. Since the membrane potential blows up in finite time, the time of the first spike emitted by the neuron will not be very sensitive to changes in the cutoff value. But the after-spike reset point (c, W(Vc) + d) will significantly change when varying Vc. The whole subsequent evolution of the system is therefore affected as soon as the second spike is emitted. Thus, the spike pattern produced depends on the cutoff value. This value is therefore a critical parameter of the quadratic model.

Indeed, small changes in the cutoff value dramatically affect the spike patterns produced. We illustrate the sensitivity to the cutoff value of the quadratic model using the same parameters as in Izhikevich (2003) in Figure 1. Figures 1e and 1f represent the stationary sequence of reset values of the adaptation variable as functions of the cutoff. Regular spiking corresponds to the convergence of this sequence and bursts to cycles. We observe a period doubling bifurcation when varying the cutoff value in Figure 1c, that results in abruptly switching from a regular spiking behavior (see Figure 1a) to a bursting behavior (see Figure 1b). More complex bifurcation structures involving chaotic patterns also appear (see Figure 1f), and in that case, infinitesimal changes in the cutoff value result in dramatic changes in the qualitative spiking behavior. Changing the cutoff in that case makes the system switch between chaotic spiking, bursts with eight, four, and eventually two spikes. Other important properties such as bistability strongly depend on the cutoff (see Figures 1c and 1d).

Besides these bifurcations, the divergence of the adaptation variable as a function of the cutoff Vc induces a slow, continuous decrease to zero of the firing rate of the neuron, whereas the firing rate converges to the related intrinsic model's firing rate for the quartic and the exponential model (see Figure 1g).

Because of this sensitivity, the cutoff value has to be carefully set in order to qualitatively and quantitatively fit experimental data sets. Since it has no clear biophysical interpretation, the meaning of the cutoff Vc in the quadratic model and the problem of its accurate evaluation have to be specifically addressed.

Eventually, from a numerical viewpoint, the unboundedness of the adaptation variable and its time derivative at the explosion times of the membrane potential makes the accurate computation of this value difficult in the quadratic adaptive model. In particular, the time step necessary to accurately estimate this value has to be very small (or to be adaptive as a function of the value of the membrane potential variable) in order to obtain the right spike pattern. This questions the efficiency of numerical simulations making use of this particular model.

4.  Conclusion

In this note, we proved that the adaptation variable of the adaptive quadratic model blew up at the times of the spikes, whereas it converged for the quartic and the adaptive exponential models. This property has some important implications that we discussed. From a theoretical point of view, we showed that the nature of the spike patterns produced undergoes bifurcations with respect to the cutoff value, and this made the system very sensitive to this parameter: small changes in the value of this parameter can deeply affect the nature of the spiking pattern. From a quantitative viewpoint, it raises the question of how to evaluate this cutoff in order to fit data sets, and from a numerical viewpoint, it has implications on the efficiency of the simulation algorithms to use. The convergence of this value for models having a faster blow-up at the times of the spike, such as the quartic or the exponential adaptive models, implies that the system presents intrinsic spiking properties that can be mathematically studied and efficiently simulated and are robust to changes in the cutoff value.

Acknowledgments

I warmly acknowledge Bruno Cessac, Olivier Faugeras, and François Grimbert for interesting discussions and suggestions.

Notes

1

Or in the case where the initial condition (v, w) is such that v is greater than the largest v-value of the fixed points (the largest solution of F(v) − bv + I = 0) and wbv.

2

For the quadratic model, we can get analytic expressions of the solutions involving the tangent function and therefore can derive an upper bound of the explosion time.

References

Badel
,
L.
,
Lefort
,
S.
,
Brette
,
R.
,
Petersen
,
C. C. H.
,
Gerstner
,
W.
, &
Richardson
,
M. J. E.
(
2008
).
Dynamic I-V curves are reliable predictors of naturalistic pyramidal-neuron voltage traces
.
Journal of Neurophysiology
,
99
,
656
66
.
Brette
,
R.
, &
Gerstner
,
W.
(
2005
).
Adaptive exponential integrate-and-fire model as an effective description of neuronal activity
.
Journal of Neurophysiology
,
94
,
3637
3642
.
Clopath
,
C.
,
Jolivet
,
R.
,
Rauch
,
A.
,
Lüscher
,
H.
, &
Gerstner
,
W.
(
2007
).
Predicting neuronal activity with simple models of the threshold type: Adaptive exponential integrate-and-fire model with two compartments
.
Neurocomputing
,
70
(
10–12
),
1668
1673
.
Fourcaud-Trocme
,
N.
,
Hansel
,
D.
,
van Vreeswijk
,
C.
, &
Brunel
,
N.
(
2003
).
How spike generation mechanisms determine the neuronal response to fluctuating inputs
.
Journal of Neuroscience
,
23
(
37
),
11628
.
Gronwall
,
T. H.
(
1919
).
Note on the derivative with respect to a parameter of the solutions of a system of differential equations
.
Annals of Mathematics
,
20
,
292
296
.
Izhikevich
,
E.
(
2003
).
Simple model of spiking neurons
.
IEEE Transactions on Neural Networks
,
14
(
6
),
1569
1572
.
Izhikevich
,
E.
(
2004
).
Which model to use for cortical spiking neurons?
IEEE Transactions on Neural Networks
,
15
(
5
),
1063
1070
.
Izhikevich
,
E.
(
2007
).
Dynamical systems in neuroscience: The geometry of excitability and bursting
.
Cambridge, MA
:
MIT Press
.
Izhikevich
,
E. M.
, &
Edelman
,
G. M.
(
2008
).
Large-scale model of mammalian thalamocortical systems
.
Proc. Natl. Acad. Sci. U.S.A.
,
105
(
9
),
3593
3598
.
Jolivet
,
R.
,
Kobayashi
,
R.
,
Rauch
,
A.
,
Naud
,
R.
,
Shinomoto
,
S.
, &
Gerstner
,
W.
(
2008
).
A benchmark test for a quantitative assessment of simple neuron models
.
Journal of Neuroscience Methods
,
169
(
2
),
417
424
.
Koch
,
C.
, &
Segev
,
I.
(Eds.). (
1998
).
Methods in neuronal modeling: From ions to networks
.
Cambridge, MA
:
MIT Press
.
Lapicque
,
L.
(
1907
).
Recherches quantitatifs sur l'excitation des nerfs traitée comme une polarisation
.
J. Physiol. Paris
,
9
,
620
635
.
Latham
,
P.
,
Richmond
,
B.
,
Nelson
,
P.
, &
Nirenberg
,
S.
(
2000
).
Intrinsic dynamics in neuronal networks. I. Theory
.
Journal of Neurophysiology
,
83
,
828
835
.
Rinzel
,
J.
, &
Ermentrout
,
B.
(
1989
).
Analysis of neural excitability and oscillations
.
Cambridge, MA
:
MIT Press
.
Stein
,
R. B.
(
1967
).
Some models of neuronal variability
.
Biophysical Journal
,
7
,
37
68
.
Touboul
,
J.
(
2008
).
Nonlinear and stochastic models in neuroscience
.
Unpublished doctoral dissertation, Ecole Polytechnique
.
Touboul
,
J.
, &
Brette
,
R.
(
2008a
).
Dynamics and bifurcations of the adaptive exponential integrate-and-fire model
.
Biological Cybernetics
,
99
(
4–5
),
319
334
.
Touboul
,
J.
, &
Brette
,
R.
(
2008b
).
Spiking dynamics of bidimensional integrate-and-fire neurons
(
Research Rep. 6531
).
Sophia Antipolis
:
INRIA
.