1. science
2. physics

CHIN. PHYS. LETT. Vol. 32, No. 4 (2015) 040201
Bursting Behavior in the Piece-Wise Linear Planar Neuron Model with Periodic
Stimulation *
JI Ying(季颖), WANG Ya-Wei(王亚伟)**
Faculty of Science, Jiangsu University, Zhenjiang 212013
(Received 12 December 2014)
A piece-wise linear planar neuron model with periodic stimulation which can mimic the behavior of bursting is
explored. The periodic bursting with three frequencies can be observed in numerical simulation. We present
an analysis of bursting in this non-smooth non-autonomous system by considering the system as a generalized
autonomous system and introduce an appropriate form of the associated generalized equilibrium solution. The
bifurcation mechanism of bursting as well as the coexistence of three frequencies is investigated in detail.
PACS: 02.30.Oz, 02.60.Cb, 05.45.Pq
DOI: 10.1088/0256-307X/32/4/040201
The biological nervous system is a very complex
information network composed of innumerable coupling neurons. Neurons can encode, transfer, and integrate information by firing activities, which are mainly
embodied in generation and processing of action potentials (APs).[1] Bursting is one of the most important firing activities of neuronal systems which can
exhibit a transition between a rest state and a spiking state.[2−5] Due to the importance of neuronal electrical activities, numerous mathematical models[6,7]
as well as some standard mathematical methods[8,9]
have been established to promote theoretical investigations. From a modeling perspective, a series of
single neuron models have been built, which can accurately reflect corresponding AP shapes for different neuronal cell types. Typically such models, being based around that of Hodgkin–Huxley,[6] are highdimensional and usually need to be analyzed by using
the singular perturbation theory. On the other hand, a
class of one-dimensional integrate-and-fire type model
is often advocated since it is easy to analyze.[10] However, both cases have disadvantages.[11] Thus planar
models possessing one voltage and one gating variable that can mimic the behavior of high-dimensional
conductance based models are of interest.[12] Inspired
by the FitzHugh–Nagumo model, a class of piece-wise
linear (PWL) models that can describe the rich dynamical behavior of many common cell types has attracted wide attention.[13] As far as theoretical analysis is concerned, many scholars have revealed the
inherent dynamical nature of bursting by some crucial bifurcations.[8,9] On the basis of these works, the
bifurcation mechanism of the bursters has been well
summarized by Izhikevich, who classified almost all
the bursters in systems with low dimensions.[14] However, up to now, most of the work has been confined
to smooth autonomous systems.[15−17] The bursting
phenomena of non-smooth non-autonomous systems
still have problems.
In this Letter, a modified two-dimensional McKean
model of a single neuron activity with a periodic drive
is explored. The McKean model is a planar piecewise linear model that can mimic the firing response
of several different cell types.[18] As a focus of theoretical research, it is usually modified to different forms.
One can be written as[19]
𝐶 𝑣˙ = 𝑓 (𝑣) − 𝐶𝛾 − 𝐶𝑤 + 𝐼,
𝑤˙ = 𝑣 − 𝛾𝑓 (𝑣) − 𝛾𝐼,
where 𝑣 stands for the
gating variable, and
⎧
⎨−𝑣,
𝑓 (𝑣) = 𝑣 − 𝑎,
⎩
1 − 𝑣,
(1)
membrane potential, 𝑤 is the
𝑣 < 𝑎/2,
𝑎/2 ≤ 𝑣 ≤ (1 + 𝑎)/2,
𝑣 > (1 + 𝑎)/2.
Here 𝐶 > 0, 𝛾 > 0, 𝑓 (𝑣) is a PWL approximation
of the cubic FitzHugh–Nagumo nonlinearity 𝑓 (𝑣) =
𝑣(1 − 𝑣)(𝑣 − 𝑎), provided that 0 < 𝑎 < 1, 𝐼 is an external drive.[19] For the neuronal systems, the stimulation may be constant or not; for example, a periodic
stimulation. Thus the external drive is in the form of
𝐼 = cos(𝜔0 𝑡) in this work.
When there is an order gap between the stimulation frequency and the natural frequency of the system, the associated dynamics may exhibit a fast-slow
effect; that is, bursting. To reveal the mechanism
of the movements, the above parameters are fixed at
𝐶 = 2.0, 𝑎 = 0.25, 𝛾 = 0.65, meanwhile 𝜔0 is set as
𝜔0 = 0.05. The periodic bursting can be observed by
numerical simulation (see Fig. 1).
During each longest period, which is the same as
the period of the external stimulation, another two
frequencies could be observed from the time series,
which are much higher than the frequency of stimulation. The combination of the oscillations associated
with these two high frequencies and the oscillation
corresponding to the stimulation forms the bursting
behavior in this system. It is noteworthy that the coexisting phenomenon of three frequencies is not very
common in bursting. In the following, we investigate
* Supported by the National Natural Science Foundation of China under Grant Nos 11302086, 11374130, 11474134 and 11402226,
the Natural Science Foundation of Jiangsu Province under Grant No BK20141296, the Post-doctoral Science Fund of China under
Grant No 2014M561574, and the Post-doctoral Science Fund of Jiangsu Province under Grant No 1302094B.
** Corresponding author. Email: [email protected]
© 2015 Chinese Physical Society and IOP Publishing Ltd
040201-1
CHIN. PHYS. LETT. Vol. 32, No. 4 (2015) 040201
the bifurcation details in the evolution of the dynamics.
Firstly, for the autonomous case, i.e., 𝐼 = 0, Eq. (1)
can be described as follows, where 𝑓 (𝑣) has the same
form
𝐶 𝑣˙ = 𝑓 (𝑣) − 𝐶𝛾 − 𝐶𝑤,
𝑤˙ = 𝑣 − 𝛾𝑓 (𝑣).
(2)
Two switching boundaries exist, denoted by Σ 1 : 𝑣 =
𝑎/2 and Σ2 : 𝑣 = (1 + 𝑎)/2. Furthermore, three
equilibria can be observed, which are expressed as
𝐸1 : {𝑣 = 0, 𝑤 = −𝛾} for 𝑣 < 𝑎/2, 𝐸2 : {𝑣 =
𝛾𝑎/(𝛾 − 1), 𝑤 = 𝛾𝑎/𝐶(𝛾 − 1) − 𝑎/𝛾 − 𝛾} for 𝑎/2 ≤ 𝑣 ≤
(1 + 𝑎)/2, 𝐸3 : {𝑣 = 𝛾/(𝛾 + 1), 𝑤 = 1/𝐶(𝛾 + 1) − 𝛾}
for 𝑣 > (1 + 𝑎)/2.
2
q/
M (0, 1.0)
q/⊲
1
Im (λ)
0.2
0.0
-0.2
-0.4
-0.6
-0.8
-1.0
Due to the fact that the vector field is non-smooth
while is still continuous, we can use the generalized
differential of Clarke to set up a generalized Jacobian matrix.[20] The generalized Jacobian of Eq. (2),
expressed by 𝐽G (𝑞) = {𝑞𝐽1(3) + (1 − 𝑞)𝐽2 , ∀𝑞 ∈ [0, 1]},
at the bifurcation point is the closed convex hull of
the Jacobians on each side of the switching boundary.
The eigenvalues of 𝐽G (𝑞), denoted by 𝜆G1,2 , are setvalued and form a path in the complex plane with 𝑞
as path parameter (see Fig. 2).
The path shown in Fig. 2 illustrates that the eigenvalues of 𝐽G (𝑞) are purely imaginary for 𝑞 = 0.5. The
path of the eigenvalues of 𝐽G (𝑞) shows that the discontinuous bifurcation point is a discontinuous Hopf
bifurcation,[16] and a new frequency, which is calculated at 𝜔𝑡1 = 1.0, is generated due to this bifurcation.
-0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0
1.0
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
10000
1.0
0.8
0.6
0.4
0.2
0.0
-0.2
-0.4
10100
q/
0
q/
-1
q/⊲
q/
-2
-0.3
-0.2
M ' (0, −1.0)
-0.1
0
0.1
0.2
0.3
Re (λ)
10200
10400
10150
10200
Fig. 2. Set of eigenvalues of the generalized Jacobian of
Eq. (2).
10600
10250
Fig. 1. Phase portrait and time series of the system.
The stability of 𝐸1 is determined by the eigenvalues of 𝐽1 , while the stabilities of 𝐸2 and 𝐸3 are determined by the eigenvalues of 𝐽2 and 𝐽3 , respectively,
where 𝐽1 and 𝐽2 , 𝐽3 are expressed with the above parameters, as follows:
⃒
⃒ ⃒
⃒
⃒ −1/𝐶 −1 ⃒ ⃒ −0.5 1 ⃒
⃒
⃒
⃒
⃒,
𝐽1 = 𝐽3 = ⃒
=
1+𝛾
0 ⃒ ⃒ 1.65 0 ⃒
with eigenvalues 𝜆1 = −0.25 + 1.26𝑖 and 𝜆2
1.26𝑖,
⃒
⃒ ⃒
⃒ 1/𝐶 −1 ⃒ ⃒ 0.5 −1
⃒
⃒=⃒
𝐽2 = ⃒
1 − 𝛾 0 ⃒ ⃒ 0.35 0
= −0.25 −
Now we turn to the non-autonomous case with
𝐼 ̸= 0, which implies that the dynamics of the system
will be influenced by the external stimulation. The
stimulation may cause the disappearance of the equilibria described above and drive the system to oscillate
periodically according to the frequency of the stimulation.
To investigate the bifurcation mechanism, we introduce the transformation as 𝜔0 𝑡 = 𝜃, and regard the
dimensionless time 𝑡 as a variable parameter, which
leads the vector field to the generalized autonomous
form described as Eq. (3), and 𝑓 (𝑣) here has the same
form that it had in Eq. (1),
𝐶 𝑣˙ = 𝑓 (𝑣) − 𝐶𝛾 − 𝐶𝑤 + cos(𝜃),
𝑤˙ = 𝑣 − 𝛾𝑓 (𝑣) − 𝛾 cos(𝜃).
Due to the linear structure of the system, though nonsmooth, we now assume the solution of Eq. (3) can be
written as
⃒
⃒
⃒,
⃒
with eigenvalues 𝜆′1 = 0.25 + 0.54𝑖 and 𝜆′2 = 0.25 −
0.54𝑖.
Obviously, the equilibria 𝐸1 and 𝐸3 are asymptotically stable focus points while 𝐸2 is an unstable
focus.
(3)
{(𝑣, 𝑤)|𝑣 = 𝑘11 sin 𝜃 + 𝑘12 cos 𝜃 + 𝑘10 ,
𝑤 = 𝑘21 sin 𝜃 + 𝑘22 cos 𝜃 + 𝑘20 }.
(4)
By substituting this solution (4) into Eq. (3), one may
obtain three types of coefficients in the solution corresponding to different conditions. The details can be
040201-2
CHIN. PHYS. LETT. Vol. 32, No. 4 (2015) 040201
described as Eq. (1) for 𝑎/2 ≤ 𝑣 ≤ (1 + 𝑎)/2,
Although the generalized autonomous system
Eq. (3) is actually a non-autonomous system, the value
of stimulation is constant with regard to any certain
time; that is, it could be looked as an autonomous system at any certain moment. The solution described in
Eq. (4) can then be regarded as a generalized equilibrium solution to Eq. (3). Due to the piecewise linearity
of the system, the stabilities of the related generalized
equilibrium solution at that moment should be discussed in the corresponding regions, respectively. For
regions I and III, the linearized matrix of the generalized equilibrium solution is expressed in the same
form as that of genuine equilibria 𝐸1 and 𝐸3 for the
autonomous case due to the piecewise linear structure
of the equations, which means the stabilities as well as
the corresponding dynamical properties of the generalized equilibrium solution are the same with those of 𝐸1
and 𝐸3 . That is, when the parameters are still fixed
at those values given above, the solution of Eq. (3) is a
stable focus for 𝑣 > (1+𝑎)/2 and 𝑣 < 𝑎/2 while an unstable focus for 𝑎/2 ≤ 𝑣 ≤ (1 + 𝑎)/2. The distribution
of these solutions could be seen in Fig. 3.
𝑘11 = 𝜔0 𝐶(𝜔02 − 1)/(𝐶 2 𝜔04 + 2𝐶 2 𝜔02 𝛾
− 2𝐶 2 𝜔02 + 𝐶 2 𝛾 2 − 2𝐶 2 𝛾 + 𝐶 2 + 𝜔02 ),
𝑘12 = − (𝐶 2 𝜔02 𝛾 + 𝐶 2 𝛾 2 − 𝐶 2 𝛾 + 𝜔02 )/(𝐶 2 𝜔04
+ 2𝐶 2 𝜔02 𝛾 − 2𝐶 2 𝜔02 + 𝐶 2 𝛾 2
− 2𝐶 2 𝛾 + 𝐶 2 + 𝜔02 ),
𝑘21 = − 𝜔0 (𝐶 2 𝜔02 𝛾 + 𝐶 2 𝛾 2 − 𝐶 2 𝛾 + 1)/(𝐶 2 𝜔04
+ 2𝐶 2 𝜔02 𝛾 − 2𝐶 2 𝜔02 + 𝐶 2 𝛾 2
− 2𝐶 2 𝛾 + 𝐶 2 + 𝜔02 ),
𝑘22 = 𝐶(𝜔02 𝛾 − 𝜔02 − 𝛾 + 1)/(𝐶 2 𝜔04 + 2𝐶 2 𝜔02 𝛾
𝑘10
− 2𝐶 2 𝜔02 + 𝐶 2 𝛾 2 − 2𝐶 2 𝛾 + 𝐶 2 + 𝜔02 ),
= 𝑎𝛾/(𝛾 − 1),
𝑘20 = − 𝑎(4𝐶𝛾 2 − 4𝐶𝛾 − 1)/𝐶(𝛾 − 1);
(2) for 𝑣 > (1 + 𝑎)/2,
𝑘11 = 𝜔0 𝐶(𝜔02 − 1)/(𝐶 2 𝜔04 − 2𝐶 2 𝜔02 𝛾
− 2𝐶 2 𝜔02 + 𝐶 2 𝛾 2 + 2𝐶 2 𝛾 + 𝐶 2 + 𝜔02 ),
0.2
0.0
-0.2
-0.4
-0.6
-0.8
-1.0
𝑘12 = − (𝐶 2 𝜔02 𝛾 − 𝐶 2 𝛾 2 − 𝐶 2 𝛾 − 𝜔02 )/(𝐶 2 𝜔04
− 2𝐶 2 𝜔02 𝛾 − 2𝐶 2 𝜔02 + 𝐶 2 𝛾 2
+ 2𝐶 2 𝛾 + 𝐶 2 + 𝜔02 ),
𝑘21 = − 𝜔0 (𝐶 2 𝜔02 𝛾 − 𝐶 2 𝛾 2 − 𝐶 2 𝛾 − 1)/(𝐶 2 𝜔04
− 2𝐶 2 𝜔02 𝛾 − 2𝐶 2 𝜔02 + 𝐶 2 𝛾 2
+ 2𝐶 2 𝛾 + 𝐶 2 + 𝜔02 ),
𝑘22 = − 𝐶(𝜔02 𝛾 + 𝜔02 − 𝛾 − 1)/(𝐶 2 𝜔04 − 2𝐶 2 𝜔02 𝛾
𝑘10
− 2𝐶 2 𝜔02 + 𝐶 2 𝛾 2 + 2𝐶 2 𝛾 + 𝐶 2 + 𝜔02 ),
= 𝛾/(𝛾 + 1),
𝑘20 = − (𝐶𝛾 2 + 𝐶𝛾 − 1)/𝐶(𝛾 + 1);
(I)
(II)
(III)
Stable
"solution"
Stable
"solution"
Unstable
"solution"
-0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0
Fig. 3. The generalized equilibrium solution to Eq. (3).
(3) for 𝑣 < 𝑎/2,
0.2
0.0
-0.2
-0.4
-0.6
-0.8
-1.0
𝑘11 = 𝜔0 𝐶(𝜔02 − 1)/(𝐶 2 𝜔04 − 2𝐶 2 𝜔02 𝛾
− 2𝐶 2 𝜔02 + 𝐶 2 𝛾 2 + 2𝐶 2 𝛾 + 𝐶 2 + 𝜔02 ),
𝑘12 = − (𝐶 2 𝜔02 𝛾 − 𝐶 2 𝛾 2 − 𝐶 2 𝛾 − 𝜔02 )/(𝐶 2 𝜔04
− 2𝐶 2 𝜔02 𝛾 − 2𝐶 2 𝜔02 + 𝐶 2 𝛾 2
+ 2𝐶 2 𝛾 + 𝐶 2 + 𝜔02 ),
𝑘21 = − 𝜔0 (𝐶 2 𝜔02 𝛾 − 𝐶 2 𝛾 2 − 𝐶 2 𝛾 − 1)/(𝐶 2 𝜔04
− 2𝐶 2 𝜔02 𝛾 − 2𝐶 2 𝜔02 + 𝐶 2 𝛾 2
+ 2𝐶 2 𝛾 + 𝐶 2 + 𝜔02 ),
𝑘22 = − 𝐶(𝜔02 𝛾 + 𝜔02 − 𝛾 − 1)/(𝐶 2 𝜔04 − 2𝐶 2 𝜔02 𝛾
𝑘10
𝑘20
− 2𝐶 2 𝜔02 + 𝐶 2 𝛾 2 + 2𝐶 2 𝛾 + 𝐶 2 + 𝜔02 ),
= 0,
= − 𝛾.
(I)
(II)
(III)
-0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0
Fig. 4. Overlap of the generalized equilibrium solution
and the real solution.
Note that the solution, which is a periodic function of
𝑡 with period 𝑇 = 2𝜋/𝜔0 , is not the real solution of
the system. The reason to introduce the form of the
solution here is that it may help to understand the
bifurcation mechanism of Eq. (3).
The phase plane is divided into three parts by the
switching boundary Σ1 and Σ2 . The generalized equilibrium solution is a stable focus in region I and III
while it is an unstable focus in region II.
We will now explain the mechanism of the periodic movement observed above. For this purpose, the
phase portrait shown in Fig. 1 is projected onto Fig. 3
040201-3
CHIN. PHYS. LETT. Vol. 32, No. 4 (2015) 040201
(seen as Fig. 4).
An arbitrary point on the trajectory could be set
as the starting point. For the case of the description,
we set it to locate at the switching boundary. On either side, the generalized equilibrium solution is an
unstable focus in region II while it is an asymptotically stable focus in region III. Thus, the trajectory
will tend to the stable attractor; namely, the solution
in region III with the angular frequency 𝜔𝑡2 = 1.26
is determined by the imaginary part of the complex
eigenvalues 𝜆1,2 , which is associated with the generalized equilibrium solution mentioned above. The anticlockwise direction of the trajectory can be demonstrated by the time series plotted in Fig. 1. The value
of 𝜔𝑡2 could also be proved by the numerical simulation; that is, the relevant frequency could be calculated at 𝜔𝑆2 = 2𝜋/𝑇2 ≈ 1.24 from Fig. 1. When the
trajectory is drawn into the stable generalized equilibrium solution in region III, it may move exactly along
the solution to the left with the stimulation angular frequency 𝜔0 = 0.05. The corresponding oscillation angular frequency obtained from simulation is
𝜔0𝑆 = 2𝜋/𝑇0 ≈ 0.050, which fits the theoretical value
perfectly. Note that the passage remaining on the solution is relatively slow, which can be called a rest
state (RS1 ). The trajectory moves along the stable
attractor until it meets the switching boundary.
Non-smooth bifurcation, i.e., the discontinuous
Hopf bifurcation, occurs when the trajectory crosses
the switching boundary Σ2 at the point P1 . The solution in region III then loses its stability via this bifurcation and the trajectory circles with the angular
frequency 𝜔𝑡1 = 1.0 determined by the pair of pure
imaginary zeros of 𝐽G (𝑞). This frequency can also be
verified by the above numerical simulation. Based on
time series shown in Fig. 1, the associated angular frequency can be computed as 𝜔𝑆1 = 2𝜋/𝑇1 ≈ 0.971. It
is very asymptotic to the value of 𝜔𝑡1 generated from
discontinuous Hopf bifurcation. Since the generalized
equilibrium solution in region II is unstable, the trajectory will quickly pass region II and tend to the stable
solution located in region I with the angular frequency
𝜔𝑡2 = 1.26 determined by the property of the generalized equilibrium solution in region I, which is the same
as that of region III. On account of this the values of
both 𝜔𝑡1 and 𝜔𝑡2 are much higher than that of 𝜔0 , the
passage related to 𝜔𝑡1 and 𝜔𝑡2 is relatively fast, which
is called spiking (SP1 ).
Until the trajectory converges to the stable generalized equilibrium solution in region I, it moves exactly along the solution, which is again similar to that
described above. It is a state that is relatively slow,
which can be called the rest state (RS2 ).
Until the trajectory runs across the switching
boundary Σ1 , discontinuous Hopf bifurcation occurring at P2 leads the instability of RS2 , resulting in
another fast passage; i.e., spiking (SP2 ). Due to the
similarity of the characteristics of generalized equilibrium solution between regions I and III, the behavior
is similar to that described above, until the trajectory
reaches the point P0 to start another period.
In summary, bursting with coexistence of three frequencies is observed for the piece-wise linear planar
neuron model with periodic stimulation. Based on the
characteristics of this non-smooth non-autonomous
system, the concept of generalized autonomous system as well as generalized equilibrium solution is introduced. Furthermore, the appropriate form of the
generalized equilibrium solution is presented to elaborate the non-smooth bifurcation mechanism of the
system, as well as the origin of the respective frequencies of oscillation. It is suggested that discontinuous
Hopf bifurcation occurring at the non-smooth boundaries leads to the behavior of the system exhibiting a
transition between a rest state and a spiking state by
generating a frequency which is much higher than that
of the stimulation, meanwhile the frequency of spiking state may change slightly due to the dynamical
property of the generalized equilibrium solution.
References
[1] Johnson S W, Seutin V and North R A 1992 Science 258
665
[2] Ivanchenko M V, Osipov G V, Shalfeev V D and Kurths J
2007 Phys. Rev. Lett. 98 108101
[3] Ji Y and Bi Q S 2011 Chin. Phys. Lett. 28 090201
[4] Wang H X, Wang Q Y and Zheng Y H 2014 Sci. Chin.
Technol. Sci. 57 872
[5] Gu H G and Chen S G 2014 Sci. Chin. Technol. Sci. 57
864
[6] Hodgkin A L and Huxley A F 1952 J. Physiol. 117 500
[7] Chay T R, Fan Y S and Lee Y S 1995 Int. J. Bifurcation
Chaos Appl. Sci. Eng. 5 595
[8] Rinzel J 1985 Ordinary Partial Differential Equations
(Berlin: Springer-Verlag)
[9] Sherman A and Rinzel J 1992 Proc. Natl. Acad. Sci. USA
89 2471
[10] Tiesinga P H E 2002 Phys. Rev. E 65 041913
[11] Coombes S 2008 SIAM J. Appl. Dyn. Syst. 7 1101
[12] Fitzhugh R 1961 Biophys. J. 1 445
[13] Yamashita Y and Torikai H 2014 IEEE T. Circuits-II: Express Briefs 61 54
[14] Izhikevich E M 2000 Int. J. Bifurcation Chaos Appl. Sci.
Eng. 10 1171
[15] Yang Z Q and Lu Q S 2008 Sci. Chin. Phys. Mech. Astron.
51 687
[16] Zhang F, Lubbe A, Lu Q S and Su J Z 2014 Discrete Continuous Dynamical Syst. Ser. S 7 1363
[17] Xie Y, Kang Y M, Liu Y and Wu Y 2014 Sci. Chin. Technol.
Sci. 57 914
[18] Tonnelier A 2003 SIAM J. Appl. Math. 63 459
[19] Jaume L, Manuel O, Manuel O and Enrique P 2013 Nonlinear Anal.: Real World Appl. 14 2002
[20] Leine R I 2006 Physica D 223 121
040201-4