化學與電子耦合的 Hindmarsh-Rose 神經細胞其產生的多狀態、層級同步現象

國 立 交 通 大 學

應用數學系

數學建模與科學計算碩士班

碩

士

論

文

化學與電子耦合的 Hindmarsh-Rose

神經細胞其產生的多狀態、層級同步現象

Multi-state and multi-stage synchronization of

Hindmarsh-Rose neurons with chemical

and electrical synapses

研 究 生：周芳竹

指導老師：莊 重 教授

化學與電子耦合的 Hindmarsh-Rose

神經細胞其產生的多狀態、層級同步現象

Multi-state and multi-stage synchronization of

Hindmarsh-Rose neurons with chemical

and electrical synapses

研 究 生：周芳竹 Student：Fang-Jhu Jhou

指導教授：莊 重 Advisor：Jonq Juang

國 立 交 通 大 學

應用數學系數學建模與科學計算碩士班

碩 士 論 文

A Thesis

Submitted to Department of Applied Mathematics College of Science,

Institute of Mathematical Modeling and Scientific Computing

National Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of

Master

In

Applied Mathematics

June 2010

Hsinchu, Taiwan, Republic of China

i

誌 謝

首先，感謝我的指導教授莊重老師，在研究所的期間所給予在課

業方面的教導。同時引領我進入做研究的世界，讓我對數學研究

工作有了更進一步的理解，也對我無論在數學以及其他的層面上

都有著深刻的啟發。另外也要感謝成長的路上所有曾經教導過我

的老師們，激發我對數學的好奇心，讓我能夠有機會踏入數學領

域。感謝所有學長姐以及同學朋友們，謝謝各位所給予的支持與

鼓勵。然而，在論文研究過程和準備口試的階段特別感謝梁育豪

學長的指導，讓我的論文得以更順利的完成。最後，我要感謝我

的家人，因為有你們在成長的路上的支持與陪伴，讓我能夠順利

的完成這篇論文。此份榮耀我要與你們分享，謹以此文表達我最

誠摯的謝意。

ii

化學與電子耦合的 Hindmarsh-Rose 神經細胞其產生

的多狀態、層級同步現象

學生：周芳竹 指導老師：莊 重 教授

國立交通大學應用數學系數學建模與科學計算碩士班

摘 要

本文研究化學與電子耦合的 Hindmarsh-Rose (HR) 神經細胞產

生的多狀態同步(混沌共存的穩定和周期/穩態同步)和多層級同步複

雜的神經網絡。我們明確的得到 HR 神經元網路的完全同步狀態的全

局穩定和局部穩定的同步區間。同時有化學和電子耦合才會有共存的

穩定多狀態同步，包括混沌同步、週期/穩態同步。相比之下耦合振

盪器系統或耦合映象格子模型只有單一狀態的同步現像。對於局部同

步，我們得到即使沒有電子耦合，不管化學耦合的連結網络是如何稀

疏，耦合神經元可以達到穩定穩態，我們預測最小的化學耦合強度和

每個神經元接收的信號成正比。此外，我們提供化學耦合的連結網络

要如何的濃密，化學突觸才會扮演增強同步化的角色來達到系統的同

步。另外我們證明 HR 耦合神經元是有界耗散，來得到全局同步的區

域，這樣的結果可以應用於一般的單神經元模型和複雜的網絡。我們

的方法很一般，例如，它可以立即應用到其他單神經元模型，如

FitzHugh-Nagumo 模型。所需的工具和概念，包括坐標變換，矩陣

測量，單調動力學和時間平均估計。

Multi-state and Multi-stage Synchronization of

Hindmarsh-Rose Neurons

with Chemical and Electrical Synapses

Student: Fang-Jhu Jhou

Advisors: Jonq Juang

Institute of Mathematical and Scientific Computing National Chiao Tung University

Hsinchu, Taiwan, R.O.C.

Abstract

The multi-state synchronization, the coexistence of stably chaotic and periodic/steady-state synchronization, and multi-stage synchronization of Hindmarsh-Rose(HR) neurons with both chemical and electrical synapses over the complex network are analytically studied. The synchronization regions for both global and local stability of the com-pletely synchronous state in such networks of HR neurons are explicitly obtained. The coexistence of the stably multi-state synchronization, including chaotic synchronization and periodic/steady-state synchronization, is provided with the presence of both chem-ical and electrchem-ical synapses in the network. This is in contrast with Coupled Oscillator Systems or Coupled Map Lattices where only single-state synchronization is found. For local synchronization, we are able to show that even without electrical coupling, the coupled neurons may reach stably steady-state synchrony regardless of how sparsely the chemically coupling networks is coupled and that the minimum chemically coupling bound as predicted is inversely proportional to the number of signals each neuron re-ceives. Moreover, we provide a measurement of how densely coupled the system should be so as to have the chemical synapses to play an enforcing role in achieving the syn-chrony of the system. We establish the bounded dissipation of the coupled HR neurons, followed by the attainment of its global synchronization region. Such a result can be applied to a general class of single neuron models and the complex networks. Our method employed here is quite general. For instance, it can be immediately applied to other single neuron models such as the FitzHugh-Nagumo model. The analytical tools and concepts needed include coordinate transformations, matrix measures, monotone dynamics and time averaging estimates.

Contents

1

Introduction

The fundamental building block of every nervous system is the neuron. There is an increasing trend towards studying the dynamical behavior of relatively large networks of neurons, and model-ing/emulating such networks is also on the rise. Neural synchronization has been suggested as par-ticularly relevant for neuronal signal transmission and coding in the brain. Brain [1, 2, 3, 4, 5, 6, 7] oscillations 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 this paper, the phenomena of multi-state and multi-stage synchronization of Hindmarsh-Rose neurons with chemical and electrical synapses over complex network are analytical studied. These are in contract with coupled oscillator systems or coupled map lattices where only single-state synchro-nization is found. It should be noted that chemical synapses or electrical synapses alone can not produce such new phenomena. Our melthod employed here is quite general. For instance, it can be applied to other single neuron models such as the FitzHugh-Nagumo model.

Many biological neuron models have been developed in the last decades for an accurate descrip-tion and predicdescrip-tion of biological phenomena. The early works of Hodgkin and Huxley surely represent a milestone. Since such a model turns out to be quite complex, simpler approximations, namely, second order systems such as the FitzHugh-Nagumo(FN), Morris-Lecar and Hindmarsh-Rose(HR) neuronal models have been proposed. However, the second order model is not able to reproduce some interesting phenomena such as terminating itself by triggering a set of stable firings. Hence, the 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. Such model of HR has turned out to be accurate in capturing both qualitative and quantitative aspects of experimental data [9, 10, 11, 12]. Furthermore, it has been shown that the HR neuron model is capable of producing major neuronal behaviors such as spiking, bursting, and chaotic regime [12, 13, 14].

There are about 1010 neurons with an approximate 1014 links between them, all packed in a

human brain. Although neurons are sparsely connected, they are within only a few synaptic steps from all other neurons and their underlying network has small-world property [8]. Neurons in a population synchronize their activity using electrical (via gap junctions) and chemical synapses with other neurons in the same population as well as with neurons from other populations. In the first case, the coupling through gap junctions is linear and directly dependent on the difference of the membrane potentials. In the second case, the coupling is pulsatile and often modeled as a static sigmoidal nonlinear input-output function with a threshold and saturation.

In this work, we study the multi-state synchronization and multi-stage synchronization in en-sembles of electrically and chemically coupled Hindmarsh-Rose neurons whose connection topology with respect to the electrical coupling is allowed to be complex including, e.g., Newman-Wattts networks, and whose couplings through chemical synapses are unidirectional from presynaptic cell to the postsynaptic cell. By multi-state synchronization, we mean the coexistence of stably chaotic and periodic/steady-state synchronization, depending on the choice of initial conditions. By varying a certain parameters, if the corresponding coupled neurons yield different types of synchronization such as chaotic, periodic or steady-state synchronization, then it is said that the system of the cou-pled neurons exhibits the multi-stage synchronization. Our main results contain the following. The regions in terms of chemically and electrically coupling strengths for both local and global stability of the completely synchronous state in such networks of HR neurons are explicitly obtained. The regions depend on the details of the topology of the electrically connected network, the second largest eigenvalue of its associated connection matrix. However, they only depend on the number of sig-nals each neuron received, independent of all other details of chemically coupling network topology. Moreover, with the presence of both chemical and electrical synapses in the network, the coexistence of the stable multi-state synchronization including chaotic synchronization and periodic/steady-state synchronization, is provided. This is in contrast with Coupled Oscillator Systems or Coupled Map Lattices where only single-state synchronization is observed. It should be noted that without the chemical synapses between neurons, the multi-state synchronization would be impossible. It is also

shown that for any kgs> 0, where k is the number of signals each neuron received and gsthe strength

of chemical synapses, there exists a minimum strength of electrical synapses so that the coupled HR neurons synchronize. Since on the synchronous manifold, the dynamics of the synchronous equation goes from chaotic (bursting) to periodic (spiking) to steady-state. Consequently, the phenomena of the multi-stage synchronization is also observed. The above described scenarios can not exist without the presence of electrical synapses between neurons. For local synchronization, we are able to show that even without the electrical coupling, the coupled neurons may reach stably steady-state synchrony regardless of how sparsely the chemically coupling network is coupled and that the minimum chemically coupling bound as predicted is inversely proportional to the number of signals each neuron received. If, in addition, the network is also electrically coupled, then the minimum electrically coupling bound to reach even stably bursting (chaotic) and multi-state synchronization can also be explicitly computed. Moreover, we provide a measurement of how densely coupled the system should be so as to have the chemical synapses to play an enforcing role in achieving the synchrony of the system. For global synchronization, we first establish the bounded dissipation of the coupled HR neurons. The synchronization region is also obtained. Such a result can be ap-plied to a general class of single neuron models and the complex networks. Our method employed here is quite general. For instance, it can be immediately applied to other single neuron models

such as the FitzHugh-Nagumo model. The analytical tools and concepts needed include coordinate transformations, matrix measures, monotone dynamics and time averaging estimates.

The most closely related works to ours are those done by Jalili [15], Kopell and Ermentrout [16] and Belykh, Lange and Hasler [17]. In [16], the single neuron model is a quadratic integrate and fire. They obtained that the chemical and electrical coupling play complementary roles in the coherence of rhythms in inhibitory networks. In [17], densely coupled bursting HR neurons with only chemical coupling was studied. They demonstrated the bound of the minimum chemical strength for obtaining the steady-state synchronization only depends on the numbers of signals each neuron receives, independent of all other details of the network topology. These two works used both numerical and analytical techniques to address local synchronization. Whereas the results in [15], though the same model as ours was studied, was numerical. The surprising phenomena of multi-state synchronization was not mentioned there.

The paper is organized as follows. Section 2 is to lay down the foundation of our paper, which includes the formulation and needed preliminaries, including coordinate transformations, matrix measures, monotone dynamics as well as time averaging estimates. The main results are contained in Section 3. Their detailed applications concerning coupled HR neurons and some concluding remarks are recorded in Section 4.

2

Formulation

The HR model may be seen either as a generalization of the FitzHugh-Nagumo equations or as a simplification of physiologically realistic model proposed by Hodgkin and Huxley. The motion of the model reads as follows:

˙x = f (x) + y − z + q,

˙y = −y − 5x2+ 1, (1)

˙z = µ(b(x − x0) − z).

Here f (x) = ax2 _{− x}3_{, x is the membrane potential, y and z are the recovery(fast) and the}

adaptation(slow) current, respectively. The roles played by the system parameters are roughly the following. q mimics the membrane input current for biological neurons; a allows one to switch be-tween bursting and spiking behaviors and to control the spiking frequency; µ controls the speed of variation of the slow variable z, and in the presence of spiking behaviors, it governs the spiking frequency, whereas in the case of bursting, it affects the number of spikes per burst; b governs adap-tation; a unitary value of b determines spiking behavior without accommodation and subthreshold adaptation, whereas around b = 4 give strong accommodation and subthreshold overshoot, or even

oscillations; x0 sets the resting potential of the system. Hereafter, the parameters are chosen and

fixed as follows: x0 = −1.6, µ = 0.01, b = 4, q = 4, and a = 2.6. The dynamics of the neuron

with such set of parameters is bursting and chaotic (see, e.g. [11]). Moreover, the dynamics on the corresponding synchronous manifold of the coupled HR neurons may generate multistability region (see equation (10) and Table 1) such as chaotic attractor, stable periodic solutions and stable fixed point.

Neuronal synaptic connections are either chemical or electrical, and chemical connections might be excitatory or inhibitory. Moreover, the electrical coupling through gap junctions is bidirectional, whereas chemical synapses is unidirectional from the presynaptic cell to the postsynaptic cell. In

fact, the current qij injected from the presynaptic cell j to the postsynaptic cell i, is a nonlinear

function of the membrane potential xj of the presynaptic cell and a linear function of the membrane

potential xi of the postsynaptic cell. The current qij has the following form

qij = gs(v − xi)p(xj) (2a)

where gsis the strength of chemical coupling and v = 2 is the synaptic reversal potential. If v > xj,

the current injected to the cell is positive and depolarizes it, thus the coupling is excitatory. On the

other hand, for v < xj, the injected current to the cell is negative and consequently hyperpolarizes

it, thus introducing inhibitory coupling. The synaptic coupling function is modeled by the sigmoidal function

p(xj) =

1

1 + exp{−λ(xj− θs)}

, (2b)

where θs = −0.25 is the threshold and λ = 10. The threshold is chosen so that every spike in the

single neuron burst can reach the threshold. In the limit λ → ∞, the above sigmoid function reduces to a Heaviside step function.

We are now in a position to consider a network of n excitatory HR neurons with bidirectional electrical coupling and unidirectional chemical coupling. The equations of motion are the following. For, i = 1, . . . , n ˙ xi= f (xi) + yi− zi+ q + σ n X j=1 gijxj− gs(xi− v) n X j=1 cijp(xj), = f (xi) + yi− zi+ q + σ n X j=1 gijxj− gsk(xi− v)p(xi) − gs(xi− v) n X j=1 c†_ijp(xj), (3) ˙ yi= −yi− 5x2i + 1, ˙ zi= µ(b(x − x0) − zi), where cij = 0 or 1, cii= 0, n X j=1

and C†= (c†_ij), and c†_ij =    −k if i = j, cij if i 6= j. (4b) Here k represents the number of signals each neuron receives. Moreover, σ is the coupling strength

for electrical synapses via gap junctions, and coupling matrix G =: (gij), is a symmetric matrix with

vanishing row sums and nonnegative off-diagonal entries. It should be noted that the symmetry of

G is a biological assumption. From the mathematical side, our analysis here is capable of treating

unsymmetrical matrix with both positive and negative off-diagonal entries. C =: (cij) is the

connec-tion matrix of the chemical coupling which is not necessarily symmetric; cij = 1 if neuron i receives

synaptic current(via chemical synapses) from neuron j, otherwise cij = 0. The matrix C†=: (c†ij)

has all row sums being zero and nonnegative off-diagonal entries.

The notions of global and local synchronization are to be distinguished. Coupled network (3) synchronizes globally if for any solution we have

max{|xi(t) − xj(t)|, |yi(t) − yj(t)|, |zi(t) − zj(t)|} → 0 as t → ∞ for any i, j ∈ {1, 2, . . . , n}. (5)

Coupled network (3) is said to be locally synchronized if there exists an ε > 0 such that for any solution with

max{|xi(0) − xj(0)|, |yi(0) − yj(0)|, |zi(0) − zj(0)|} < ε,

we have (5) holds true.

To isolate the synchronous manifold, we need to make a coordinate transformation. Specifically, let xi = xi− x1, yi = yi− y1, zi= zi− z1, i = 2, ..., n. Then for i = 2, ..., n,

˙ xi = f (xi) − f (x1) + yi− zi+ σ n X j=1 (gij− g1j)xj − gsk [p(xi)(xi− x1) + (x1− v)(p(xi) − p(x1))] − gs(xi− v) n X j=1 (c†_{ij − c}†_1j)p(xj) − gs(xi− x1) n X j=1 c†_1jp(xj) = f′(si)xi+ yi− zi+ σ n X j=2 (gij− g1j)xj− gsk[p(xi)xi+ (x1− v)p′(ui)xi] − gs(xi− v) n X j=2 (c†_{ij − c}†_1j)p′(uj)xj− gsxi n X j=2 c†_1jp′(uj)xj. (6a)

The last equality above is justified by using the Mean Value Theorem and the fact that all row sums

of matrices G and C† are zero. Moreover,

˙

y_{i = −yi− 5(x}i+ x1)xi, (6b)

and

˙

The following notations are needed to set up the vector-matrix form of (6a-6c). Let x= (x1, . . . , xn)T, y = (y1, . . . , yn)T, z = (z1, . . . , zn)T, x= (x2, . . . , xn)T, y = (y2, . . . , yn)T, z = (z2, . . . , zn)T. Set e = √1_n(1, 1, ..., 1)T_, E1=         −1 1 0 _{· · · 0} −1 0 1 . .. ... .. . ... . .. ... 0 −1 0 · · · 0 1         (n−1)×n , A =   E1 eT   , and E† 1= ET1(E1ET1)−1. (7a) Then A−1_{= (E}†

1, e). For any matrix K ∈ Rn×nwhose row sums are all equal to zero, we have that

AKx= AKA−1Ax=   E1KE † 1 0 ∗ 0   Ax =   E1KE † 1 0 ∗ 0      x n P i=1 xi √ n    =   E1KE † 1x +   =:   Kx +   . (7b)

Therefore, in vector-matrix form, (6a-6c) become

˙x =nD1(t) + σG − gs[kD2(t) + k(x1− v)D3(t) + D4(t)C†D3(t) + D5(t)C∗D3(t)] o x_{+ y − z} =: H(t)x + y − z (8a) ˙y = −5D6(t)x − y, (8b) and ˙z = µbx − µz. (8c) Here, D1 = diag(f′(s2(t)), ..., f′(sn(t))), D2 = diag(p(x2(t)), ..., p(xn(t))), D3 = diag(p′(u2(t)), ..., p′(un(t))), D4= diag(x2(t) − v, ..., xn(t) − v), D5 = diag(x2(t) − x1(t), ..., xn(t) − x1(t)), D6= diag(x1(t) + x2(t), ..., x1(t) + xn(t)), (8d) C∗ =      c12 c13 . . . c1n .. . ... . .. ... c12 c13 . . . c1n      ,

Gand C† _{are defined as in (7b). From here on, an (n − 1) × n full-rank matrix with all its row sums}

being zero is to be termed as coordinate transformation. Matrix E1, as given in (7a), is an example

To study the global synchronization of (3) is then equivalent to showing that the origin of (8a-8c)

is globally, asymptotically stable. For local synchronization of (3), Di(t), i = 1, ..., 6, are reduced to

D1(t) = f′(x(t))I, D2(t) = p(x(t))I, D3(t) = p′(x(t))I,

D4(t) = (x(t) − v)I, D5(t) = 0, D6(t) = 2x(t)I (9)

where x(t) lies on the synchronous manifold of (3). Specifically, x(t), y(t) and z(t) satisfy the synchronous equation

˙x = f (x) + y − z + q − kgs(x − v)p(x),

˙y = −y − 5x2+ 1, (10)

˙z = µ(b(x − x0) − z).

In our derivation of synchronization of system (3), we need the concepts of matrix measures, a function being of type K, which generates a monotone dynamics, and a Lyapunov order number of the system of linear differential equations. For completeness and ease of the references, we also recall the definitions of the above described concepts and their properties [20, 19, 18].

Definition 1. _{([18]) Let || · ||}i be an induced matrix norm on Rn×n. The matrix measure of matrix

Kon Rn_{×nis defined to be µ}

i(K) = lim

ǫ_→0+

||I + ǫK||i− 1

ǫ .

Lemma 1. _{([18]) Let || · ||}k be an induced k-norm on Rn×n, where k = 1, 2, ∞. Then each of matrix

measure µk(K), k = 1, 2, ∞, of matrix K = (kij) on Rn×n is, respectively,

µ_∞(K) = max i {kii+ X j_6=i |kij|}, µ1(K) = max j {kjj+ X i6=j |kij|}, and µ2(K) = λmax(KT + K)/2.

Here λmax(K) is the maximum eigenvalue of K.

Theorem 1. ([18]) (i) µi(αA) = αµi(A), ∀ α ≥ 0 (ii) µi(A + B) ≤ µi(A) + µi(B). (iii) If λ is an

eigenvalue of A, then Reλ ≤ µ(A). (iv) Consider the differential equation ˙x(t) = K(t)x(t) + v(t),

t ≥ 0, where x(t) ∈ Rn_{, K(t) ∈ R}n×n_{, and K(t), v(t) are piecewise-continuous. Let || · ||i be a norm}

on Rn_{, and || · ||}

i, µi denote, respectively, the corresponding induced norm and matrix measure on

Rn×n. Then whenever t ≥ t0 ≥ 0, we have

||x(t)||i ≤ ||x(t0)||iexp Z t t0 µi(K(s))ds  + Z t t0 exp Z t s µi(K(τ ))dτ  ||v(s)||ids.

Let Rn

+ = {x = (x1, x2, . . . , xn)T ∈ Rn : xi ≥ 0, i = 1, . . . , n} be the nonnegative cone. Let

a_{, b ∈ R}n_{. We write a ≤ b if b − a ∈ R}n₊.

Definition 2. We say that a function f = (f1, . . . , fn) : D ⊂ Rn→ Rnis of type K on D if, for each

i, fi(a) ≤ fi(b) whenever a = (a1, . . . , an) and b = (b1, . . . , bn) are in D with a ≤ b and ai = bi.

The following theorem amounts to saying that a vector field being of type K is a sufficient condition to generate a monotone dynamics.

Theorem 2. ([19]) Let f (t, x) be of type K on Rn _{for each fixed t and let x(t) be a solution of}

˙x(t) = f (t, x) on [a, b]. Let z(t) be continuous on [a, b] and satisfy Dlz(t) ≤ f(t, z). Then z(t) ≤ x(t)

for a ≤ t ≤ b provided that z(a) ≤ x(a).

Consider a function y(t) for t ≥ 0. A number τ is called a Lyapunov order number [20] for y(t)

if, for every ǫ > 0, there exist positive constants c1,ǫ and c2,ǫ such that

||y(t)|| ≤ c1,ǫe(τ +ǫ)t for all large t,

||y(t)|| ≥ c2,ǫe(τ −ǫ)t for some arbitrary large t.

Consider linear system of differential equations in the homogeneous case

y′ = A(t)y. (11a)

Here A(t) is a d×d matrix. Clearly, the nontrivial solutions of (11a) have d Lyapunov order numbers

τ1, . . . , τd. Let τmax = max{τ1, . . . , τd}. Then τmax is called the Lyapunov order number of the

system. The Lyapunov order number of linear system of differential equations in the inhomogeneous case

y′ = A(t)y + f (t) (11b)

is to be defined similarly.

Proposition 1. _{([20]) (i) If y(t) 6= 0 for all large t, then the Lyapunov order number for y(t) is}

equal to lim

t→∞

ln ||y(t)||

t .

(ii) If A(t) is continuous for t ≥ 0, then a sufficient condition for every nontrivial solution y(t) of (11a) to possess an order number τ is that

Rt

0||A(s)||ds

t be bounded, in which case, τ ≤ limt→∞

Rt

0µ2(A(s))ds

t .

3

Main Results

We begin with the study of the bounded dissipation of coupled system (3), which is needed in proving the global synchronization of (3). The following coupled system, which is slightly more general than

that of in (3), is considered. For i = 1, . . . , n, ˙ xi = h1(xi) + b2yi+ b3zi+ σ n X j=1 gijxj− gs(xi− v) n X j=1 cijp(xj), ˙ yi = −c2yi+ h2(xi), (12a) ˙ zi = −c3zi+ h3(xi).

Suppose (x(t), y(t), z(t)) lies on the synchronous manifold of (12a), then it satisfies the synchronous equation

˙x = h1(x) + b2y + b3z − kgs(x − v)p(x),

˙y = −c2y + h2(x), (12b)

˙z = −c3z + h3(x).

Here k and p(x) are defined as in (4a) and (2b). Clearly, if (xc, yc, zc) is an equilibrium of (12b),

then gs=  c2c3h1(xc) + b2c3h2(xc) + b3c2h3(xc) c2c3k(xc− v)p(xc)  , (12c)

provided xc 6= v. To obtain the bounded dissipation of system (12a), we first prove that xc is

bounded for all gs ≥ 0.

Lemma 2. Let g1(x) = h1(x) + b_c2₂h2(x) +b_c3₃h3(x). Assume that g1(v) < 0, and lim

x→−∞g1(x) = ∞.

Then there exists an equilibrium (xc, yc, zc) of (12b) so that xc ∈ (d1, v) for some constant d1,

independent of gs≥ 0.

Proof. From the assumption on g1(x), there is a d1 with d1 < v such that g1(d1) > 0. Let g2(x) =

(x − v)p(x). Since g2(v) = 0 and g2(x) < 0 whenever x < v, there is an intersection point xc of

y = g1(x) and y = kgsg2(x) whose x-coordinate lies in (d1, v) for all gs ≥ 0. It implies g1(xc) −

kgs(xc− v)p(xc) = 0. Consequently, (xc,_c1₂h2(xc),_c1₃h3(xc)) is an equilibrium of (12a).

Theorem 3. Let (xc, yc, zc) be an equilibrium of (12b). Assume that

(i) There exists a k1 > 0 such that lim

x_→±∞ k14cjxh1(x) k2₁b2_jx2_{+ 2k} 1bjxhj(x) + h2_j(x) ≤ −1 and lim |x|→∞|k1bjx + hj(x)| = ∞ where j = 2, 3. (13a) (ii) c2, c3 > 0, (13b) and

(iii) xc is bounded by a constant, which is independent of gs≥ 0. (13c)

Suppose the matrix measure µ2(G) of G = (gij) is nonpositive and C = (cij) is a nonnegative matrix

in the componentwise sense. Let p(x) be a nonnegative, bounded function. Then the corresponding

Proof. Consider an energy function V of the form V (x, y, z) = k1

2

Pn

i=1(xi − xc)2 + 1₂Pni=1(yi−

yc)2+1₂Pn_i=1(zi− zc)2. Then the time derivative of V along (12a) has the following form

˙ V =k1 n X i=1 (xi− xc) ˙xi+ n X i=1 (yi− yc) ˙yi+ n X i=1 (zi− zc) ˙zi = − n X i=1 c2  yi− (h2(xi) + k1b2(xi− xc) + c2yc) 2c2 2 − n X i=1 c3  zi− (h3(xi) + k1b3(xi− xc) + c3zc) 2c3 2 + n X i=1  k1(xi− xc)h1(xi) − ych2(xi) − zch3(xi) + (h2(xi) + k1b2(xi− xc) + c2yc)2 4c2 + (h3(xi) + k1b3(xi− xc) + c3zc)2 4c3  + k1σxTGx− gsk1 n X i=1   (xi− xc)(xi− v) n X j=1 cijp(xj)    = : I1+ I2+ I3+ I4+ I5, where I4= k1σxTGx≤ 0. (14)

It follows from (13a-13c) that I3→ −∞ whenever kxk → ∞. Moreover,

I5 ≤ −gsk1 n X i=1   ( xc+ v 2 − xc)( xc+ v 2 − v) n X j=1 cijp(xj)    ≤ n(k1k)(v − xc )2 4  c2c3h1(xc) + b2c3h2(xc) + b3c2h3(xc) c2c3k(xc− v)p(xc)  = −nk1v − xc 4  c2c3h1(xc) + b2c3h2(xc) + b3c2h3(xc) c2c3p(xc)  . (15)

Therefore, if either k y k or k z k→ ∞, then I1+ I2 + I5 → −∞. We then conclude that if the

energy level of V is sufficiently large, then ˙V along the solution trajectory of (12a), whenever lying

in the outside the ellipsoid V (x, y, z) = c, is strictly negative for any σ ≥ 0 and any gs ≥ 0.

Tremendous progress (see e.g.,[21] and the work cited therein) has been made in the theory of

global synchronization of coupled chaotic systems, i.e., gs= 0. To obtain such global synchronization,

one needs to assume the bounded dissipation of the coupled chaotic systems. However, to the best of our knowledge, no general theorem until now is developed to check if the coupled systems are bounded dissipative. Moreover, our approach is quite general. For instance, if the individual oscillator is governed by the FitzHugh-Nagumo equation, then (13a-13b) are satisfied.

If the single neuron is governed by HR model, then c2 = 1, c3 = µ, b2 = 1, b3 = −1, h1(x) =

ax2_{− x}3_{+ q, h}

2(x) = −5x2+ 1 and h3(x) = µb(x − x0). For j = 2, it is easy to see that the first

limit in (13a) tends to −4k1

25 as x tends to ±∞. For j = 3, the first limit in (13a) tends to −∞ as x

tends to ±∞. The second limit in (13a) is obviously satisfied. The assumptions of Theorem 3 are

then satisfied. Hence, the coupled system is bounded dissipative for all σ ≥ 0 and for all gs≥ 0. It

also can be easily checked that if the single neuron is satisfied by FitzHugh-Nagumo equations, then the corresponding coupled system still enjoys the property of being bounded dissipative.

The following propositions, which, among other things, makes use of time average estimates, plays a critical step in obtaining our main results.

Proposition 2. _{Suppose ξ(t), η(t) and ζ(t) are nonnegative functions on [0, ∞) satisfying the} fol-lowing inequalities ξ(t) ≤ a0φ1(t, 0) + Z t 0 φ1(t, s)(a4(s)η(s) + a5(s)ζ(s))ds, (16a) η(t) ≤ b0φ2(t, 0) + Z t 0 φ2(t, s)(a6(s)ξ(s) + a7(s)ζ(s))ds, (16b) ζ(t) ≤ c0φ3(t, 0) + Z t 0 φ3(t, s)(a8(s)ξ(s) + a9(s)η(s))ds. (16c)

Here ai(t), i = 4, 5, ..., 9, are nonnegative functions on [0, ∞) and φi(t, s) = e

Rt

sai(τ )dτ, i = 1, 2, 3.

Then ξ(t), η(t), and ζ(t) converge to zero exponentially provided that lim

t→∞ Rt 0µ2(A(s))ds t ≤ −r, for some r > 0, where A(t) =      a1(t) a4(t) a5(t) a6(t) a2(t) a7(t) a8(t) a9(t) a3(t)      , (17)

Proof. Let ξ(t), η(t) and ζ(t) satisfy the following equation.

˙ξ = a1(t)ξ + a4(t)η + a5(t)ζ, ξ(0) = ξ(0),

˙η = a6(t)ξ + a2(t)η + a7(t)ζ, η(0) = η(0),

˙ζ = a8(t)ξ + a9(t)η + a3(t)ζ, ζ(0) = ζ(0).

It is easily checked that the above system is of type K. Hence, ξ(t) ≥ ξ(t), η(t) ≥ η(t), and ζ(t) ≥ ζ(t), ∀t ≥ 0. Using Proposition 1-(ii), and Theorem 2, we see that the proposition holds as claimed.

Proposition 3. _{(i) Let ˙x = A(t)x. Here A(t) is an n × n matrix. Suppose lim}

t→∞

A(t) = A. Then

x(t) converges to zero exponentially provided that µ2(A) < 0. (ii) If, in addition, the eigenvalues of

A are distinct and their real parts are negative, then x(t) converges to zero exponentially.

Proof. To prove the first assertion of the proposition, we see that A(t) = A + (A(t) − A), and so,

µ2(A(t)) ≤ µ2(A)+µ2(A(t)−A) = µ(A)+λmax

h (A(t)−A)+(A(t)−A)T 2 i ≤ µ2(A)+µ1 h (A(t)−A)+(A(t)−A)T 2 i We have used Theorem 1-(iii) to justify the last inequality. To complete the proof of the first part of the proposition, one needs to show that lim

t→∞

Rt

0µ1((A(s)−A)+(A(s)−A)T)ds

t ≤ 0, which is the case

provided that lim

t→∞

Rt

0|aij(s)−aij|ds

t = 0, for all i. To see this, let ε > 0 be given, there exists a tε such

that |aij(s) − aij| < ε whenever t ≥ tε. And so,

Rt 0 |aij(s) − aij|ds ≤ Rtε 0 |aij(s) − aij|ds + ε(t − tε). Hence, lim t→∞ Rt 0|aij(s)−aij|ds

t = 0. The first assertion of the proposition now follows. To conclude the

diagonal matrix. Thus P      ˙ξ ˙η ˙ζ     = PAP −1_{+ P(A(t) − A)P}−1
_P      ξ η ζ      . Now, lim t→∞ Rt

0µ2(PAP−1+ P(A(s) − A)P−1)ds

t ≤ λmax+ limt→∞

Rt

0 µ2(P(A(s) − A)P−1)ds

t = λmax,

where λmax is the largest real part of the eigenvalues of A. We have used the similar techniques,

as given in the proof of the first part of the proposition, to justify the above equality. We have completed the proof of proposition.

In view of (8a), it is apparent that the more negative the matrix measure of H(t) is, the easier the origin of the system (8a-8c) is to be made asymptotically stable. However, the choice of a coordinate

transformation [21] will greatly influence how negative the matrix measure of G and C†, as appeared

in (8a), could be. In the earlier works, the choice of coordinate transformations is either E1 or

        1 −1 0 _{· · ·} 0 0 1 ₋₁ . .. ... .. . . .. ... ... 0 0 · · · 0 1 ₋₁         .

The drawback for the above choices is that even if G is the diffusive matrix with periodic boundary conditions, the corresponding matrix measure of G is positive whenever n > 7 (see Table 1 of [21]). Therefore, a better choice of transformation is needed. Let

C_{= {C ∈ R}(n−1)×n _{: C is full rank, and all its row sums are zeros}.}

and O ⊂ C be such that

O₌   C∈ C :   C eT   is orthogonal   .

Proposition 4. The following hold for any E1 and E2 ∈ O.

(i) E2E†1E1 = E2. (ii) E1E†2E2 = E1. (iii) (E2E†1)G(E2E†1)−1 = E2GE†2 =: eG. (iv) σ( eG) =

σ(G) = σ(G) − {0}, where σ(G) is the spectrum of G. Proof. We have, via (7a), that the inverse of A =

  E1 eT   is A−1 ₌  _E† 1, e  . Hence I = A−1A =  E†₁, e   E1 eT   = E†

proposition now follows since the row sums of E2 are all zeros. The proof of the second

asser-tion of the proposiasser-tion is similar, and thus omitted. To complete the proof of the third asserasser-tion

of the proposition, we first note, via (i), that (E2E†1)−1 = E1E†2. Hence (E2E†1)G(E2E†1)−1 =

E2E†1E1GE†1E1E2†= E2G(I − eeT)E†2 = E2GE†2. The last assertion of the proposition follows from

the fact that   E1 eT   G   E1 eT   −1 =   E1 eT   G E†₁, e  =   G 0 ∗ 0   , where E1 ∈ O.

Proposition 5. ([21]) Assume that all eigenvalues of the coupling matrix G have nonpositive real

parts. Then infC_∈Cµ2(CGC†) ≥ Re λ2(G). Here Re λ2(G) is the second largest real part of

eigen-values of G. If, in addition, G is symmetric, then the above equality can be achieved by choosing any C in O.

Proposition 6. ([21]) Assume that G is a node-balancing matrix, i.e., its row sums and column

sums are equal. Then

µ2(G) = λ2(

G_{+ G}T

2 ), (18)

whenever all eigenvalues of G + GT _{are nonpositive.}

We are now in a position to make a better coordination transformation. Letting

e x= E2      x1 .. . xn     , ey= E2      y1 .. . yn     , ez = E2      z1 .. . zn     ,

where E2∈ D, multiplying E2E†1 on both sides of equation (8a), we get

˙ex = E2E†1H(t)(E2E1†)−1ex+ey− ez.

We have used Proposition 4-(ii) to justify the above equation. Applying Proposition 4-(iii), we have that

˙ex = eH(t)_ex+y_{e− ez,} (19a)

where e H(t) =  e D1(t) + σ eG+ gs[−k eD2(t) − k(x1(t) − v) eD3(t) − eD4(t) fC†De3(t) − eD5(t) fC∗De3(t)]  (19b) =:nDe1(t) + σ eG+ gsHe1(t) o .

Here eDi(t) = E2E†1Di(t)(E2E†1)−1, and all other terms on the right-hand side of (19c) with e on

measure of eG remains negative regardless the number of neurons involved. The dynamics of the

motion fory_e and _{ez is, respectively,}

˙ey = −5 eD6(t)ex− ey, (19c)

and

˙ez = µbex − µez. (19d)

Theorem 4. (i) Assume that µ2( eH(t)) ≤ λ < 0 and lim sup{5k eD6(t)k} = d. If −λ > d + b, then

the origin is globally, asymptotically stable with respect to equations (8a-8c).

(ii) Assume that µ2(G) is nonpositive. Let λ2 be the second largest eigenvalue of 1₂(G + GT). Let

d1 = lim sup µ2( eD1(t)) and h1 = lim sup µ2( eH1(t)). Then the coupled HR system (3) is globally

synchronized provided that

(−λ2)σ + (−h1)gs> d + b + d1. (20)

Proof. Using the inequality in Theorem 1-(iv), we get kex_{(t)k ≤ ke}x0keλt+

Z t

0

eλ(t−s)_(key_{(s)k + kez(s)k)ds.} (21a) Moreover key_{(t)k ≤ ke}y0ke−t+ Z t 0 e−(t−s)_(dkex_(s)k)ds, (21b) and kez(t)k ≤ kez0ke−µt+ Z t 0 e−µ(t−s)_(µbkex_(s)k)ds. (21c) Applying Proposition 2, we see that the first part of assertion of Theorem 4 holds true provided the real parts of the eigenvalues of

B=      λ 1 1 d ₋₁ 0 µb 0 _−µ      (22)

are negative. Now, the characteristic polynomial q(x) of B is q(x) = x3 _{+ a}

2x2+ a1x + a0, where

a2= −λ + µ + 1, a1 = −λ − d − µb − λµ + µ and a0 = −λµ − dµ − bµ. The Routh-Hurwitz Criterion

asserts that all roots of q(x) have negative reals parts if and only if a0 > 0, a2 > 0, and a1a2 > a0.

A direct computation yields that a1a2> a0 provided that a0> 0, which is equivalent to −λ > d + b.

The first part of the theorem thus follows.

From (19b), we have µ2( eH(t)) ≤ µ2( eD1(t)) + σµ2( eG) + gsµ2( eH1(t)) ≤ d1+ λ2σ + h1gs. Hence,

Theorem 5. Let x(t) satisfies synchronous equation (10). Let r2 be the second largest eigenvalue of C†, as defined in (4b). Let α =: 1 + r2 k and hkgs,α(x) = f′(x) + kgs  −p(x) + (v − x)p′(x)α. (23a) Set sup x hkgs,α(x) =:    (hα)kgs if kgs6= 0, d1 if kgs= 0, (23b)

where hα is a constant and d1 = max

x∈R f

′_{(x) ≈ 2.253. Let d = sup}

x(t)10|x(t)| ≤ 20. Then coupled HR

system (3) is locally synchronized provided that

(−λ2)σ + (−hα)kgs> d + b = 20 + b, for kgs> 0, and

−λ2σ > d + b + d1 ≈ 26.253, for kgs= 0. (23c)

Proof. Since we consider excitatory HR neurons, x(t) < v = 2 for all t. From (9) and (19b), we then

have µ2( eH(t)) ≤ λ2σ + (−hα)gsk. By Proposition 2, the assertion of the theorem holds true.

If the network is globally coupled, then k = n − 1, r2 = −n, and so α is negative. It can be

computed that the denser the network is coupled, the smaller α is. Hence, α is an indicator of how densely coupled the system is. Note also that −1 ≤ α < 1. We shall call α the density of the coupling network.

4

Applications

In the section, we shall focus on obtaining synchronization theories for coupled HR neurons (3). To this end, one needs to have information on the dynamics of synchronous equation (10). Its dynamics is to be provided numerically. The theoretical study of (10) is to be addressed in a forthcoming paper.

For q = 4, a = 2.6, x0= −1.6, µ = 0.01, b = 4 and v = 2, the single HR neuron model, i.e., gs= 0, is

capable of producing major neuronal behaviors such as spiking, bursting, and chaotic regimes. (see

e.g., [14]). Furthermore, such neuron is excitatory, i.e., x(t) < v = 2 for all t ≥ 0. We shall treat kgs

as a bifurcation parameter. The corresponding dynamical behavior of equation (10) is summarized in Table 1. On the synchronous manifold, the solution trajectory x(t) of (10), depending on initial

conditions and kgs, may settle into various stable states. Figure 1 provides the maximum Lyapunov

exponent (MLE) of synchronous equation (10) versus kgs. There is a set of initial conditions with

positive measure for which their corresponding MLE is positive for 0 ≤ kgs ≤ 0.84. On the same

range of parameters kgs, there is also a set of initial conditions with positive measure for which its

corresponding MLE is negative. For instance, if 0.809 ≤ kgs ≤ 0.813, then there are sets of initial

solution (see Fig.2) and chaotic attractor (see Fig.3), respectively. Specifically, let (xc, yc, zc) be

the steady state of (10), and let Cr = {(x, y, z) : |x − xc| < r, |y − yc| < r, and |z − zc| < r},

and Ir = {(x, y, z) : |x − xc| > r, |y − yc| > r, and |z − zc| > r}. In fact, our numerical results

suggest that the following hold. Pick, for instance, kgs = 0.812. If the initial condition (x0, y0, z0)

is randomly chosen form C0.02 (resp., I1), then its trajectory converges to a periodic orbit (resp., a

chaotic attractor) (see Figs. 2,3). Similarly, for kgs ∈ [0.814, 0.84], synchronous equation (10) also

exhibits rich dynamics showing the coexistence of stable multi states. Moreover, if kgs ≥ 0.814,

the numerical results suggest that the corresponding steady state is locally stable. If one performs

linearized stability at the steady-state (xc, yc, zc), then one sees that (xc, yc, zc) is stable whenever

kgs≥ 0.814 (see Fig.4). Such analysis of linearized stability provided some supportive evidence for

the validity of Table 1.

gsk gsk <0.808 0.809 ≤ gsk≤ 0.813 0.814 ≤ gsk≤ 0.84 gsk≥ 0.87

Chaotic attractor Chaotic attractor

Types Chaotic attractor Stable steady state Stable periodic solution Stable steady state

Table 1: The dynamics of synchronous equation (10) with various range of kgs. The

multi-stability of (10) is observed with kgs∈ [0.809, 0.84].

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 0.01 0.02 0.03 0.04 X: 0.84 Y: 0.0014 k⋅g_s M L E

Figure 1: The maximum Lyapunov exponent (MLE) of synchronous equation (10) is com-puted for various kgs. For 0 ≤ kgs ≤ 0.84, MLE> 0 for a set of initial conditions with

positive measure. For kgs ∈ [0.809, 0.84], there is also a set of initial conditions for which its

corresponding MLE< 0.

In summary, the numerical results suggest that on the synchronous manifold, for kgs small, the

chaotic behavior of single HR persists. For kgs in an intermediate range, the multistability of

equa-tion (10) occurs. Depending on initial condiequa-tions, the coexistence of multi stability states including a chaotic attractor and a stable periodic solution/a stable fixed point could be observed. When

kgs becomes large, equation (10) has a globally asymptotically stable fixed point. Such complex

dynamical behavior of synchronous equation (10) leads to the possibility of stable multi-state

0.01 0.02 0.03 0.04 0.05 0.99 0.995 1 1.005 6.502 6.504 6.506 6.508 6.51 x y z (x(t),y(t),z(t))

Figure 2: The solution trajectory with randomly chosen initial conditions form C0.02converges

to a stable periodic orbit. Here kgs= 0.812.

−3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −35 −30 −25 −20 −15 −10 −5 0 5 2 3 4 5 6 7 8 x y z (x(t),y(t),z(t))

Figure 3: The solution trajectory with randomly chosen initial conditions form I1 converges

to a stable chaotic attractor. Here kgs = 0.812.

that the corresponding synchronous equation leads to a chaotic solution, then the associated coupled HR neurons (3) achieves stably chaotic synchronization. Likewise, we define stably periodic syn-chronization and stably steady-state synsyn-chronization accordingly. As we can see, via Table 1, that for 0.809 ≤ k ≤ 0.84, the coexistence of stable multi-state synchronization of coupled HR neurons (3) could be observed.

(I) Local synchronization

We next turn our attentions to the local synchronization of coupled HR neurons, from which much more information can be extracted.

0.5 0.6 0.7 0.813 0.9 −0.5 0 1 2 3 3.5 λm a x ( ¯ A) kgs

Figure 4: The maximum eigenvalue of the linearized operator with respect to the synchronous equation (10) is computed for various kgs.

10 20 30 40 0 0.5 0.6044 0.87 1 1.5 −λ2σ k gs Sync. Region 26.253

Figure 5: The shaded area is the synchronization region satisfied by (25) and kgs≥ 0.87.

(a) Neurons with only chemical synapse.

In [17], a local steady-state synchronization of bursting neurons with no electrical coupling is studied without providing mathematical details. Moreover, their approach fails to see if synchro-nization of neurons can be achieved when the networks are intermediately and sparsely coupled. Their major contribution was to prove that the bound for achieving synchronization of HR neurons depends only on the number k of signals each neuron received, and is independent of all other de-tails of the network topology. Form (23a) and (23c), it is clear that the smaller the density α of the network is, the greater chance coupled HR neurons (3) gets synchronized. In the following, we shall prove that the system of coupled HR neurons (3) achieves steady-state synchronization regardless

0 100 200 300 400 500 600 700 1000 2000 −0.0005 0.0005 0.002 t x1 −x 2 0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 0.02 0.04 0.06 t x1

Figure 6: The time series of x1(t)−x2(t) and x1(t) The graphs demonstrate the stable periodic

synchronization. Here σ = 30, kgs= 0.812, initial = [0.26459e − 1 + r, 0.996499 + r, 6.5058 +

r,_{0.26459e − 1 − r, 0.996499 − r, 6.5058 − r] and r = 0.001.} 0 50 100 150 200 250 300 350 400 1000 2000 −0.2 0 0.2 1.945 t x1 −x 2 0 200 400 600 800 1000 1200 1400 1600 1800 2000 −4 −2 0 2 t x1

Figure 7: The time series of x1(t)−x2(t) and x1(t). The graphs demonstrate the stable chaotic

synchronization. Here σ = 30, kgs = 0.812, initial= [0.26459e − 1 + r, 0.996499 + r, 6.5058 +

r,_{0.26459e − 1 − r, 0.996499 − r, 6.5058 − r] and r = 1.}

To this end, we begin with the study of equations (8a-8c) and (9) with σ = 0. A coordination transformation, such as the one performed in the proof of Proposition 3-(ii), is applied on those equations, their resulting equations can be written as follows.

d dt      ξi ηi ζi     = Ai(t)      ξi ηi ζi     , where Ai(t) =      f′_{(x(t)) + hi(x(t)) 1 −1} −10x(t) −1 0 µb 0 _−µ     , i = 2, ..., n.

0 0.2 0.4 0.6 0.67 0.8 1 0.5 0.604 0.87 1 α k gs

Figure 8: Turning points of hα.

Here hi(x(t)) = kgs



−p(x(t)) + (v − x(t))p′(x(t))(1 +ri

k))



, where ri, i = 2, ..., n, are eigenvalues

of C† with r2 ≥ r3 ≥ ... ≥ rn. The problem of showing synchronization of system (8a-8c) and (9)

with σ = 0 is then equivalent to proving that (ξi, ηi, ζi), i = 2, ..., n, converge to zero. Upon using

Proposition 3, we conclude that for kgs≥ 0.814, (x(t), y(t), z(t)) converges to zero with a certain set

of initial conditions with positive measure provided that λmax(A2) < 0, where

A2 =      f′_{(xc) + h2(xc) 1 −1} −10xc −1 0 µb 0 _−µ      .

Note that A2is also the linearized matrix of (10) at (xc, yc, zc). Form Fig.4, we see that λmax(A2) < 0

whenever kgs≥ 0.814. Consequently, we have the following conclusion.

Coupled HR neurons (3) achieves the steady-state local synchronization whenever kgs ≥ 0.87

regardless how sparsely the network is coupled and only depends on the number of signals each neu-ron received. We are unable to prove that the existence of chaotic synchneu-ronization. However, for

these ranges of kgs we can show, as demonstrated earlier, that coupled neurons (3) may achieve

the steady-state synchronization on a set of initial conditions near (xc, yc, zc). Numerically, we also

observed that coupled two HR neurons achieves synchronization, whenever kgs ≥ 0.809. No

syn-chronization occurs for kgs ≤ 0.808. For kgs ∈ [0.809, 0.84], chaotic synchronization can only be

detected whenever initial conditions of two neurons are identical. Otherwise, for kgs ∈ [0.809, 0.813]

(resp., [0.814, 0.84]) stably periodic (resp., steady-state) synchronization is discovered whenever ini-tial conditions of two neurons are distinct. It should be noted that without the presence of electrical synapses the generic multi-state synchronization is not possible.

(b) Neural synchronization with only electrical synapse.

For gs = 0, we obtain stable chaotic synchronization with σ > d+b+d_−λ2 1 =: σmin. Consider, for

instance, a ring of 2l-nearest neighbor mutually coupled networks, the predicted minimum electrical

coupling strength σmin is computed with the number n of neurons and l being given. The results

are listed in Table 2. Note that in such case

λ2 = −4 l X i=1 sin2 iπ n = sin (l + 1₂)2π_{n
− sin}π_n sinπ n − 2l. (24) n 11 102_{+ 1 10}3_{+ 1 10}4_{+ 1}

l = 1, the nearest-neighbor coupling. 82.69 6790 _{6.67 × 10}5 _{6.66 × 10}7 l= [n−1

4 ] 17.66 1.4 0.144 0.0144

l = n−1

2 , the globally coupled network. 2.4 0.26 0.0262 0.0026

Table 2: The table gives the minimum electrical coupling strength σmin. For instance, with

n= 102_{+ 1, l =}n−1 4



= 25, the predicted minimum electrical coupling strength is σmin = 1.4.

(c) Neural synchronization with both electrical and chemical synapses.

In this subsection, networks with both electrical and chemical connections are considered. To

this end, we first observe that r2≤ 0, x(t) < v and p′(x(t)) ≥ 0 for all t. So hα, defined in (23b), is

increasing in α. If

(−λ2)σ + (−h1)kgs> d + b, (25)

then (23c) is also satisfied. The synchronization region satisfying (25) and kgs ≥ 0.87 is

demon-strated in Fig.5. That is to say, if (λ2σ, kgs) is chosen from the shaded region in Fig. 5, then

multi state or single state synchronization can be realized depending on the range of kgs and

ini-tial conditions. Consider, for instance, coupled two HR neurons. Let i = 1, 2 and kgs = 0.812. If

(xi0, yi0, zi0) ∈ C0.02(resp., I1) and are distinct, then the stably periodic (resp., chaotic)

synchroniza-tion occurs (see Figs. 6,7). We further observe that there exists a t1such that h1< 0 (resp., h1> 0)

whenever kgs∈ [0, t1) =: I1 (resp., kgs ∈ (t1, 0.87] =: I2) (see Fig.5). Here t1 ≈ 0.644. Hence, both

chemical and electrical synapses enforce the synchronization phenomena whenever kgs ∈ I1. For

kgs∈ I2, the chemical synapses play a dragging role to reach synchrony. To synchronize, the

electri-cal synapses have to be strong enough to suppress the dragging force created by chemielectri-cal synapses.

Such t1 is called a turning point for h1. We are then led to compute turning points for hα (see

Fig 8). For α ≤ 0.67 the corresponding turning points are kgs = 0.87. Hence, for 0 ≤ kgs < 0.87,

if the density of the coupling network is at least 0.67, then chemical synapses can also enforce the

predicted σ, σmin, so that local synchronization can be achieved. This, in turns, gives that a

sta-ble multi-state synchronization, including stasta-ble chaotic, periodic state/steady-state, or multi-state synchronization can be achieved.

To summarize, a synchronization region is obtained in Fig. 5. Particularly, multi-state

syn-chronization of coupled HR neurons can be realized whenever kgs ∈ [0.809, 0.84] and (−λ2σ, kgs)

lies in the synchronization region. Numerically, it should be mentioned that no generic multi state synchronization can be achieved without the presence of electrical synapses. It is also obvious that the system of the coupled neurons with only electrical synapses is incapable of producing multi-state synchronization. Hence, multi-state synchronization can only be achieved with both presence of

elec-trical and chemical synapses. Furthermore, for 0 ≤ kgs< 0.87, if the density of the coupling network

is at least 0.67, then chemical synapses can enforce the synchrony of the system. (II) Global synchronization

To study global synchronization of HR neurons, we shall explicitly find a rectangle region where coupled system (3) is bounded dissipation.

Proposition 7. Consider HR neurons (3) with a = 2.6, b = 4, q = 4, µ = 0.01 and x0 = −1.6.

Let (xc, yc, zc) be the equilibrium of (10). Let the energy function V be as the form considered in

Theorem 3 with k1 = 7. Then ˙V (x, y, z) < 0 whenever (x, y, z) lies in the outside of the rectangle

region

S = {(x, y, z) : |xi| ≤ m1, |yi| ≤ m2, |zi| ≤ m3, ∀i = 1, . . . , n}

where m1 = 53√4n, m2 = ₂5m21+72m1+ 16 + 1375

√_{n, and m}

3 = 348m1+ 711 + 1375√n. Moreover,

xi, ∀i = 1, . . . , n, will eventually stay in the ball Bec(0), where the radius ec is given in the following

ec =

√

3n · (max{m1, m2, m3} + 19) + 2.

Proof. It can be seen that coupled HR neurons (3) is of the form of (12a) with b2 = 1, b3 = −1,

c2 = 1, c3 = µ, h1(x) = ax2− x3+ q, h2(x) = 1 − 5x2, and h3(x) = µb(x − x0). Let the energy

function V be the same as in the proof of Theorem 3 with k1 = 7, i.e.,

V (x, y, z) = 7 2 n X i=1 (xi− xc)2+ 1 2 n X i=1 (yi− yc)2+ 1 2 n X i=1 (zi− zc)2.

Then from equations (14), (15), we have ˙ V ≤ − n X i=1  yi−1 − 5x 2 i + 7(xi− xc) + yc 2 2 − n X i=1 µ  zi− µb(xi− x0) − 7(xi− xc) + µzc 2µ 2 + n X i=1  7(xi− xc) · (ax2i − x3i + q) − yc(1 − 5x2i) − zcµb(xi− x0) +(1 − 5x 2 i + 7(xi− xc) + yc)2 4 +(µb(xi− x0) − 7(xi− xc) + µzc) 2 4µ  + n X i=1 7−(ax 2 c − x3c + yc− zc+ q) 4p(xc) (v − xc ) =: n X i=1 J1(xi, yi) + n X i=1 J2(xi, zi) + n X i=1 J3(xi) + n X i=1 J4.

Since xc ∈ [−1₂, 2] for all gs ≥ 0, it can be computed directly that J4 ≤ 20. Now J3(x),

J3(x) = − 3 4x 4₊ 7 10+ 7xc  x3+  122079 100 − 7 10xc+ 5 2yc  x2 +  2307 250 − 4921 2 xc+ 7 2yc− 88 25zc  x + " −28xc− yc− 8 125zc+ 1 4(1 − 7xc+ yc) 2_{+ 25}  8 125 + 7xc+ 1 100zc 2# < −1₄x4_{+ 1890000 ≤ −}1 4(x 4_{− 53}4_).

Hence, ˙_{V < 0, whenever |x}i| ≥ 53√4n =: m1 for some i = 1, . . . , n. If |yi| ≥ 5₂m21+ 72m1+ 16 +

1375√n =: m2, and |xi| ≤ m1 for some i = 1, . . . , n, then J1(xi, yi) ≤ −1890020n. Consequently,

˙

V < 0. Hence, if |yi| ≥ m2, for some i = 1, . . . , n, then ˙V < 0. Similarly, if |zi| ≥ 348m1+ 711 +

1375√n =: m3, for some i = 1, . . . , n, then ˙V < 0. Thus, ˙V (x, y, z) < 0 whenever (x, y, z) lies in

the outside of the rectangle region

S = {(x, y, z) : |xi| ≤ m1, |yi| ≤ m2, |zi| ≤ m3, ∀i = 1, . . . , n}.

Moreover, it can be proved easily that the region of V ≤ Vc, where

Vc =

21n

2 (max{m1, m2, m3} + 19)

2_,

contains the rectangle region S. It then follows that |xi(t)| ≤ ec for i = 1, . . . , n and large t. where

ec =

√

3n · (max{m1, m2, m3} + 19) + 2.

We have just completed the proof of the proposition.

Using Proposition 4-(iii), we see that

max µ2( eD1(t)) ≤ max k eD1(t) k≤ e max

x f ′_{(x) =} 169 75 e, max µ2(− eD2(t)) ≤ e max x (−p(x)) ≤ 0. µ2  (−k(x1(t) − v)) eD3(t)  ≤ ek(ec+ v)ep, µ2  − eD4(t) fC†De3(t)  = µ2  (E2E†1)(−D4C†D3)(E2E1)−1  ≤ e(ec+ v)ep k C†k, µ2  − eD5(t) fC∗De3(t)  ≤ 2eecep k C∗ k . where ep = max x p

′_{(x). Combining the above, we arrive at the following conclusion.}

HR neurons (3) achieves global synchronization provided that

(−λ2)σ −



ek(ec + v)ep k C†k +2eecepr3k C∗k



gs> d + b + d1. (26)

5

Conclusions

Synchronization of coupled HR neurons over complex networks with chemical and electrical synapses is analytical studied. Particularly, multi-state and multi-stage synchronizations are observed with the presence of both of chemical and electrical synapses. A measurement for the density of the network is introduced to ensure that chemical synapses play a positive effect on the synchronization of the system of coupled neurons. We conclude this work by mentioning the possible future work. It would be interesting to analytically study the rich dynamical behavior of synchronous equation (10). Numerically, one sees that coupled HR neurons are capable of producing multi-stage synchronization even without the help of electrical synapses. It is worthwhile to give a rigorous proof. The inequality

(26) is unsatisfactory since the coefficient of gs in (26) is negative. This inequality suggests that the

chemical synapses play a dragging role in the process of reaching synchrony. And only the electrical synapses play positive effect on the synchronization of the system. It is certainly an interesting task

to overcome some technical difficulties to obtain an inequality of (26) type with the coefficient of gs

being positive.

References

[1] W. Singer, Neuronal synchrony: a versatile code for the definition of relations?, Neuron 24, 49 (1999).

[2] P. J. Uhlhaas and W. Singer, Neural synchrony in brain disorders: Relevance for cognitive dysfunctions and pathophysiology, Neuron 52, 155 (2006).

[3] S. Neuenschwander and W. Singer, Long-range synchronization of oscillatory light responses in the cat retina and lateral geniculate nucleus, Nature (London) 379, 728 (1996).

[4] L. Glass, Synchronization and rhythmic processes in physiology, Nature (London) 410, 277 (2001).

[5] C. M. Gray, P. Konig, A. K. Engel, and W. Singer, Oscillatory responses in cat visual cortex ex-hibit inter-columnar synchronization which reflects global stimulus properties, Nature (London) 338, 334 (1989).

[6] R. Dzakpasu and M. ˙Zochowski, Changes in synchrony and phase synchrony between individual

neurons in normal and epileptic brain, Physica D 208, 115 (2005).

[7] C. J. Stam, B. F. Jones, G. Nolte, M. Breakspear, and Ph. Scheltens, Small-world networks and functional connectivity in Alzheimer’s disease, Cereb. Cortex 17, 92 (2006).

[8] G. Buzs´aki, Rhythms of the brain, Oxford University Press, N.Y. (2006)

[9] 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 17, 043128 (2007).

[10] E. M. Izhikevich, Which model to use for cortical spiking neurons? IEEE Trans. Neural Netw. 15, 1063 (2004).

[11] M. Storace, D. Linaro, and E.d Lange, The Hindmarsh-Rose neuron model: Bifurcation analysis and piecewise-linear approximations, Chaos 18, 033128 (2008).

[12] E. de Lange and M. Hasler, Predicting single spikes and spike patterns with the Hindmarsh-Rose model, Biol. Cybern. 99, 349-360 (2008).

[13] J. L. Hindmarsh and R. M. Rose, A model of neuronal bursting using three coupled first order differential equations, Proc. R. Soc. Lond. B 221, 87-102 (1984).

[14] G. Innocenti and R. Genesio, On the dynamics of chaotic spiking-bursting transition in the Hindmarsh-Rose neuron, Chaos 19, 023124 (2009).

[15] M. Jalili, Synchronizing Hindmarsh-Rose neurons over Newman-Watts networks, Chaos 19, 033103 (2009).

[16] N. Kopell and B. Ermentrout, Chemical and electrical synapses perform complementary roles in the synchronization of interneuronal networks, PNAS vol.101, no.43 15482-15487(2004). [17] I Belykh, E.d. Lange, and M. Hasler, What matters in the Network Topology, PRL 94,

188101(2005).

[18] M. Vidyasagar, Nonlinear Systems Analysis, Prentice Hall Inter., Inc, 1st ed. (1978).

[19] S. B. Hsu, Ordinary Differential Equations with Applications, Hackensack, NJ, London: World Scientific (2006).

[20] P. Hartman, Ordinary Differential Equations, Birkh¨anser, Boston, 2nd ed. (1982).

[21] J. Juang and Y. H. Liang, Coordinate transformation and matrix measure approach for syn-chronization of complex networks, Chaos 19, 033131 (2009).