|Home | About | Journals | Submit | Contact Us | Français|
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.
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.
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
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 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:
and the spike rules:
where 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
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.
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: 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:
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 (the neuron fires) at point Pk = (h,yk), then according to the spike rules (4)–(5), the trajectory jumps to the point due to the impulsive effect (see Fig. 2). After that, the trajectory evolves according to (2)–(3), and it could reach again the section at a point Pk+1 = (h,yk+1). Since 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:
and some straightforward calculations yield:
that defines an implicit map for the adaptive variable at the firing times. Assuming that ∀k, yk < a + h2 and setting L = (2a + h2 + q2−b)(h2−q2), H = a + h2, Q = p−a−q2, then we can define the map f
that gives the successive values of the adaptive variable at the firing times.
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:
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:
which yields the following two fixed points:
Moreover, according to (7) we have
where Df represents the derivative of f. The function Df is decreasing if
since c > 0. In addition, by defining 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 and a necessary and sufficient condition for the existence of yc is
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
By construction, for all y 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
and its neighborhood Br(y*).
According to map (7), we have which yields two possible solutions: and The first solution is identically equal to its corresponding fixed point, hence the possibility of existence of a snap-back repeller ym Br(y*) is relative to the second one, noted as y1. Obviously we have y1 < y*.
Assume that the following conditions are satisfied:
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 ]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 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 [yz,y*] such that y2 = f(y3). Again if y3 Br(y*), then we found m = 3.
Otherwise, let us firstly consider With f(y2) = y1, f(y*) = y*, we get 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 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:
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 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
Finally, since , hence if
we have Dfm(yk) ≠ 0. It is worth noting that if there exist some points such that , 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.
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:
and the spike rules
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:
Numerically we check that inequality (8) is satisfied and straightforward calculations yield the following values for the variables introduced in the proof:
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 Br(y*), condition (i) of Theorem 1 is fulfilled.
It is easy to check that, for a y4≈12.6150 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.
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).
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.
Gang Zheng, Phone: +33-1-30736666, Fax: +33-1-30736627, Email: rf.aesne@gnehZ.gnaG.
Arnaud Tonnelier, Email: rf.seplairni@reilennoT.duanrA.