PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of cogneurospringer.comThis journalToc AlertsSubmit OnlineOpen Choice
 
Cogn Neurodyn. 2009 September; 3(3): 197–204.
Published online 2008 November 6. doi:  10.1007/s11571-008-9069-6
PMCID: PMC2727159

Chaotic solutions in the quadratic integrate-and-fire neuron with adaptation

Abstract

The quadratic integrate-and-fire (QIF) model with adaptation is commonly used as an elementary neuronal model that reproduces the main characteristics of real neurons. In this paper, we introduce a QIF neuron with a nonlinear adaptive current. This model reproduces the neuron-computational features of real neurons and is analytically tractable. It is shown that under a constant current input chaotic firing is possible. In contrast to previous study the neuron is not sinusoidally forced. We show that the spike-triggered adaptation is a key parameter to understand how chaos is generated.

Keywords: QIF, Chaos, Poincaré map, Lyapunov exponent

Introduction

Neurons transform incoming stimuli into train of all-or-one electric events known as spikes. Precise spike times are believed to play a fundamental role in the encoding of information in the brain. Recent experimental evidences have accumulated suggesting that precise spike-time coding is used in various neuronal systems (VanRullen et al. 2005). In this view, the spike time reliability is a key to understand the basis of the neural code (Mainen and Sejnowski 1995). It has been shown that the responses of neurons to time varying stimuli have high reliability with a precision of 1 ms or less (Mainen and Sejnowski 1995; Berry et al. 1997). On the other hand, constant stimuli elicits spike trains with high variability (de Ruyter van Steveninck et al. 1997). Determining the variability of spike timing can provide fundamental insights into the strategies used in neural systems to represent and transmit information. However the origin of the neuron’s response variability is still debated. It has been primarily suggested that such variability in discharge times is due to accumulating noise. However, it could be due as well to a chaotic behavior (Hayashi et al. 1982; Aihara et al. 1984; Chacron et al. 2004). In fact, the responses of a neuron to given stimuli over repeated trials could be different because of intrinsic noise but also because of the chaotic behavior of the neuronal dynamics. It remains unclear if these different sources of variability coexist, what is the main cause of non-reproducibility in neuronal responses under constant stimuli and what are the consequences on the neural code.

The reproducibility of spikes has been largely investigated in theoretical studies. Simplified phenomenological neuron models of integrate-and-fire type are commonly used to explore the dynamical properties of neurons and particularly the robustness of neuron responses (Tiesinga 2002; Brette and Guigon 2003). Chaotic solutions have been investigated for the leaky integrate-and-fire neuron (Coombes 1999; Chacron et al. 2004). The time-varying injected current combined with a periodic modulation of parameters are the minimal requirements to obtain chaotic solutions.

The traditional leaky integrate-and-fire neuron has some limitations and drawbacks: it has an unrealistic behavior close to the threshold and cannot reproduce accurately the frequency response of real neurons. An extension to quadratic integrate-and-fire (QIF) neurons allows to describe the nonlinear spike-generating currents of biological neurons and captures the so-called type I firing behavior near the bifurcation (Ermentrout 1996). In this paper we show that adaptive QIF neuron can exhibit chaotic solutions under constant input current. This result corroborates the non-reliability of spike times observed experimentally when a constant current is injected into a neuron (Mainen and Sejnowski 1995; de Ruyter van Steveninck et al. 1997).

The paper is organized as follows. In the next section the model is presented. We derive in section “Poincaré map” the expression of the Poincaré map associated with the firing of a spike and use the Marotto theorem in section “Chaotic behavior” to show that the adaptive QIF can exhibit chaotic dynamics. In fifth section “Numerical simulations” are made to support our analysis and to highlight the dynamical behavior of the model.

The model

Simplified neuron models are widely used to describe the electrical activity of biological spiking neurons. Among them, integrate-and-fire models have become very popular probably due to their capabilities to accurately reproduce the dynamics of real neurons through a low-dimensional dynamical system that can be easily tuned by varying few parameters. The QIF neuron has seen increasing interest in recent years, primarily because it reproduces the properties of detailed conductance-based neurons near the threshold. In the QIF model, the membrane potential follows

equation M1
1

where x is the membrane potential and a is a constant input current. Because of the quadratic term in (1) the membrane potential can escape to infinity in finite time. We thus introduce h the cutoff of a spike that defines the firing of a spike. When x(t) reaches h a spike is produced and x(.) is reset to a subthreshold potential q. The quadratic neuron model presents two distinct regimes. When a < 0, there are two fixed points. The stable one defines the resting state equation M2 and the unstable one defines an ’effective’ threshold below which trajectories tend toward the peak of the spike in finite time. When a > 0 the neuron fires regularly. In this regime, the quadratic model is a canonical representation of firing behavior of so-called type-I neurons close to bifurcation (Ermentrout 1996; Ermentrout and Kopell 1986) and is related to the so-called θ neuron (Gutkin and Ermentrout 1998). The addition of a second variable, y, allows inclusion of adaptation, subthreshold resonance (Richardson et al. 2003) and captures the bursting behavior of neurons (Izhikevich 2000). We consider an adaptive quadratic integrate-and-fire model given by the following differential equations:

equation M3
2
equation M4
3

and the spike rules:

equation M5
4
equation M6
5

where equation M7 Parameters h and q are respectively the peak of a spike and the reset value. Parameters b, c and p describe the adaptive current and we assume c ≥ 0. Note that y can be seen as a recovery variable. The function τ(x) models the voltage-dependence of the adaptive variable time-constant. For a fixed voltage value x, the variable y approaches the equilibrium b with a time constant τ(x). In the following we consider τ(x) = τ/x where τ is a positive constant. This choice approximates the voltage dependence of the time constant of the adaptive variable commonly used in detailed models: for large value of x the integration time of y is shorter. Therefore the time evolution of the adaptive current is given by

equation M8

where the factor 2 is used for convenience (see below) and can be removed (change of variables).

After a spike has been triggered the membrane potential restarts at the reset value q and the variable y is multiplied by a constant c and is increased by an amount p which account for spike-triggered adaptation. The differential equations (2)–(3) model the subthreshold dynamics and the spike initiation. When x attains the peak, the impulsive effect occurs according to (4)–(5). The time evolution is illustrated in Fig. 1 for parameter values where the model exhibits regular spiking. Simulations (not shown) indicate that classical neuro-computational features of neurons (Izhikevich 2004) including bistability, resonator and bursting (see later) can be reproduced by this model. Note that the proposed model differs from the spiking model of Izhikevich (2003) because of the non-linearity in the equation for the adaptive current. Despite extensive studies on nonlinear integrate-and-fire models with adaptive current, including the adaptive quadratic (Ermentrout et al. 2001; Izhikevich 2007), the adaptive exponential (Brette and Gerstner 2005) or related models (Touboul 2008), chaotic solutions have not been reported, to our knowledge.

Fig. 1
Illustrative plots of the time evolution for (a) the membrane voltage x(t) and (b) the adaptive variable y(t) in the regular spiking regime. Parameters are a = 2, τ = 15, b = 1, q = 1, ...

Poincaré map

Due to the non smooth nature of integrate-and-fire models, the standard theory of dynamical systems does not apply directly. This problem was overcome in previous studies through the use of integrate-and-fire models where the map of successive firing times can be studied directly (Brette 2008) or analytically derived (Coombes 1999; Chacron et al. 2004). However, because of the nonlinearity of the subthreshold dynamics of the QIF it is no longer possible here to derive an analytical expression for this map. In order to analyze the behavior of system (2)–(5), we introduce the following Poincaré section: equation M9 It is easy to show that the differential system (2)–(3) is conservative. In fact, if we define the energy of the system as follows:

equation M10
6

we have dE/dt = 0 along the trajectories of the system provided that no firing occurs. Let us now consider a trajectory that reaches the section equation M11 (the neuron fires) at point Pk = (h,yk), then according to the spike rules (4)–(5), the trajectory jumps to the point equation M12 due to the impulsive effect (see Fig. 2). After that, the trajectory evolves according to (2)–(3), and it could reach again the section equation M13 at a point Pk+1 = (h,yk+1). Since equation M14 and Pk+1 are both located on the trajectory of system (2)–(3) where no impulsive effect occurs, these two points satisfy relation (6). Consequently, we obtain:

equation M15

and some straightforward calculations yield:

equation M16

that defines an implicit map for the adaptive variable at the firing times. Assuming that ∀k, yk < a + h2 and setting L = (2a + h2 + q2b)(h2q2), H = a + h2, Q = paq2, then we can define the map f

equation M17
7

that gives the successive values of the adaptive variable at the firing times.

Fig. 2
A piece of trajectory containing a jump (spike) in the phase plane (x,y). The vertical dashed lines are the reset line x = q and the Poincaré section that corresponds to the threshold line x = h

Chaotic behavior

In this section, with the help of Marotto’s Theorem (Marotto 1978, 2004) we prove that the adaptive QIF neuron (2)–(5) can exhibit chaotic behavior for specific values of parameters. The main idea of Marotto’s theorem is to seek a snap-back repeller. Marotto showed that the presence of a snap-back repeller is a sufficient criterion for the existence of chaos. Let us first recall the basic Marotto’s Theorem for one-dimensional system:

Theorem 1 (Marotto 1978) Lety*be a fixed point of the mapf. This fixed point is said to be a snap-back repeller if the following conditions are satisfied:

  1. fis differentiable in a neighborhoodBr(y*) of the fixed pointy*with radiusr > 0 and the eigenvalue ofDf(y) is strictly larger than one in absolute value for ally [set membership] Br(y*).
  2. There exists a pointym [set membership] Br(y*), withym ≠ y*such that for some positive integerm, fm(ym)  = y*, andfmis differentiable atymwithDfm(ym) ≠ 0.

A snap-back repeller arises when the fixed point is repelling while there is a trajectory starting in the repelling neighborhood of the fixed point that goes back to the fixed point, i.e. the repelling fixed point has an associated homoclinic orbit. Note that the existence of a snap-back repeller for one-dimensional system is closely related to the existence of a periodic point with period 3. In fact, it can be shown that the existence of a snap-back repeller of f is equivalent to the existence of a point of period 3 for the map fn for some positive integer n (Li and Yorke 1975).

Using Theorem 1 we will prove the following

Proposition 1 For appropriate values of parameters the mapfis chaotic in the sense of Marotto.

In the remainder of this section we prove that the map f fulfills the conditions of Marotto theorem. In order to identify a snap-back repeller for the one dimensional map (7), we derive the eigenvalue of the system and show that, for values of the state variable close to the fixed point, the corresponding eigenvalues lie outside the unit circle, whereas for values of the state variable that are sufficiently far from the fixed point, the eigenvalues are within the unit circle.

Let us firstly check whether condition (i) in Theorem 1 is satisfied. From (7) it is easy to show that the map has a fixed point y* = f(y*) if the following inequality is fulfilled:

equation M18
8

which yields the following two fixed points:

equation M19
9

Moreover, according to (7) we have

equation M20

where Df represents the derivative of f. The function Df is decreasing if

equation M21
10

since c > 0. In addition, by defining equation M22 it is easy to check that f(y) is increasing if y < yz and is decreasing when y > yz.

Let yc be the solution of the equality: Df(yc) = −1. Since equation M23 and equation M24 a necessary and sufficient condition for the existence of yc is

equation M25
11

Moreover since Df(yc) < Df(yz) and Df(y) is decreasing when y > yz, we have yc > yz.

Let us now consider a neighborhood of y*, noted as Br(y*), with

equation M26
12

By construction, for all y [set membership] Br(y*), Df(y) is always strictly larger than one in absolute values provided that y* > yc > yz. Therefore, condition (i) of Theorem 1 is satisfied.

On the other hand, in order to fulfill condition (ii) of Theorem 1, we seek a point ym ≠ y* in the neighborhood Br(y*), such that ym−1 = f(ym), ..., y1 = f(y2) = fm−1(ym), and y* = f(y1) = fm(ym), with some positive integer m, and Df(yk) ≠ 0, for 1 ≤ k ≤ m.

Let us consider the following fixed point

equation M27
13

and its neighborhood Br(y*).

According to map (7), we have equation M28 which yields two possible solutions: equation M29 and equation M30 The first solution is identically equal to its corresponding fixed point, hence the possibility of existence of a snap-back repeller ym [set membership] Br(y*) is relative to the second one, noted as y1. Obviously we have y1 < y*.

Assume that the following conditions are satisfied:

equation M31
14

Since f(yA) = yB and f(y*) = y*, and f(y) is a decreasing continuous function over [y*,yA], therefore there always exists a point y2 [set membership] ]y*,yA], such that y1 = f(y2). If such a y2 is located in Br(y*), we have m = 2. Otherwise, consider the interval ]y1,yz[, since y* = f(y1) and yA = f(yz), and f(y) is increasing over [y1,yz], it is always possible to find a equation M32 such that y2 = f(y3). It should be noted that there exists another possible y3 satisfying y2 = f(y3), which is located over [yz,y*]. Notice that y* = f(y*), yA = f(yz), and f(y) is decreasing over [yz, y*], thus we can find the second y3 [set membership] [yz,y*] such that y2 = f(y3). Again if y3 [set membership] Br(y*), then we found m = 3.

Otherwise, let us firstly consider equation M33 With f(y2) = y1, f(y*) = y*, we get equation M34 such that f(y4) = y3. If this y4 is not yet located in Br(y*), according to f(y3) = y2, f(yz) = yA, f(y*) = y*, since f(y) is increasing over [y3,yz] but decreasing over [yz,y*], hence it always exists a equation M35 such that f(y5) = y4. By induction, if the found points are not yet entered into Br(y*), the above procedure generates the sequence: yB < y1 < y3 < yz < y5 < ··· < y* < ··· < y4 < y2 < yA < H, which tends towards the equilibrium point y*.

In order to prove the convergence of the above sequence to the equilibrium point y*, let us analyze the right side of the sequence (i.e., even subscript), which in fact can be written as follows:

equation M36

with k ≥ 1. When k → ∞, we have y2k = y2k+2, noted as ylim. Obviously ylim is one solution of equation: ylim = f2(ylim), which in fact has two solutions under the constrain y* < y2k+2 < y2k < H. It is easy to check that equilibrium points of (9) are these two solutions. However, since our analysis is focused only on the equilibrium point of (13), thus we have ylim = y*. It should be noticed that, for the left side of the sequence (odd subscript), we can obtain the same result through a similar process, which implies the sequence converges towards the equilibrium point y* via both sides of y*. Consequently, we can always find a point ym [set membership] Br(y*), such that fm(ym) = y* with some positive integer m.

With the same procedure, it is easy to adapt the above proof for the second possible value of equation M37

Finally, since equation M38, hence if

equation M39
15

we have Dfm(yk) ≠ 0. It is worth noting that if there exist some points such that equation M40, then they correspond to a periodic orbit of (7).

In summary, provided that the technical conditions: (8), (10)–(11), (14)–(15) are fulfilled, it is possible to determine a neighborhood Br(y*) around the fixed point (13) with r defined in (12) such that conditions (i) and (ii) are satisfied. Therefore, Theorem 1 holds and this ends the proof of Proposition 1.

Numerical simulations

Besides the theoretical analysis, numerical simulations are made in order to show that system (2)–(5) can exhibit chaotic dynamics. In particular, we show that there exist parameter values for which the technical assumptions used in our proof are satisfied.

We consider the following QIF neuron with adaptation:

equation M41
16

and the spike rules

equation M42
17

i.e. we choose the following set of parameter values: a = 6, b = 2, τ = 1, p = −0.2, q = 10 and h = 20. We explore the influence of parameter c on the asymptotic behavior of the model and we investigate the existence of adaptation-induced chaotic dynamics using different values for c.

For c = 13.8, we have L = 153,000, H = 406 and Q = −106.2 that give the corresponding Poincaré map:

equation M43
18

Numerically we check that inequality (8) is satisfied and straightforward calculations yield the following values for the variables introduced in the proof:

equation M44

It is easy to check that conditions (10)–(11) and (14)–(15) are all satisfied. Moreover we choose r = 1.5 < y*yc that defines the neighborhood Br(y*)  = ]9.9434, 12.9434[ (see Fig. 3a). Since we have Df(9.9434) = −1.0909 and Df(12.9434) = −2.5123, thus for all yk [set membership] Br(y*), condition (i) of Theorem 1 is fulfilled.

Fig. 3
Snap-back repeller of the QIF neuron with adaptation. The numerical investigations for the existence of a snap-back repeller are done in (a) and (b), respectively. a The derivative of the Poincaré map (18). The successive iterations yk are represented ...

It is easy to check that, for a y4≈12.6150 [set membership] Br(y*) and y4 ≠ y*, we have y3 = f(y4) ≈ 9.0005, y2 = f(y3) ≈ 14.4336, y1 = f(y2) ≈ 3.9479 and y* = f(y1) (see Fig. 3b). Moreover Df4 (y4) ≈ −8.6461 ≠ 0. Hence condition (ii) of Theorem 1 is satisfied. Consequently, the point y* is a snap-back repeller and the map (7) exhibits chaotic behavior.

Besides the existence of chaos, the scenario leading to such dynamics as c is varied is illustrated in Fig. 4. The phase portrait of system (16)–(17) is shown together with the membrane potential as a function of time. When c = 10, the system tends towards an asymptotically stable periodic orbit (See Fig. 4a, b). In Fig. 4c, d, for c = 13.9, a stable period 3 solution of the Poincaré map occurs. A chaotic behavior can be generated by choosing c = 13.8 as shown in Fig. 4e, f.

Fig. 4
Transition to chaos through a cycle of period 3 in the QIF neuron with adaptation. The phase portrait are given in the left column and the x dynamics towards time are depicted in the right column. a, b Asymptotically stable periodic orbit for system with ...

To gain deeper understanding of the chaotic dynamics, the bifurcation diagram obtained while varying parameter c is shown Fig. 5a, while Fig. 5b shows the leading Lyapounov exponent of the system. The lines in the bifurcation diagram correspond to periodic spiking and the thick parts are associated with regimes where the adaptive variable takes either a countable but very large number of values or an uncountable value. It is not clear that in the bifurcation diagram the regimes where the adaptive variable takes a very large number of values are associated with periodic orbits with large period or chaotic regimes. One of the usual ways to test for chaos is to compute system’s largest lyapunov exponent. The Lyapunov exponent of a dynamical system is a quantity that characterizes the rate of separation of infinitesimally close trajectories. The general idea is to follow two nearby orbits and to calculate their average logarithmic rate of separation. A positive Lyapunov exponent indicates chaos. According to expression (18), it is easy to calculate the largest lyapunov exponent and the result is shown Fig. 5b. The range of positive exponents suggests chaotic regimes within these regions. Moreover the bifurcation diagram reveals that the system exhibits transition to chaos through period doubling bifurcation cascade. Because system (18) is chaotic, then its corresponding continuous system is also chaotic (a non-chaotic trajectory of a continuous system, for instance a periodic orbit, will never generate a chaotic Poincaré map).

Fig. 5
a Adaptive variable at the firing time, yk, as a function of c illustrating the bifurcation diagram of our model. b The corresponding variation of the Lyapunov exponent with respect to parameter c. Simulation was made with the initial condition: x0 = 5, ...

Conclusion

Periodically forced integrate-and-fire neurons have been proposed to produce chaotic firing (Chacron et al. 2003, 2004). It is expected that chaos in spiking neuronal models can also be generated through an applied constant current in agreement with the experimentally observed neural spike trains variability in this condition. This paper presents a simple nonlinear integrate-and-fire neuron which can exhibit a chaotic behavior under a constant current. The model is the quadratic integrate-and-fire neuron with a nonlinear adaptive current that takes into account the voltage dependence of the adaptive-variable time constant. The situation significantly differs from those in previous studies that report chaos in integrate-and-fire models for several reasons: (i) there is no periodic forcing (ii) threshold or resetting voltage are not periodically modulated and (iii) the integrate-and-fire model is nonlinear. We proved mathematically the existence of chaotic solutions provided that parameters of the adaptive current are suitably chosen. Numerical simulations are given including the largest Lyapunov exponent and the bifurcation diagram.

A further development of this work is the study of entrainment in sinusoidally driven adaptive QIF neuron. It is expected that the model exhibits specific dynamics in response to periodic forcing that do not exist in LIF models. Further research is needed to clarify this point.

Contributor Information

Gang Zheng, Phone: +33-1-30736666, Fax: +33-1-30736627, rf.aesne@gnehZ.gnaG.

Arnaud Tonnelier, rf.seplairni@reilennoT.duanrA.

References


  • Aihara K, Matsumoto G, Ikegaya Y (1984) Periodic and non-periodic responses of a periodically forced hodgkin-huxley oscillator. J Theor Biol 109:249–269 [PubMed]

  • Berry MJ, Warland DK, Meister M (1997) The structure and precision of retinal spike trains. Proc Natl Acad Sci USA 94:5411–5416 [PubMed]

  • Brette R (2008) The cauchy problem for one-dimensional spiking neuron models. Cogn Neurodyn 2:21–27 [PMC free article] [PubMed]

  • Brette R, Gerstner W (2005) Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J Neurophysiol 94:3637–3642 [PubMed]

  • Brette R, Guigon E (2003) Reliability of spike timing is a general property of spiking model neurons. Neural Comput 15:279–308 [PubMed]

  • Chacron MJ, Pakdaman K, Longtin A (2003) Interspike interval correlations, memory, adaptation, and refractoriness in a leaky integrate-and-fire model with threshold fatigue. Neural Comput 15:253–278 [PubMed]

  • Chacron MJ, Longtin A, Pakdaman K (2004) Chaotic firing in the sinusoidally forced leaky integrate-and-fire model with threshold fatigue. Physica D 192:138–160

  • Coombes S (1999) Liapunov exponents and mode-locked solutions for integrate-and-fire dynamical systems. Phys Lett A 255(1–2):49–57

  • de Ruyter van Steveninck RR, Lewen GD, Strong SP, Koberle R, Bialek W (1997) Reproducibility and variability in neural spike trains. Science 275:1805–1808 [PubMed]

  • Ermentrout B (1996) Type I membranes, phase resetting curves, and synchrony. Neural Comput 8:979–1001 [PubMed]

  • Ermentrout B, Kopell N (1986) Parabolic bursting in an excitable system coupled with a slow oscillation. SIAM J Appl Math 46:233–253

  • Ermentrout B, Pascal M, Gutkin B (2001) The effects of spike frequency adaptation and negative feedback on the synchronization of neural oscillators. Neural Comput 13:1285–1310 [PubMed]

  • Gutkin B, Ermentrout B (1998) Dynamics of membrane excitability determine interspike interval variability: A link between spike generation mechanisms and cortical spike train statistics. Neural Comput 10:1047–1065 [PubMed]

  • Hayashi H, Ishizuka S, Ohta M, Hirakawa K (1982) Chaotic behavior in the onchidium giant neuron under sinusoidal forcing. Phys Lett A 88:435–438

  • Izhikevich EM (2000) Neural excitability, spiking, and bursting. Int J Bifurcat Chaos 10:1171–1266

  • Izhikevich EM (2003) Simple model of spiking neurons. IEEE Trans Neural Netw 14:1569–1572 [PubMed]

  • Izhikevich EM (2004) Which model to use for cortical spiking neurons? IEEE Trans Neural Netw 15:1063–1070 [PubMed]

  • Izhikevich EM (2007) Dynamical systems in neuroscience: the geometry of excitability and bursting. MIT, Cambridge, USA

  • Li T-Y, Yorke JA (1975) Period three implies chaos. Am Math Mon 82:985–992

  • Mainen ZF, Sejnowski TJ (1995) Reliability of spike timing in neocortical neurons. Science 268:1503–1506 [PubMed]

  • Marotto FR (1978) Snap-back repellers imply chaos in Rn. J Math Anal Appl 63:199–223

  • Marotto FR (2004) On redefining a snap-back repeller. Chaos Solitons Fractals 25:25–28

  • Richardson M, Brunel N, Hakim V (2003) From subthreshold to firing-rate resonance. J Neurophysiol 89:2538–2554 [PubMed]

  • Tiesinga PHE (2002) Precision and reliability of periodically and quasiperiodically driven integrate-and-fire neurons. Phys Rev E 65:41913 [PubMed]

  • Touboul J (2008) Bifurcation analysis of a general class of non-linear integrate and fire neurons. SIAM Appl Math 4:1045–1079

  • VanRullen R, Guyonneau R, Thorpe SJ (2005) Spike times make sense. Trends Neurosci 28:1–4 [PubMed]

Articles from Cognitive Neurodynamics are provided here courtesy of Springer Science+Business Media B.V.