國 立 交 通 大 學
應用數學所
碩
士 論 文
神經網路上的同步方程之
多個穩定狀況的分析
Multistate Stability of Synchronous Equations
in Hindmarsh-Rose Networks
研
究 生： 黃 于 哲
指導教授：
莊 重 教授
華 民 國 一 ○ 三 年 七 月
神經網路上的同步方程之
多個穩定狀況的分析
Multistate Stability of Synchronous Equations
in Hindmarsh-Rose Networks
研
究 生：黃 于 哲 Student：Yu-Jhe Huang
指導教授：
莊 重 Advisor： Jonq Juang
國
立 交 通 大 學
應
用 數 學 所
碩
士 論 文
A Thesis
Submitted to Department of Applied Mathematics
College of Science
National Chiao Tung University
In Partial Fulfillment of the Requirements
for the Degree of
Master
in
Applied Mathematics
July 2014
Hsinchu, Taiwan, Republic of China
神經網路上的同步方程之
多個穩定狀況的分析
學生：黃 于 哲
指導教授：莊 重
國立交通大學應用數學學系﹙研究所﹚碩士班
摘
要
在這篇論文中，我們利用幾何奇異擾動理論研究來自
Hindmarsh-Rose 網路之同步化方程的多重穩定狀態。
我們的主要
結果如下
: 首先，我們給出同步化 Hindmarsh-Rose 方程多重穩定狀
態的解釋，例如我們能下結論說在爆裂
(bursting)解與具有 canard 現
象的週期解能共存。
其次，我們可以充分了解從初始狀態至穩定狀
態的過程。
最後，我們可識別出穩定狀態的吸引範圍。 這些都說
明了用幾何奇異擾動理論理解實際生物系統的全域動態性質是相當
有用的。
關鍵字
: 多重穩定狀態、同步化方程、幾何奇異擾動理論。
iMultistate Stability of Synchronous Equations
in Hindmarsh-Rose Networks
Student: Yu-Jhe Huang
Advisors: Jonq Juang
Department of Applied Mathematics
National Chiao Tung University
ABSTRACT
In this thesis, geometric singular perturbation theory is applied to
investigate multistate stability of synchronous equations derived from
Hindmarsh-Rose Networks. Our main results contain the following. First,
explanation of multistability of the synchronous Hindmarsh-Rose equation
can be given. For instance, we are able to conclude among other things that
a bursting solution and a periodic solution with canard explosion can
coexistence. The transition from initial states toward stable states can be
fully predicted. Finally, the attraction region with respect to each stable
state can be identified. This illustrates the power of using singular
perturbation theory to understand the global dynamical properties of
realistic biological systems.
Key words: multistate, synchronous equations, geometric singular
perturbation theory.
誌 謝
這篇論文的完成，仰賴許多人的幫助與支持，其中
以指導老師莊重教授及梁育豪學長對我這個工作的完
成幫助最多，莊重老師對我的指導使我對此工作有明
確的目標及想法，而梁育豪學長則是幫助我論文中數
值相關的部分和一些應注意的細節。 還有我必須感
謝同為莊老師學生的黃俊銘學長、林辰燁學長以及謝
瑞洸，他們平常或多或少對我的幫助及鼓勵。
在我研究所兩年中，我也必須感謝下列的一些人，
同屆的研究生郭柏均、陳學民、江培華、周敬淯以及
同間研究室的黃建順、蔡志奇、蔡睿翊，他們對我在
研究所的生活有很大的影響。
最後我在此感謝我的父母，他們讓我可以無後顧之
憂的作研究。
黃于哲
西元 2014 年 7 月
iiiTable of Contents
Chinese Abstract
………
i
English Abstract
………
ii
Acknowledgement
………
iii
Table of Contents
………
iv
1.
Introduction………
1
2.
Preliminaries………
2
2.1
Geometric Singular Perturbation Theory
2
2.1.1
Basic Setting………
2
2.1.2
Fenichel’s Theorems………
3
2.2
Hopf Bifurcation………
7
2.2.1
Basic assumptions………
8
2.2.2
Example………
8
2.2.3
Center manifold reduction………
9
2.2.4
Preparation for statement of the
Poincaré–Andronov–Hopf Bifurcation
theorem………
9
2.2.5
ThePoincaré–Andronov–Hopf
Bifurcation theorem………
9
3.
Synchronous Equations in
Hindmarsh-Rose Networks………
10
4.
Dynamics of Synchronous Equations…
11
4.1
Critical manifolds………
12
4.2
ynamics on the fast subsystem (17a)……
16
4.2.1
𝑘𝑘𝑘𝑘
𝑠𝑠= 0.8………
16
4.2.2
𝑘𝑘𝑘𝑘
𝑠𝑠= 0.81………
19
4.2.3
𝑘𝑘𝑘𝑘
𝑠𝑠= 0.83………
22
4.2.4
𝑘𝑘𝑘𝑘
𝑠𝑠= 0.88………
25
4.3
Numerical results for various of
𝑘𝑘𝑘𝑘
𝑠𝑠…
28
4.4
Analysis and conclusions………
33
References
………
34
1 Introduction
The fundamental building block of every nervous system is the neuron. The be-havior of even a single cell can be quite complicated [1], [2] and [3]. For instance, an individual single cell may fire repetitive action potentials or bursts of action potentials that are followed by a silent phase of near quiescent behavior [4]. In the last decades, many biological neuron models have been proposed for an accurate description and prediction of biological phenomena. The pioneering work in such direction is due to Hodgkin and Huxley. To simplify such a model, simpler approximations, namely, the second order systems such as the FitzHugh-Nagumo (FN) and Morris-Lecar neuron models have been proposed. However, the second order models are not able to re-produce some interesting phenomena such as terminating themselves by triggering a set of stable firings. Hence, the Hindmarsh-Rose (HR) model was added a third dynamical component, whose role is to tune the above subsystem over the mono-and bistability regions in order to activate or terminate the neuronal response. The third order system of HR has turned out to be accurate in capturing both qualitative and quantitative aspects of experimental data [5]–[8]. Furthermore, major neuronal behaviors such as spiking, bursting, and chaotic regime have been produced by such HR model [1], [8] and [9].
Example of population rhythms include synchronous behavior, in which easy all in the network fires at the same time, and clustering [10], [11], [12], [13] and [14], in which the entire population of cells breaks up into subpopulations or blocks; every cell within a single block fires synchronously and different blocks are desynchronyed from each other. Neural synchronization has been suggested as particularly relevant for neuronal signal transmission and coding in the brain. Brain [15]–[22] oscilla-tions that are ubiquitous phenomena in all brain areas eventually get into synchrony and consequently allow the brain to process various tasks from cognitive to motor tasks. Indeed, it is hypothesized that synchronous brain activity is the most likely mechanism for many cognitive functions such as attention, feature binding, learning, development, and memory function.
In [23], a new phenomena of the multistate and multistage synchronization of HR neurons with excitatory chemical and electrical synapses over the complex network are analytical studied. In this thesis, we are to use geometric singular perturbation theory to understand such phenomena of multistate synchronization. Specifically, we shall consider the HR networks restricted on its synchronous manifold. The re-sulting equation is to be termed synchronous equation with respect to HR networks. Our main results contain the following. First, explanation of multistability of the synchronous Hindmarsh-Rose equation can be given. For instance, we are able to conclude among other things that a bursting solution and a periodic solution with canard explosion can coexistence. The transition from initial states toward stable states can be fully predicted. Finally, the attraction region with respect to each
sta-ble state can be identified. This illustrates the power of using singular perturbation theory to understand the global dynamical properties of realistic biological systems.
The thesis is organized as follows. In Section 2, we state some basic facts concern-ing geometric sconcern-ingular perturbation theory and Hopf bifurcation. The synchronous equation of HR networks and some known properties, observed by computer simula-tions done in [23], are stated in Section 3. The main results of the thesis is contained in Section 4.
2 Preliminaries
For ease of reference, we introduce the well-known geometric singular perturbation theory. For more details we refer to [24] and [25].
2.1 Geometric Singular Perturbation Theory
2.1.1 Basic Setting
The equations under consideration are of the form {
x′ = f (x, y, ε),
y′ = εg(x, y, ε), (1)
where ′ = d
dt, x ∈ R
n, y ∈ Rm and ε ∈ R is a small parameter (0 < ε ≪ 1) which
gives the equations a singular character.
The following statement is one of the basic assumptions:
Hypothesis 1. The functions f and g are both assumed to be Cr on a set U × I,
where U ⊂ Rn+m is open and I is an open interval, containing 0 and r≥ 1.
From (1), let s = εt, then equation (1) is transformed to the following form {
ε ˙x = f (x, y, ε),
˙
y = g(x, y, ε), (2)
where ˙ = d
ds. The time scale given by s is said to be slow whereas that for t is fast.
In fact, as long as ε ̸= 0 the two systems are equivalent. Thus we call (1) the fast
system and (2) the slow system.
Letting ε−→ 0, that the limiting systems of (1) and (2) are, respectively, {
x′ = f (x, y, 0),
and {
f (x, y, 0) = 0,
˙
y = g(x, y, 0). (4)
Let N ={ (x, y) ∈ U | f(x, y, 0) = 0}. In view of system (3), each point of N is a fixed point and the system (3) has nontrivial dynamical behaviors on U\N. On the other hand, system (4), defined only on N , usually has nontrivial dynamical behaviors on
N . We next define the concept of a set being normally hyperbolic
Definition 1. Let M0 ⊂ N. The set M0 is said to be normally hyperbolic if
Dxf (˜x, ˜y, 0) has exactly n eigenvalues λ with Re(λ)̸= 0 for any (˜x, ˜y) ∈ M0.
Hypothesis 2. Let M0 ⊂ N. The set M0 is assumed to be a compact manifold and
is normally hyperbolic.
The set M0 will be called the critical manifold.
Since from Hypotheses 2, we have that Dxf (˜x, ˜y, 0) is invertible for any (˜x, ˜y) ∈
M0. By the Implicit Function Theorem, x can locally be solved for y, that is, x =
h0(y) for some smooth function h0. Since the set which is a graph of some smooth
function is easier to manipulate, we want to decompose M0 into a finite number of
its subset that is a graph of some smooth function and satisfies Hypotheses 2. (i.e.
M0 = M0(1) ∪ M (2) 0 ∪ · · · ∪ M (k) 0 where M (i)
0 is a graph of some smooth function and
satisfies Hypotheses 2 for i=1,2,. . . , k.
In order to simplify discussion, we consider that M0 is a graph of function satisfying
the following two hypotheses.
Hypothesis 3. M0 is normally hyperbolic and there are a compact subset K ofRm
and a function h0 defined on K such that M
0 ={ (x, y) ∈ U | x = h0(y)}.
Hypothesis 4. h0 is a smooth function on K and K is simply connected domain
whose boundary is an (m− 1)-dimensional C∞ submanifold.
In addition, we also assume that the following hypothesis holds.
Hypothesis 5. Dxf (˜x, ˜y, 0) has ℓ eigenvalues λ with Re(λ) < 0 and k(= n − ℓ)
eigenvalues λ with Re(λ) > 0 for any (˜x, ˜y) ∈ M0.
2.1.2 Fenichel’s Theorems
We give some definitions and some notations in order to state Fenichel’s Theorems. Notation. The notation x· t is used to denote the application of a flow after time t to the initial condition x. Similarly, V · t denote the application of a flow after time
t to a set V , and x· [t1, t2] is the resulting trajectory if the flow is applied over the
Definition 2. A set M is locally invariant under the flow from (1) if there exists a neighbourhood V of M so that for all x∈ M, x · [0, t] ⊂ V implies that x · [0, t] ⊂ M, similarly with [0, t] replaced by [t, 0] when t < 0.
Theorem 1. (Fenichel’s First Theorem, see e.g., [24])
Under the Hypotheses 1 and 2. Then for ε > 0 and sufficiently small, there exists a manifold Mε that lies within O(ε) of M0 and is diffeomorphic to M0. Moreover
it is locally invariant under the flow of system (1), and Cr, including in ε, for any
r < +∞.
The manifold Mε will be called the slow manifold.
Theorem 2. (Fenichel’s First Theorem for graph version, [24])
Under the Hypotheses 1, 3and 4. Then for ε > 0 and sufficiently small, there exists a function hε, defined on K, so that M
ε = { (x, y) | x = hε(y), y ∈ K} is locally
invariant under the flow of system (1). Moreover hε is Cr, for any r < +∞, jointly
in y and ε.
Remark. The diffeomorphism between Mε and M0 follows easily in this formulation
through the diffeomorphism of the graph to K. Flow on Mε:
y′ = εg(hε(y), y, ε). (5)
In the alternative slow scaling we recast (5) as ˙
y = g(hε(y), y, ε). (6)
which has distinct advantage that a limit exists as ε → 0, given by g(h0(y), y, 0)
which naturally describes a flow on the critical manifold M0, and is exactly the second
equation in (2). Using this theorem and this resulting equation (6), the problem of studying (1), at least on Mε is reduced to a regular perturbation problem.
Example. Consider the system ˙u1 = u2 ˙u2 = u3 ε ˙u3 = u4 ε ˙u4 = u5 ε ˙u5 = u6
ε ˙u6 =−Au5− u3− cu2− f(u1),
where A > 0, f (ϕ) = ϕ(ϕ− a)(1 − ϕ), a < 1 2, c: constant. Let x = u3 u4 u5 u6 , y = [ u1 u2 ]
The critical manifold M0 can be taken as any compact subset of {u4 = u5 = u6 =
0, u3 =−cu2− f(u1)} and will shall be large enough to contain any of the dynamics
of interest.
The eigenvalues of linearization at any point of M0, other than the double eigenvalue
at 0 are seen to be solutions of quartic equation µ4+ Aµ2+ 1 = 0 which are not pure
imaginary if 0 < A < 2.
The equation for the slow flow on the critical manifold M0are given by
{
˙u1 = u2,
˙u2 =−cu2− f(u1).
The slow manifold Mε, which exists by virtue of theorem 1,is given by the equations
(u3, u4, u5, u6) = hε(u1, u2) = (−cu2 − f(u1), 0, 0, 0) + O(ε)
and the equations on Mε are
{
˙u1 = u2,
˙u2 =−cu2− f(u1) + O(ε).
We know nothing of the flow off the slow manifold and this must now be addressed. The slow manifold possesses accompanying stable and unstable manifolds that are perturbations of the corresponding manifolds when ε = 0.
Theorem 3. (Fenichel’s Second Theorem, see e.g., [25])
Under the Hypotheses 1 and 2. Then for ε > 0 and sufficiently small, there exists manifolds Ws(M
ε) and Wu(Mε) that are O(ε) close and are diffeomorphic to Ws(M0)
and Wu(M
0) respectively, and they are each locally invariant under the flow of system
(1), and Cr, including in ε, for any r <∞.
We want to state Fenichel’s Second Theorem for graph version but we need a detailed analysis and more notations. So we postpone this after Fenichel’s third Theorem. The following is talking about a motivation for Fenichel’s third Theorem. First, we observe that
Ws(M0) = ∪ v0∈M0 Ws(v0), Wu(M0) = ∪ v0∈M0 Wu(v0)
In other words, the manifolds Ws(v
0) and Wu(v0) form collections of fibers for
Ws(M
0) and Wu(M0), respectively, with base points v0 ∈ M0.
Second, a natural question to ask now, is whether the individual stable and unstable manifolds Ws(v
0) and Wu(v0) perturb to analogous objects ? —— Fenichel’s third
Theorem answers this question. In order to avoid difficulties we restrict ourselves to a neighbourhood D of Mεin which the linear terms of (1) are dominant, and consider
only trajectories in Wu(M
ε) that have not left D in forward time (yet), and
trajecto-ries in Ws(M
ε) that have not left D in backward time. To facilitate this discussion,
we need a definition.
Definition 3. The forward evolution of a set V ⊂ D restricted to D is given by the
set
Figure 1: A diagram for Fenichel’s Third Theorem. Now, we sate Fenichel’s Third Theorem.
Theorem 4. [Fenichel’s Third Theorem]1
Under the Hypotheses 1, 2 and 5. Then if ε > 0 and sufficiently small then for every vε ∈ Mε there are an ℓ-dimensional manifold Ws(vε) ⊂ Ws(Mε) and an
k-dimensional manifold Wu(v
ε) ⊂ Wu(Mε), that are O(ε) close and diffeomorphic to
Ws(v
0) and Wu(v0) respectively.
Moreover, they are Cr for any r, including in v and ε.
The families { Ws(v
ε)| vε ∈ Mε} and { Wu(vε)| vε∈ Mε} are invariant in the sense
that
Ws(vε)·Dt ⊂ Ws(vε· t)
if vε· s ∈ D for all s ∈ [0, t] (t > 0) and
Wu(vε)·Dt ⊂ Wu(vε· t)
if vε· s ∈ D for all s ∈ [t, 0] (t < 0).
The idea is, that the order in which the flow after time t is applied to a base point and the fiber of a base point is constructed does not matter, as depicted in this figure (1).
In the unperturbed setting (1) with ε = 0, the decay in forward time of points in
Ws(M
0) to M0 is clearly to the base point v0 of their fiber, where the decay rate as
t→ ∞ is exponential, since all associated eigenvalues have nonzero real part.
The fibers give a very useful matching between the points in Ws(M
ε) and partners
they have in Mε. One can then see that the decay of points in Ws(Mε) to Mε is
actually to the base point of the fiber, this gives a decay result with ”asymptotic phase”; similarly for points in Wu(M
ε).
Preparation for statement of Fenichel’s Second Theorem for graph version Suppose that M0satisfies the Hypotheses1,3and4. Without loss of generality, we
can assume that h0(y) = 0 for all y∈ K. Indeed, we can replace x by ˆx = x − h0(y)
and recompute the equations. For each y ∈ K, there are subspaces S(y) and U(y), corresponding, respectively, to stable and unstable eigenvalues. Since the eigenvalues are bounded uniformly away from the imaginary axis over K, the dimensions of S(y) and U (y) are independent of y. Let dim S(y) = ℓ and dim U (y) = k. Since K is simply connected by hypothesis4, we can smoothly choose bases for S(y) and U (y). Changing the coordinates to in terms of these new bases, we can set x = (a, b), where a∈ Rℓ and b ∈ Rk, so that our equations have the form
a′ = A(y)a + F1(x, y, ε), b′ = B(y)b + F2(x, y, ε), y′ = εg(x, y, ε), (7)
where σ(A(y)) ⊂ {λ | Re(λ) < 0} and σ(B(y)) ⊂ {λ | Re(λ) > 0}. Both F1 and
F2 consist of any higher order terms in x, y for each ε; to be precise, we have the
estimates |Fi| ≤ γ(|x| + ε), i = 1, 2 and γ can be taken to be as small as desired by
restricting to a set with|a| and |b| small.
With this notation established, we can determine Ws(M
ε) and Wu(Mε) as graphs
and give the following restatement of theorem 3.
Theorem 5 (Fenichel’s Second Theorem for graph version). 2
If ε > 0, but sufficiently small, then ,for some ∆ > 0,
(a) there is a function a = hs(b, y, ε) defined for y ∈ K and |b| ≤ ∆, so that
the graph Ws(M
ε) = {(a, b, y) | a = hs(b, y, ε)} is locally invariant under (7).
Moreover, hs(b, y, ε) is Cr in (b, y, ε) for any r < +∞.
(b) there is a function b = hu(a, y, ε) defined for y ∈ K and |b| ≤ ∆, so that
the graph Wu(M
ε) = {(a, b, y) | b = hu(a, y, ε)} is locally invariant under (7).
Moreover, hu(a, y, ε) is Cr in (a, y, ε) for any r < +∞.
These theorems also apply when ε = 0 and provide the stable and unstable man-ifolds of the known critical manifold, the existence of which is also guaranteed by the usual stable and unstable manifold theorems at critical points. These latter two theorems then assert that these manifolds, Ws(M
ε) and Wu(Mε), are perturbed from
Ws(M
0) and Wu(M0) respectively. Theorem 1 can be concluded from Theorem 3 by
taking the intersection of Ws(M
ε) with Wu(Mε). Locally, the Implicit Function
The-orem gives the intersection as a graph, and these functions can be patched together since K is a compact set. Moreover, we need only give the stable manifold, as that of the unstable manifold follows immediately by a reversal of time.
2.2 Hopf Bifurcation
In this subsection, we review Hopf Bifurcation. For more details we refer to [26].
2.2.1 Basic assumptions
First, we consider a continuous-time system depending on a parameter.
˙x = f (x, µ), x∈ Rn, µ∈ R, (8)
where f is smooth with respect to both x and µ. If x = x0be a hyperbolic equilibrium
in the system for µ = µ0then under a small parameter variation the equilibrium moves
slightly but remains hyperbolic. In general, there are only two ways in which the hyperbolicity condition can be violated. Either a simple real eigenvalue approaches zero and we have λ1 = 0, or a pair of simple complex eigenvalues reaches the imaginary
axis and we have λ1,2 =±iω0, ω0 > 0 for some value of the parameter.
Definition 4. We say that a fixed point of the system (8) has the Hopf (or Poincaré–
Andronov–Hopf ) bifurcation at µ = µ0 if a fixed point of the system (8) at µ = µ0 has
λ1,2 =±iω0, ω0 > 0 and a small perturbation at µ = µ0 results in stable or unstable
limit cycle from the fixed point.
We suppose that the system (8) satisfy the following assumptions.
Assumption 1. Assume that (0, 0)∈ Rn×R is the fixed point of (8), i.e., f (0, 0) = 0.
Assumption 2. Assume that Dxf (0, 0) has two purely imaginary eigenvalues with
the remaining n− 2 eigenvalues having nonzero real parts. 2.2.2 Example
Example. Consider the following planar system depending on one parameter µ: {
˙x =−y + µx − x(x2+ y2), ˙
y = x + µy− y(x2+ y2), (9) This system has the equilibrium (0, 0) for all µ with the Jacobian matrix
A =
[
µ −1
1 µ
]
having eigenvalues λ1,2 = µ± i. For µ = 0, the linearized equations of system (9)
have a center, and as µ varies from negative to positive, the fixed point changes from a stable focus to an unstable focus. Thus the system (9) has a Hopf bifurcation at
µ = 0.
In polar coordinates, {
˙r = µr− r3,
˙
θ = 1.
We see that at µ = 0 the system has a stable focus at origin and for µ > 0 the system has a stable limit cycle
γµ(t) =
√
2.2.3 Center manifold reduction
Under assumptions 1 and 2, by the center manifold theorem, we know that the orbit structure near (x, µ) = (0, 0) is determined by the vector field (8) restricted to the center manifold. This restriction gives us a one-parameter family of vector fields on a two-dimensional center manifold. On the center manifold the vector field (8) has the following form
[ ˙x1 ˙x2 ] = [ Re λ(µ) −Im λ(µ) Im λ(µ) Re λ(µ) ] [ x1 x2 ] + [ g1(x1, x2, µ) g2(x1, x2, µ) ] , (10)
(x1, x2, µ)∈ R×R×R and where g1 and g2 are nonlinear in x1 and x2 and λ(µ), λ(µ)
are the eigenvalues of the vector field linearized about the fixed point at the origin. Let λ(µ) = α(µ) + iω(µ) then, by assumption, we have α(0) = 0 and ω(0)̸= 0. 2.2.4 Preparation for statement of the Poincaré–Andronov–Hopf
Bifur-cation theorem
We transform (10) into normal form. The normal form was found to be ˙x = α(µ)x− ω(µ)y + (a(µ)x − b(µ)y)(x2+ y2) +O(|x|5,|y|5), ˙
y = ω(µ)x + α(µ)y + (b(µ)x + a(µ)y)(x2+ y2) +O(|x|5,|y|5). (11) In polar coordinates (11) is given by
˙r = α(µ)r + a(µ)r3+O(r5), ˙
θ = ω(µ) + b(µ)r2+O(r4). (12) Because we are interested in the dynamics near µ = 0, it is natural to Taylor expand the coefficients in (12) about µ = 0. Equation (12) becomes
˙r = α′(0)µr + a(0)r3+O(µ2r, µr3, r5), ˙
θ = ω(0) + ω′(0)µ + b(0)r2+O(µ2, µr2, r4). (13) where ′ denotes differentiation with respect to µ and we have used the fact that
α(0) = 0.
2.2.5 The Poincaré–Andronov–Hopf Bifurcation theorem
Theorem 6. [The Poincaré–Andronov–Hopf Bifurcation Theorem] 3
Consider the full normal form (13).
Then, for µ sufficiently small, Case 1, Case 2, Case 3, and Case 4 described below hold.
Case 1: α′(0) > 0, a(0) > 0.
In this case the origin is an unstable fixed point for µ > 0 and an asymptotically stable fixed point for µ < 0, with an unstable periodic orbit for µ < 0.
Case 2: α′(0) > 0, a(0) < 0.
In this case the origin is an asymptotically stable fixed point for µ < 0 and an unstable fixed point for µ > 0, with an asymptotically stable periodic orbit for µ > 0.
Case 3: α′(0) < 0, a(0) > 0.
In this case the origin is an unstable fixed point for µ < 0 and an asymptotically stable fixed point for µ > 0, with an unstable periodic orbit for µ > 0.
Case 4: α′(0) < 0, a(0) < 0.
In this case the origin is an asymptotically stable fixed point for µ < 0 and an unstable fixed point for µ > 0, with an asymptotically stable periodic orbit for µ < 0.
3 Synchronous Equations in Hindmarsh-Rose
Net-works
To derive synchronous equations in Hindmarsh-Rose Networks, we begin with the introduction of the Hindmarsh-Rose equation which governs the dynamics of single neuron. The equation has the following form.
˙x = f (x) + y− z + q, ˙ y =−y − 5x2+ 1, ˙z = µ(b(x− x0)− z), (14) where f (x) = ax2 − x3 and a, b, q, x
0 and µ are parameters. Moreover, µ is small
parameter (0 < µ ≪ 1). We are now in a position to consider a network of excita-tory HR neurons with bidirectional electrical coupling and unidirectional excitaexcita-tory chemical coupling. The equations of motion are the following. For i = 1, . . . , n,
˙xi = f (xi) + yi− zi+ q + σ n ∑ j=1 gijxj − gs(xi− v) n ∑ j=1 cijp(xj), ˙ yi =−yi− 5x2i + 1, ˙zi = µ(b(xi− x0)− zi), (15) where G =: (gij), n ∑ j=1 gij = 0 for all i, C =: (cij), cij = 0 or 1, cii= 0, n ∑ j=1 cij = k for all i, D =: (dij), and dij = { −k if i = j, cijif i ̸= j,
1
1 + e−λ(xj−θs) is the chemically synaptic coupling function, k represents the number
of chemical signals each neuron receives, and σ is the coupling strength for electrical synapses via gap junctions.
We next describe the synchronous equation of HR network (15). On the synchronous manifold { (x1, y1, z1) = (x2, y2, z2) = · · · = (xn, yn, zn) = (x, y, z)}, dynamics of HR
network is governed by the following equations: ˙x = f (x) + y− z + q − kgs(x− v)p(x) ˙ y =−y − 5x2+ 1, ˙z = µ(b(x− x0)− z), (16) where f (x) = ax2−x3, p(x) = 1
1+e−λ(x−θs) and a, b, q, x0, µ, v, λ and θsare parameters.
Moreover, µ is small parameter (0 < µ≪ 1).
In [23], it was reported numerically that synchronous equation (16) generate some various mutisatable states as the parameter kgs varies. Specifically, with the choice
of the parameters a = 2.6, b = 4, q = 4, x0 = −1.6, v = 4, λ = 10, θs = −0.25
and the small parameter µ = 0.01, they are able to observe the following mutistable states with various range of kgs, which are summarized in Table1.
Table 1: The dynamics of synchronous equation (16) with various range of kgs. The
multi-stability of (16) is observed with kgs∈ [0.809, 0.85].
kgs kgs< 0.808 0.809≤ kgs≤ 0.813 0.814≤ kgs≤ 0.85 kgs≥ 0.87
Stable regular bursting Stable regular bursting
Types Stable regular bursting Stable steady state
Stable periodic solution Stable steady state
4 Dynamics of Synchronous Equations
In view of equation of (1), we first note that the corresponding fast and slow variables in (16) are, respectively, x = (x, y)∈ R2 and y = z ∈ R. Consequently, the
fast limiting system of (16) has the following form ˙x = 2.6x2 − x3 + y− z + 4 − kgs(x− 2) 1 1 + e−10(x+0.25), ˙ y =−y − 5x2+ 1, (17a) ˙z = 0, (17b)
where z is regarded as a parameter. To apply geometric singular perturbation theory to equation (16), we shall construct a suitable critical manifold as well as the local and global geometric properties of the fast limiting system (17).
4.1 Critical manifolds
In this subsection, we shall choose critical manifolds with respect to equation (16). Let Nkgs ={(x, y, z) ∈ R 3 | 2.6x2−x3+y−z+4−kg s(x−2) 1 1 + e−10(x+0.25) = 0 and−y−5x 2+1 = 0}.
For each fixed z, Nkgs contains the set of the equilibriums of the fast subsystem Eq.
(17a). The geometric singular perturbation theory needs to define a critical manifold from Nkgs which satisfies Hypotheses 1, 2and 5of Eq. (16).
Let (x, y, z)∈ Nkgs then y = 1− 5x2 z = 2.6x2− x3+ (1− 5x2) + 4− kgs(x− 2) 1 1 + e−10(x+0.25) =−2.4x2− x3+ 5− kgs(x− 2) 1 1 + e−10(x+0.25) Let (x0, y0, z0)∈ Nkgs and let ˜x = x− x0, ˜y = y− y0, ˜z = z− z0.
˙˜x = ˙x = 2.6x2− x3+ y− z + 4 − kg s(x− 2) 1 1 + e−10(x+0.25) = 2.6(˜x + x0)2− (˜x + x0)3 + (˜y + y0)− (˜z + z0) + 4− kgs((˜x + x0)− 2) 1 1 + e−10((˜x+x0)+0.25) = 2.6(˜x2+ 2x0x + x˜ 20)− (˜x 3+ 3x2 0x + 3x˜ 0x˜2+ x30) + (˜y + y0)− (˜z + z0) + 4− kgs(˜x + x0− 2) 1 1 + e−10((˜x+x0)+0.25) =−˜x3+ (2.6− 3x0)˜x2+ (5.2x0− 3x20)˜x + ˜y− ˜z + (2.6x 2 0− x 3 0+ y0− z0 + 4) − kgs(˜x + x0− 2) 1 1 + e−10((˜x+x0)+0.25)
˙˜y = ˙y =−y − 5x2+ 1 =−(˜y + y
0)− 5(˜x + x0)2+ 1
=−˜y − 5˜x2− 10x0x + (˜ −y0− 5x20+ 1) =−˜y − 5˜x
2− 10x 0x˜
˙˜z = ˙z = 0.
To determine the stability and bifurcation points of a critical manifold which is a subset of Nkgs, consider the linearized equation of above system
˙˜x ˙˜y ˙˜z = 5.2x0− 3x20− kgs(1+e−10(x0+0.25))+10kgs(x0−2)e−10(x0+0.25) (1+e−10(x0+0.25))2 1 −1 −10x0 −1 0 0 0 0 ˜ x ˜ y ˜ z
det 5.2x0− 3x20− kgs(1+e−10(x0+0.25))+10kgs(x0−2)e−10(x0+0.25) (1+e−10(x0+0.25))2 − λ 1 −1 −10x0 −1 − λ 0 0 0 −λ =−λ det 5.2x0− 3x20− kgs(1+e−10(x0+0.25))+10kgs(x0−2)e−10(x0+0.25) (1+e−10(x0+0.25))2 − λ 1 −10x0 −1 − λ =−λ[λ2+ (3x20− 5.2x0+ kgs(1 + e−10(x0+0.25)) + 10kgs(x0− 2)e−10(x0+0.25) (1 + e−10(x0+0.25))2 + 1)λ +(3x20+ 4.8x0+ kgs(1 + e−10(x0+0.25)) + 10kgs(x0− 2)e−10(x0+0.25) (1 + e−10(x0+0.25))2 )]
Now, we consider the characteristic equation of the linearized equation that is im-portant to determine the stability and bifurcation points of the corresponding critical manifold. Let Dkgs(x0, λ) = λ 2 + (3x20− 5.2x0+ kgs(1 + e−10(x0+0.25)) + 10kgs(x0− 2)e−10(x0+0.25) (1 + e−10(x0+0.25))2 + 1)λ +(3x20+ 4.8x0 + kgs(1 + e−10(x0+0.25)) + 10kgs(x0− 2)e−10(x0+0.25) (1 + e−10(x0+0.25))2 ), Skgs(x0) = 3x 2 0− 5.2x0+ kgs(1 + e−10(x0+0.25)) + 10kgs(x0− 2)e−10(x0+0.25) (1 + e−10(x0+0.25))2 + 1 and Pkgs(x0) = 3x 2 0+ 4.8x0+ kgs(1 + e−10(x0+0.25)) + 10kgs(x0− 2)e−10(x0+0.25) (1 + e−10(x0+0.25))2
then the sum of the roots of Dkgs(x0, λ) is−Skgs(x0) and the product of the roots of Dkgs(x0, λ) is Pkgs(x0) for any x0 ∈ Nkgs.
From the Figure2-5, we have R+−for x5 < x0 < x6. C++ for x3 < x0 < x4. C−− for x7 < x0 < x3 or for x4 < x0 < x8.
R−−for x0 < x5, for x6 < x0 < x7 or for x8 < x0.
The number of characteristic roots with positive or negative real parts of Dkgs(x0, λ) =
0 with x0 in different regions, where R−−(R++) denotes two negative(positive) real
roots, R+− denotes one positive and one negative real roots, C−−(C++) denotes two
−2 −1 0 1 2 3 4 −10 −8 −6 −4 −2 0 2 4 6 8 10 x 0 −S(x 0 ) and P(x 0 ) x 5= −1.6 x 5 x 1= −0.36347 x 1 x 2= −0.074521 x 2 x 6= 0.026356 x6 x 3= 0.47396 x 3 x 4= 1.2554 x 4 x 7= 0.037626 x7 x 8= 2.9775 x 8 S(x i)=0 for i=1,2,3,4 P(x i)=0 for i=5,6 (S2−P)(x i)=0 for i=7,8 kg s = 0.8 −S(x 0) P(x0) S2−4P
Figure 2: −Skgs(x0) and Pkgs(x0) for kgs = 0.8.
−2 −1 0 1 2 3 4 −10 −8 −6 −4 −2 0 2 4 6 8 10 x 0 −S(x 0 ) and P(x 0 ) x 5= −1.6 x 5 x 1= −0.36529 x 1 x 2= −0.072958 x 2 x 6= 0.026473 x 6 x 3= 0.47836 x 3 x 4= 1.2511 x 4 x 7= 0.037649 x 7 x 8= 2.9766 x 8 S(x i)=0 for i=1,2,3,4 P(x i)=0 for i=5,6 (S2−P)(xi)=0 for i=7,8 kgs = 0.81 −S(x 0) P(x 0) S2−4P
−2 −1 0 1 2 3 4 −10 −8 −6 −4 −2 0 2 4 6 8 10 x 0 −S(x 0 ) and P(x 0 ) x 5= −1.6 x 5 x 1= −0.36882 x 1 x 2= −0.069959 x 2 x 6= 0.026702 x6 x 3= 0.48731 x 3 x 4= 1.2423 x 4 x 7= 0.037694 x 7 x 8= 2.9747 x 8 S(x i)=0 for i=1,2,3,4 P(x i)=0 for i=5,6 (S2−P)(x i)=0 for i=7,8 kg s = 0.83 −S(x 0) P(x0) S2−4P
Figure 4: −Skgs(x0) and Pkgs(x0) for kgs = 0.83.
−2 −1 0 1 2 3 4 −10 −8 −6 −4 −2 0 2 4 6 8 10 x 0 −S(x 0 ) and P(x 0 ) x 5= −1.6 x 5 x 1= −0.37704 x 1 x2= −0.063116 x2 x 6= 0.027245 x6 x 3= 0.51063 x3 x4= 1.2194 x 4 x 7= 0.037803 x 7 x8= 2.9701 x 8 S(x i)=0 for i=1,2,3,4 P(x i)=0 for i=5,6 (S2−P)(x i)=0 for i=7,8 kg s = 0.88 −S(x 0) P(x 0) S2−4P
Let M0(kgs; δ) ={ (x, y, z) ∈ R3 | 2.6x2− x3+ y− z + 4 − kgs(x− 2) 1 1 + e−10(x+0.25) = 0, −y − 5x 2 + 1 = 0, x∈ [−2, 2.5] but x /∈ 6 ∪ i=3 B(xi; δ)}
be a critical manifold which satisfies the Hypotheses1,2and 5of Eq. (16) from Nkgs
for some 0 < δ ≪ 1 then we can apply Geometric Singular Perturbation Theory to Eq. (16).
4.2 Dynamics on the fast subsystem (
17a
)
In this subsection, we consider the dynamics on the fast subsystem (17a) with fixed kgs and treat z as the bifurcation parameter in [−5, 8].
4.2.1 kgs= 0.8
First, we consider kgs= 0.8. When z varies, we discover twelve types of dynamical
behaviors in numerical simulation. The tables (2– 4) are the range of these types. Table 2: The dynamics of the fast subsystem (17a) with various range of z. x4 ≈
1.2554(z≈ −0.1653), x5 ≈ −1.6(z ≈ 2.9520).
z −5 ≤ z < −0.1653 −0.1653 < z < 2.9520 2.9520 < z≤ 3.2 3.3≤ z ≤ 3.4
Types Type I (see Figure (6(a))) Type II (see Figure (6(b))) Type III (see Figure (6(c))) Type IV (see Figure (6(d))
Table 3: The dynamics of the fast subsystem (17a) with various range of z. 4.2 <
zHB1 < 4.3.
z 3.5≤ z ≤ 4.1 zHB1 4.3≤ z ≤ 4.8 4.9≤ z ≤ 5.2
Table 4: The dynamics of the fast subsystem (17a) with various range of z. 5.2 <
zHB2 < 5.3, x3 ≈ 0.47396(z ≈ 5.5744) and x6 ≈ 0.026356(z ≈ 6.4836).
z zHB2 5.3≤ z < 5.5744 5.744 < z < 6.4836 6.4836 < z≤ 8
Types Type IX (see Figure (7(c)) Type X (see Figure (7(d)) Type XI (see Figure (7(e)) Type XII (see Figure (7(f))
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(a) Type I; z = −0.5 has only one stable fixed
point −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(b) Type II; z = 2 has only one unstable fixed point and only one stable periodic orbit
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(c) Type III; z = 3.2 has three fixed point and only one stable periodic orbit
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (d) Type IV; z = 3.4. −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (e) Type V; z = 3.8 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(f) Type VI; zHB1 is a parameter for Homoclinic
Bifurcation, where 4.2 < zHB1< 4.3.
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(a) Type VII; z = 4.6.
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(b) Type VIII; z = 5 has two periodic orbits.
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(c) Type IX; zHB2 is a parameter for Homoclinic
Bifurcation, where 5.2 < zHB2< 5.3. −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (d) Type X; z = 5.4. −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 −35 −30 −25 −20 −15 −10 −5 0 5 x y
(e) Type XI; z = 6.
−3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 −35 −30 −25 −20 −15 −10 −5 0 5 x y (f) Type XII; z = 7.
4.2.2 kgs= 0.81
Second, we consider kgs = 0.81. When z varies, we discover thirteen types of
dynamical behaviors in numerical simulation. The tables (5 – 8) are the range of these types.
Table 5: The dynamics of the fast subsystem (17a) with various range of z. x4 ≈
1.2511(z≈ −0.1081), x5 ≈ −1.6(z ≈ 2.9520).
z −5 ≤ z < −0.1081 −0.1081 < z < 2.9520 2.9520 < z≤ 3.31 3.32≤ z ≤ 3.46
Types Type I (see Figure (8(a))) Type II (see Figure (8(b))) Type III (see Figure (8(c))) Type IV (see Figure (8(d))
Table 6: The dynamics of the fast subsystem (17a) with various range of z. 4.2 <
zHB1 < 4.3.
z 3.47≤ z ≤ 3.56 3.6≤ z ≤ 4.2 zHB1 4.3≤ z ≤ 4.8
Types Type V (see Figure (8(e)) Type VI (see Figure (8(f)) Type VII (see Figure (9(a)) Type VIII (see Figure (9(b))
Table 7: The dynamics of the fast subsystem (17a) with various range of z. x3 ≈
0.47836(z≈ 5.5730), x6 ≈ 0.026473(z ≈ 6.5021) and 5.2 < zHB2 < 5.3.
z 4.9≤ z ≤ 5.2 zHB2 5.3≤ z < 5.5730 5.5730 < z < 6.5021
Types Type IX (see Figure (9(c)) Type X (see Figure (9(d)) Type XI (see Figure (9(e)) Type XII (see Figure (9(f))
Table 8: The dynamics of the fast subsystem (17a) with various range of z.
z 6.5021 < z≤ 8
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(a) Type I; z = −0.5 has only one stable fixed
point −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(b) Type II; z = 2 has only one unstable fixed point and only one stable periodic orbit
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(c) Type III; z = 3.2 has three fixed point and only one stable periodic orbit
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (d) Type IV;z = 3.4 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (e) Type V; z = 3.5 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (f) Type VI; z = 3.8
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(a) Type VII; zHB1 is a parameter for Homoclinic
Bifurcation, where 4.2 < zHB1< 4.3 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (b) Type VIII; z = 4.6 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (c) Type IX; z = 5 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(d) Type X; zHB2 is a parameter for Homoclinic
Bifurcation, where 5.2 < zHB2< 5.3 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(e) Type XI; z = 5.4
−3 −2 −1 0 1 2 3 −35 −30 −25 −20 −15 −10 −5 0 5 10 x y (f) Type XII; z = 6 −3 −2 −1 0 1 2 3 −35 −30 −25 −20 −15 −10 −5 0 5 10 x y (g) Type XIII; z = 7
4.2.3 kgs= 0.83
Third, we consider kgs = 0.83. When z varies, we discover twelve types of
dy-namical behaviors in numerical simulation. The tables (9–11) are the range of these types.
Table 9: The dynamics of the fast subsystem (17a) with various range of z. x4 ≈
1.2423(z≈ 0.0075), x5 ≈ −1.6(z ≈ 2.9520).
z −5 ≤ z < 0.0075 0.0075 < z < 2.9520 2.9520 < z≤ 3.4 3.5≤ z ≤ 3.8
Types Type I (see Figure (10(a))) Type II (see Figure (10(b))) Type III (see Figure (10(c))) Type IV (see Figure (10(d))
Table 10: The dynamics of the fast subsystem (17a) with various range of z. 4.3 <
zHB1 < 4.4.
z 3.9≤ z ≤ 4.3 zHB1 4.4≤ z ≤ 4.6 4.7≤ z ≤ 5.2
Types Type V (see Figure (10(e)) Type VI (see Figure (10(f)) Type VII (see Figure (11(a)) Type VIII (see Figure (11(b))
Table 11: The dynamics of the fast subsystem (17a) with various range of z. 5.2 <
zHB2 < 5.3, x3 ≈ 0.48731(z ≈ 5.5691) and x6 ≈ 0.026702(z ≈ 6.5393).
z zHB2 5.3≤ z < 5.5691 5.5691 < z < 6.5393 6.5393 < z≤ 8
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(a) Type I; z = −0.5 has only one stable fixed
point −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(b) Type II; z = 2 has only one unstable fixed point and only one stable periodic orbit
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(c) Type III; z = 3.2 has three fixed point and only one stable periodic orbit
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (d) Type IV; z = 3.5 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (e) Type V; z = 4 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x
(f) Type VI; zHB1 is a parameter for Homoclinic
Bifurcation, where 4.3 < zHB1< 4.4
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(a) Type VII; z = 4.6
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (b) Type VIII; z = 5 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(c) Type IX; zHB2 is a parameter for Homoclinic
Bifurcation, where 5.2 < zHB2< 5.3 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (d) Type X; z = 5.4 −3 −2 −1 0 1 2 3 −35 −30 −25 −20 −15 −10 −5 0 5 10 x y
(e) Type XI; z = 6
−3 −2 −1 0 1 2 3 −35 −30 −25 −20 −15 −10 −5 0 5 10 x y (f) Type XII; z = 7
4.2.4 kgs= 0.88
Finally, we consider kgs= 0.8. When z varies, we discover eleven types of
dynam-ical behaviors in numerdynam-ical simulation. The tables (12 – 14) are the range of these types.
Table 12: The dynamics of the fast subsystem (17a) with various range of z. x4 ≈
1.2194(z≈ 0.3047), x5 ≈ −1.6(z ≈ 2.9520).
z −5 ≤ z < 0.3047 0.3047 < z < 2.9520 2.9520 < z≤ 3.5 3.6≤ z ≤ 3.7
Types Type I (see Figure (12(a))) Type II (see Figure (12(b))) Type III (see Figure (12(c))) Type IV (see Figure (12(d))
Table 13: The dynamics of the fast subsystem (17a) with various range of z. 4.3 <
zHB1 < 4.4, 5.1 < zHB2< 5.2.
z 3.8≤ z ≤ 4.3 zHB1 4.4≤ z ≤ 5.1 zHB2
Types Type V (see Figure (12(e)) Type VI (see Figure (12(f)) Type VII (see Figure (13(a)) Type VIII (see Figure (13(b))
Table 14: The dynamics of the fast subsystem (17a) with various range of z. x3 ≈
0.51063(z≈ 5.5511) and x6 ≈ 0.027245(z ≈ 6.6321).
z 5.2≤ z < 5.5511 5.5511 < z < 6.6321 6.6321 < z <≤ 8
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(a) Type I; z = −0.5 has only one stable fixed
point −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(b) Type II; z = 2 has only one unstable fixed point and only one stable periodic orbit
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(c) Type III; z = 3.2 has three fixed point and only one stable periodic orbit
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (d) Type IV; z = 3.6 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (e) Type V; z = 4 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(f) Type VI; zHB1 is a parameter for Homoclinic
Bifurcation, where 4.3 < zHB1< 4.4
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(a) Type VII; z = 4.6
−3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y
(b) Type VIII; zHB2is a parameter for Homoclinic
Bifurcation, where 5.1 < zHB2< 5.2 −3 −2 −1 0 1 2 3 −30 −25 −20 −15 −10 −5 0 5 10 x y (c) Type IX; z = 5.4 −3 −2 −1 0 1 2 3 −35 −30 −25 −20 −15 −10 −5 0 5 10 x y (d) Type X; z = 6 −3 −2 −1 0 1 2 3 −35 −30 −25 −20 −15 −10 −5 0 5 10 x y
(e) Type XI; z = 7
4.3 Numerical results for various of kg
sIn this subsection, we show some pictures which associates with a solution orbit of Eq. (16).
(I) First we consider the case that kgs< 0.808. From table (1), the attractor is unique
regular bursting. In numerical simulation, we choose kgs= 0.8 with two distinct
ini-tial values x0(i) = (x(i)0 , y0(i), z0(i)), i = 1, 2 (See Figures 14 and 15), where z0(1) in Type I (see Figure (6(a)) and z0(2) in Type VII (see Figure (7(a)).
(II) Second we consider the case that 0.809 ≤ kgs ≤ 0.813. From table (1), the
attractors have only two types regular bursting and periodic solution. In numerical simulation, we choose kgs = 0.81 with two distinct initial values x
(i) 0 = (x (i) 0 , y (i) 0 , z (i) 0 ),
i = 1, 2 (See Figures 16and 17), where z(1)0 in Type XI (see Figure (9(e)) and z0(2) in Type I (see Figure (8(a)).
(III) Third we consider the case that 0.814≤ kgs ≤ 0.85. From table (1), the
attrac-tors have only two types regular bursting and steady state. In numerical simulation, we choose kgs= 0.83 with three distinct initial values x
(i) 0 = (x (i) 0 , y (i) 0 , z (i) 0 ), i = 1, 2, 3
(See Figures18, 19and 20), where z0(1) in Type V (see Figure (10(e)), z0(2) in Type I (see Figure (10(a)) and z0(3) in Type XII (see Figure (11(f)).
(IV) Finally we consider the case that kgs ≥ 0.87. From table (1), the attractor is
unique steady state. In numerical simulation, we choose kgs = 0.88 with two distinct
initial values x0(i) = (x(i)0 , y(i)0 , z0(i)), i = 1, 2 (See Figures 21 and 22), where z0(1) in Type VII (see Figure (13(a)) and z0(2) in Type I (see Figure (12(a)).
Note:
1. for kgs = 0.8, we have x3 ≈ 0.47396, x4 ≈ 1.2554, x5 ≈ −1.6, x6 ≈ 0.026356,
x7 ≈ 0.037626 and the unstable fixed point is roughly (0.020838, 0.99783, 6.4834).
2. for kgs = 0.81, we have x3 ≈ 0.47836, x4 ≈ 1.2511, x5 ≈ −1.6, x6 ≈ 0.026473,
x7 ≈ 0.037649 and the unstable fixed point is roughly (0.02553, 0.99674, 6.5021).
3. for kgs = 0.83, we have x3 ≈ 0.48731, x4 ≈ 1.2423, x5 ≈ −1.6, x6 ≈ 0.026702,
x7 ≈ 0.037694 and the unstable fixed point is roughly (0.034704, 0.99398, 6.5388).
4. for kgs = 0.88, we have x3 ≈ 0.51063, x4 ≈ 1.2194, x5 ≈ −1.6, x6 ≈ 0.027245,
0 50 100 150 200 250 300 350 400 −1 0 1 t x 0 50 100 150 200 250 300 350 400 −15 −10 −5 0 t y 0 50 100 150 200 250 300 350 400 −2 0 2 4 t z
(a) Time Series
−5 0 5 10 −4 −2 0 2 −20 −15 −10 −5 0 5 10 x z y (b) Phase portrait −15 −10 −5 0 5 10 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 z x 0.47396 1.2554 −1.6 0.026356 0.037626 (c) Projection on z− x plane −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −30 −25 −20 −15 −10 −5 0 5 x y 0.47396 1.2554 −1.6 0.026356 0.037626 (d) Projection on x− y plane
Figure 14: For kgs = 0.8, a solution orbit of Eq. (16) with initial values x (1) 0 = 1.48, y0(1) =−9.952, z(1)0 =−3.0828. 0 100 200 300 400 500 600 −2 −1 0 1 t x 0 100 200 300 400 500 600 −20 −10 0 t y 0 100 200 300 400 500 600 4 6 t z
(a) Time Series
−1 0 1 2 3 4 5 6 7 8 9 −4 −2 0 2 −40 −30 −20 −10 0 10 z x y (b) Phase portrait 1 2 3 4 5 6 7 8 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 z x 0.47396 1.2554 −1.6 0.026356 0.037626 (c) Projection on z− x plane −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −40 −35 −30 −25 −20 −15 −10 −5 0 5 x y 0.47396 1.2554 −1.6 0.026356 0.037626 (d) Projection on x− y plane
Figure 15: For kgs = 0.8, a solution orbit of Eq. (16) with initial values x (2)
0 = 0.7,
0 100 200 300 400 500 600 0.01 0.02 0.03 0.04 0.05 t x 0 100 200 300 400 500 600 0.9880.99 0.992 0.994 0.996 0.998 t y 0 100 200 300 400 500 600 6.495 6.5 t z
(a) Time Series
6.494 6.496 6.498 6.5 6.502 6.504 6.506 0 0.01 0.02 0.03 0.04 0.05 0.06 0.985 0.99 0.995 1 1.005 z x y (b) Phase portrait 6.494 6.496 6.498 6.5 6.502 6.504 6.506 6.508 −0.01 0 0.01 0.02 0.03 0.04 0.05 z x 0.026473 0.037649 (c) Projection on z− x plane 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.985 0.99 0.995 1 x y 0.026473 0.037649 (d) Projection on x− y plane
Figure 16: For kgs = 0.81, a solution orbit of Eq. (16) with initial values x (1) 0 = 0.04, y0(1) = 0.992, z(1)0 = 6.5009. 0 100 200 300 400 500 600 −1 0 1 t x 0 100 200 300 400 500 600 −15 −10 −5 0 t y 0 100 200 300 400 500 600 −2 0 2 4 t z
(a) Time Series
−15 −10 −5 0 5 10 −4 −2 0 2 −40 −30 −20 −10 0 10 z x y (b) Phase portrait −4 −2 0 2 4 6 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 z x 0.47836 1.2511 −1.6 0.026473 0.037649 (c) Projection on z− x plane −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −40 −35 −30 −25 −20 −15 −10 −5 0 5 x y 0.47836 1.2511 −1.6 0.026473 0.037649 (d) Projection on x− y plane
Figure 17: For kgs = 0.81, a solution orbit of Eq. (16) with initial values x(2)0 = 1.48,
0 100 200 300 400 500 600 0.2 0.4 0.6 0.8 t x 0 100 200 300 400 500 600 −2 −1 0 t y 0 100 200 300 400 500 600 4 5 6 t z
(a) Time Series
0 1 2 3 4 5 6 7 −2 −1 0 1 2 −15 −10 −5 0 5 z x y (b) Phase portrait 3.5 4 4.5 5 5.5 6 6.5 7 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 z x 0.48731 0.026702 0.037694 (c) Projection on z− x plane −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 x y 0.48731 0.026702 0.037694 (d) Projection on x− y plane
Figure 18: For kgs = 0.83, a solution orbit of Eq. (16) with initial values x (1) 0 = 0.8, y0(1) =−2.2, z0(1) = 3.948. 0 100 200 300 400 500 600 −1 0 1 t x 0 100 200 300 400 500 600 −15 −10 −5 0 t y 0 100 200 300 400 500 600 −2 0 2 4 t z
(a) Time Series
−15 −10 −5 0 5 10 −4 −2 0 2 −40 −30 −20 −10 0 10 z x y (b) Phase portrait −4 −2 0 2 4 6 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 z x 0.48731 1.2423 −1.6 0.026702 0.037694 (c) Projection on z− x plane −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −40 −35 −30 −25 −20 −15 −10 −5 0 5 x y 0.48731 1.2423 −1.6 0.026702 0.037694 (d) Projection on x− y plane
Figure 19: For kgs = 0.83, a solution orbit of Eq. (16) with initial values x(2)0 = 1.48,
0 100 200 300 400 500 600 700 800 0.02 0.04 t x 0 100 200 300 400 500 600 700 800 0.9880.99 0.992 0.994 0.996 0.998 t y 0 100 200 300 400 500 600 700 800 6.535 6.54 t z
(a) Time Series
6.53 6.535 6.54 6.545 0 0.02 0.04 0.06 0.08 0.1 0.985 0.99 0.995 1 z x y (b) Phase portrait 6.530 6.532 6.534 6.536 6.538 6.54 6.542 6.544 0.01 0.02 0.03 0.04 0.05 0.06 z x 0.026702 0.037694 (c) Projection on z− x plane 0 0.01 0.02 0.03 0.04 0.05 0.06 0.985 0.99 0.995 1 1.005 x y 0.026702 0.037694 (d) Projection on x− y plane
Figure 20: For kgs = 0.83, a solution orbit of Eq. (16) with initial values x (3) 0 = 0.039603, y(3)0 = 0.99412, z0(3) = 6.5427. 0 100 200 300 400 500 600 −2 −1 0 1 t x 0 100 200 300 400 500 600 −20 −10 0 t y 0 100 200 300 400 500 600 4 6 t z
(a) Time Series
2 3 4 5 6 7 8 −5 0 5 −30 −25 −20 −15 −10 −5 0 5 10 z x y (b) Phase portrait 0 1 2 3 4 5 6 7 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 z x 0.51063 1.2194 −1.6 0.027245 0.037803 (c) Projection on z− x plane −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −40 −35 −30 −25 −20 −15 −10 −5 0 5 x y 0.51063 1.2194 −1.6 0.027245 0.037803 (d) Projection on x− y plane
Figure 21: For kgs = 0.88, a solution orbit of Eq. (16) with initial values x(1)0 =−0.38,
0 100 200 300 400 500 600 0 0.5 1 1.5 t x 0 100 200 300 400 500 600 −8 −6 −4 −20 t y 0 100 200 300 400 500 600 −2 0 2 4 6 t z
(a) Time Series
−15 −10 −5 0 5 10 −4 −2 0 2 −40 −30 −20 −10 0 10 z x y (b) Phase portrait −4 −2 0 2 4 6 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 z x 0.51063 1.2194 −1.6 0.027245 0.037803 (c) Projection on z− x plane −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 −15 −10 −5 0 5 x y 0.51063 1.2194 0.027245 0.037803 (d) Projection on x− y plane
Figure 22: For kgs = 0.88, a solution orbit of Eq. (16) with initial values x (2)
0 = 1.48,
y0(2) =−9.952, z(2)0 =−3.0412.
4.4 Analysis and conclusions
First, by The Poincaré–Andronov–Hopf Bifurcation Theorem, we have a Hopf bifurcation at kgs ≈ 0.8131 and there exists stable periodic solution to kgs < 0.8131
(see Figure (23)).
Second, from subsection 4.1, we have
M0(kgs; δ) ={ (x, y, z) ∈ R3 | 2.6x2− x3+ y− z + 4 − kgs(x− 2) 1 1 + e−10(x+0.25) = 0, −y − 5x 2 + 1 = 0, x∈ [−2, 2.5] but x /∈ 6 ∪ i=3 B(xi; δ)}
is a critical manifold which satisfies the Hypotheses1, 2and 5of Eq. (16) from Nkgs
for some chosen 0 < δ ≪ 1 (e.g. δ = 0.0001) then we can apply Fenichel’s Second Theorem (theorem (3)) to Eq. (16). From subsection 4.2, we have the stable and unstable manifold of the equilibriums of the fast subsystem Eq. (36). So we know what are Ws(M
0) and Wu(M0). Thus by Fenichel’s Second Theorem (Theorem (3)),
there exists manifolds Ws(M
ε) and Wu(Mε) that are O(ε) close and are diffeomorphic
to Ws(M
0) and Wu(M0) respectively, and they are each locally invariant under the
that usually involve bursting can be explained (Figure14,15,17,18, 19, 21and 22). This illustrates the power of using Singular perturbation theory to understand the global dynamical properties of realistic biological systems.
0.8 0.82 0.84 0.86 0.88 −1 −0.5 0 kg s X1 0.8 0.82 0.84 0.86 0.88 −1 0 1 kg s Y1 0.8 0.82 0.84 0.86 0.88 −0.5 0 0.5 kg s X2 0.8 0.82 0.84 0.86 0.88 0 0.5 1 kg s Y2 0.8 0.82 0.84 0.86 0.88 −0.5 0 0.5 kgs X3 0.8 0.82 0.84 0.86 0.88 −0.6 −0.4 −0.2 kgs Y3
Figure 23: The eigenvalues (λ1 = x1 + iy1, λ2 = x2 + iy2 and λ3 = x3 + iy3) of
equilibrium of (16) to kgs. Reλ2 = 0 and Reλ3 = 0 at kgs ≈ 0.8131.
References
[1] J.L. Hindmarsh and R.M. Rose, “A model of neuronal bursting using three coupled first order differential equations,” Proc. Roy Soc. London Ser. B 221 (1984), 87-102.
[2] D. Terman, “Chaotic spikes arising from a model for bursting in excitable membranes,” SIAM J. Appl. Math. 51 (1991), 1418-1450.
[3] D. Terman, “The transition from bursting to continuous spiking in an excitable membrane model,” J. Nonlinear Sci. 2 (1992), 133-182.
[4] J. Rinzel and G.B. Ermentrout, “Analysis of neural excitability and oscillations Methods in Neuronal Modeling: From Ions to Networks,” 2nd edn., C. Koch and I. Segev, eds, MIT Press, Cambridge, MA (1998), 251-291.
[5] G. Innocenti, A. Morelli, R. Genesio, and A. Torcini, “Dynamical phases of the Hindmarsh-Rose neuronal model: Studies of the transition from bursting to spiking chaos,” Chaos, vol. 17, no. 4, p. 043128, 2007.
[6] E. M. Izhikevich, “Which model to use for cortical spiking neurons?,” IEEE Trans. Neural Netw., vol. 15, no. 5, pp. 1063–1070, Sep. 2004.
[7] M. Storace, D. Linaro, and E. de Lange, “The Hindmarsh-Rose neuron model: Bifurcation analysis and piecewise-linear approximations,” Chaos, vol. 18, no. 3, p. 033128, 2008.
[8] E. de Lange and M. Hasler, “Predicting single spikes and spike patterns with the Hindmarsh-Rose model,” Biol. Cybern., vol. 99, no. 4, pp. 349–360, 2008.
[9] G. Innocenti and R. Genesio, “On the dynamics of chaotic spiking-bursting transition in the Hindmarsh-Rose neuron,” Chaos, vol. 19, no. 2, p. 023124, 2009.
[10] D. Golomb and J. Rinzel, “Clustering in globally coupled inhibitory neurons,” Phys. D 72 (1994), 259-282.
[11] D. Golomb, X.-J. Wang and J. Rinzel, “Synchronization properties of spindle oscillations in a thalamic reticular nucleus model,” J. Neurophysiol. 72 (1994), 1109-1126.
[12] N. Kopell and G. LeMasson, “Rhythmogenesis, amplitude modulation and multiplexing in a cortical architecture,” Proc. Nat. Acad. Sci. USA 91 (1994), 10586-10590.
[13] I. Belykh and M. Hasler, “Nonlinear dynamics: New directions,” (Springer, 2013 (in press)) Chap. Patterns of synchrony in neuronal networks: the role of synaptic inputs.
[14] Jonq Juang and Yu-Hao Liang, “Cluster Synchronization in Networks of Neurons with Chemical Synapses,” Chaos, vol. 24, no. 1, p. 013110, 2014.
[15] W. Singer, “Neuronal synchrony: A versatile code for the definition of relations?,” Neuron , vol. 24, no. 1, pp. 49–65, Sep. 1999.
[16] P. J. Uhlhaas and W. Singer, “Neural synchrony in brain disorders: Relevance for cognitive dysfunctions and pathophysiology,” Neuron, vol. 52, no. 1, pp. 155–168, Oct. 2006.
[17] S. Neuenschwander and W. Singer, “Long-range synchronization of oscillatory light responses in the cat retina and lateral geniculate nucleus,” Nature, vol. 379, pp. 728–732, 1996.
[18] L. Glass, “Synchronization and rhythmic processes in physiology,” Nature, vol. 410, no. 6825, pp. 277–284, 2001.
[19] C. M. Gray, P. Konig, A. K. Engel, and W. Singer, “Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization which reflects global stimulus properties,” Nature, vol. 338, pp. 334–337, 1989.
[20] R. Dzakpasu and M. Zochowski, “Discriminating differing types of synchrony in neural systems,” Phys. D, vol. 208, no. 1–2, pp. 115–122,2005.
[21] G. Buzsaki, Rhythms of the Brain. New York: Oxford Univ. Press,2006.
[22] C. J. Stam, B. F. Jones, G. Nolte, M. Breakspear, and S. Ph, “Small-world networks and functional connectivity in alzheimer’s disease,” Cereb. Cortex, vol. 17, no. 1, pp. 92–99, 2006. [23] Jonq Juang, Yu-Hao Liang and Fang-Jhu Jhou, “Multistate and Multistage Synchronization of Hindmarsh-Rose Neurons With Excitatory Chemical and Electrical Synapses,” Circuits and Systems I: Regular Papers, IEEE Transactions on (Volume:59 , Issue: 6 ), June 2012, 1335 -1347
[24] Christopher K. R. T. Jones, Geometric singular perturbation theory. In: Johnson R (ed) Dy-namical systems, Montecatibi Terme, Lecture Notes in Mathematics, vol 1609. Springer, Berlin, pp 44–118 (1995)
[25] Geertje Hek, “Geometric singular perturbation theory in biological practice,” J. Math. Biol. (2010) 60:347–386, DOI 10.1007/s00285-009-0266-7
[26] Stephen Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Texts in applied mathematics, vol. 2, Springer, Second Edition, 2003
[27] Yuri A. Kuznetsov, Elements of Applied Bifurcation Theory, Applied Mathematical Sciences, vol. 112, Springer, Third Edition, 2004.
[28] Y. G. Zheng and Z. H. Wang “Time-delay effect on the bursting of the synchronized state of coupled Hindmarsh-Rose neurons,” Chaos, vol. 22, no. 4, p. 043127, 2012.
