3.2 Synchronization in Hindmarsh-Rose Neurons with Chemical and Elec-
3.2.2 Synchronization of the Hindmarsh-Rose Neurons
The HR model was obtained by biological consideration over the response to stimuli of a real neuronal cell. The motion of the model reads as follows:
˙x = f (x) + y − z + q,
˙y = −y − 5x2+ 1, (3.12)
˙z = µ(b(x − x0) − z).
Here f (x) = ax2− x3, 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 between 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 adaptation; 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 regular bursting (see, e.g., [87]).
Moreover, the dynamics on the corresponding synchronous manifold of the coupled HR neurons may generate multistability region (see equation (3.16) and Table 3.1) containing a stable regular bursting, a stable periodic solution and a stable fixed point.
Neuronal synaptic connections are either chemical or electrical, and chemical con-nections might be excitatory or inhibitory. Moreover, the electrical coupling through gap junctions is bidirectional, whereas the chemical synapse 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), (3.13a) where gs is the strength of chemical coupling and v is the synaptic reversal potential. If xi < v, the current injected to the cell is positive and depolarizes it, thus the coupling is excitatory. On the other hand, for xi > v, the injected current to the cell is negative and consequently hyperpolarizes it, thus introducing inhibitory coupling. In this thesis, we numerically choose v = 2 so that xi(t) < v for all t, thus the synapse is depolarizing (excitatory). It is certainly interesting to justify that such choice of v is always possible.
The chemically synaptic coupling function is modeled by the sigmoidal function
p(xj) = 1
1 + exp{−λ(xj− θs)}, (3.13b) 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 m excitatory HR neurons with bidirectional electrical coupling and unidirectional excitatory chemical coupling. The equations of motion are the following. For, i = 1, . . . , m,
Here k represents the number of chemical signals each neuron receives. Moreover, d is the coupling strength for electrical synapses via gap junctions, and coupling matrix G 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 matrices with both positive and negative off-diagonal entries. C is the connection 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 S has all row sums being zero and nonnegative off-diagonal entries.
We next describe the synchronous equation of HR network (3.14). On the
syn-chronous manifold, its dynamics is governed by the following equations
˙x = f (x) + y − z + q − kgs(x − v)p(x),
˙y = −y − 5x2+ 1, (3.16)
˙z = µ(b(x − x0) − z).
To study local synchronization, we begin with the derivation of the variational equation of (3.14) along the synchronous manifold M. The equation is
˙ui = f′(x(t))ui+ vi− wi+
"
d Xm
j=1
gijuj
#
− kgsp(x(t))ui
−
"
gs(x(t) − v)p′(x(t)) Xm j=1
cijuj
# ,
˙vi = −vi− 10x(t)ui, (3.17)
˙
wi = µ(bui− wi),
where x(t) lies on the synchronous manifold of (3.14) and satisfies equation (3.16).
In vector-matrix form, (3.17) becomes
˙
u= {[f′(x(t)) − kgsp(x(t))] I + dG
− gs(x(t) − v)p′(x(t))C}u + v − w,
= {[f′(x(t)) − kgsp(x(t)) − kgs(x(t) − v)p′(x(t))] I
+ dG − gs(x(t) − v)p′(x(t))S}u + v − w, (3.18a)
˙v = −10x(t)u − v, (3.18b)
˙
w = µbu − µw. (3.18c)
To study synchronized HR neurons (3.14), we first apply a coordinate transforma-tion, developing in Subsection 2.3.1, on G so that the resulting matrix has a negative matrix measure as possible. The structure of linear system (3.18) is then explored so that the theory of some monotone dynamics and time averaging estimates described in Section 2.1 can be applied to make the linear system asymptotically stable.
Let E be chosen from O given in Definition 2.3.1 and define A as in it. Let
Then their motions of dynamics read
˙¯x = {[f′(x(t)) − kgsp(x(t)) − kgs(x(t) − v)p′(x(t))] I + dG + gs(v − x(t))p′(x(t))S} ¯x+ ¯y− ¯z,
:= fH(t) ¯x+ ¯y− ¯z. (3.20a)
˙¯y = −10x(t)¯x − ¯y, (3.20b)
˙¯z = µb ¯x− µ ¯z. (3.20c)
Instead of calculating the transverse Lyapunov exponents of the corresponding variational equation (3.17) of equation (3.14), we would prove directly that the origin of (3.20) is asymptotically, exponentially stable. As a consequence, all transverse Lyapunov exponents of (3.14) are negative.
Sufficient conditions to obtain the synchronization of coupled HR system (3.14) are stated precisely in the following.
Theorem 3.2.1. (i) Assume x(t) satisfies synchronous equation (3.16). Let λ2 be the second largest eigenvalue of coupling matrix G. Let r2 = µ2(S), the matrix measure of S with respect to 2-norm. Here S is defined in (3.19). Set α =: −1 −rk2 and
where hα is a constant and d1 = max
x∈R f′(x) ≈ 2.253. Let d0 = sup
x(t) 10|x(t)| ≤ 20. Then coupled HR system (3.14) is locally synchronized provided that
(−λ2)d + (−hα)kgs > 24 > d0+ b, for kgs > 0, and
−λ2d > 26.253 > d0+ b + d1, for kgs = 0. (3.21c) (ii) Assume that lim
t→∞x(t) = xc. Let lim
t→∞H(t) := ff Hc. Here fH(t) is defined in (3.20a).
Then system (3.14) is locally synchronized if all real parts of eigenvalues of
Hfc I −I
−10xcI −I 0
µbI 0 −µI
=: Hc
are negative.
Proof. To obtain local synchronization of (3.14), we study equations (3.20). Note that for excitatory HR neurons, x(t) < v = 2 for all t. Clearly, µ2( fH(t)) ≤ λ2d+hαkgs =: λ.
Here hα is defined in (3.21b). Then by equation (3.20), we have
Dlk ¯x(t)k ≤ λk ¯x(t)k + k ¯y(t)k + k ¯z(t)k, (3.22a) Dlk ¯y(t)k ≤ d0k ¯x(t)k − k ¯y(t)k, (3.22b) Dlk ¯z(t)k ≤ µbk ¯x(t)k − µk ¯z(t)k. (3.22c) Applying Proposition 2.1.1-(ii), we see that the first part of the assertion of the theorem holds true provided the real parts of the eigenvalues of
B=
λ 1 1
d0 −1 0
µb 0 −µ
(3.23)
are negative. Indeed, the Routh-Hurwitz Criterion asserts that it occurs whenever
−λ > d0+ b. So the first assertion of the Theorem holds true.
The last assertion of the Theorem is a direct consequence of Proposition 2.1.2.
If the steady-state synchronization is considered, then some easier verifiable con-ditions than those stated in Theorem 3.2.1-(ii) can be obtained.
Corollary 3.2.1. Let d = 0. Assume lim
t→∞x(t) = xc. Then system (3.14) without electrical coupling is locally synchronized if the real parts of the eigenvalues of A are negative, where
Proof. Note that Hc with d = 0 has the following form.
Applying Proposition 2.1.3, we have that system (3.14) is locally synchronized provided
that have all its eigenvalues with negative real parts.
Define matrix A(x) as
Then it can be proved by applying the Routh-Hurwitz Criterion that for any y < x ≤ 0, if all eigenvalues of A(x) have positive real parts, then so do those of A(y).
Upon using the above observation and the fact that ¯λi, the real parts of eigenval-ues of S, are negative, we conclude that the assertion of the Corollary holds true.
Corollary 3.2.2. Let C be a node-balancing matrix, i.e., its row sums and column sums are equal. Assume lim
t→∞x(t) = xc. Then system (3.14) is locally synchronized if all real parts of the eigenvalues of A, as given in (3.24), are negative.
Proof. As in the proof of Corollary 3.2.1, it suffices to show that all real parts of the eigenvalues of dG + gs(v − xc)p′(xc)S are negative. However, by Theorem 2.1.1, we have λmax(dG + gs(v − xc)p′(xc)S) ≤ µ2(dG + gs(v − xc)p′(xc)S) ≤ dµ2(G) + gs(v − xc)p′(xc)µ2(S) ≤ 0. Thus, the proof of the Corollary is completed.
Remark 3.2.1. (i) To acquire synchronization of coupled networks, the second largest eigenvalue of the coupling matrix plays an inescapable and decisive role. Indeed, in certain cases, such as the system is fully coupled, the necessary and sufficient condition [83] with k = 0 for local synchronization is
hmax+ dλ2 < 0.
Here hmax is the largest Lyapunov exponent of the individual oscillator. In most of interesting networks, λ2 becomes closer to the origin from the left as the number of oscillators grows. Hence, it takes greater coupling strengths to synchronize the larger system. In other cases, such as the coupled map lattices (1.2), the system exhibits the size instability phenomena, that is, the system with the number of nodes greater than a certain critical value loses its synchrony regardless how strong the coupling strength is. Such size instability is induced by the competition between a certain eigenvalues, including λ2, of the coupling matrix, we will study it in the next chapter.
(ii) If the connection C is symmetric, then r2 is the second largest eigenvalue of S. It is easy to see that if the connection network is all-to-all coupled, then k = m − 1, r2 = −m, and so α = m−11 ≤ 1. It can be computed that the denser the network is coupled, the larger α 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.
We also mention that the computation cost to verify the synchronous conditions (3.21c) or (3.24) is very little as compared to that of computing second Lyapunov exponent of the network. Specifically, if HR system (3.14) is both electrically and chemically coupled, one needs to check the inequality (3.21c) to see if the system is synchronized. To check the steady-state synchronization, one only needs to verify the sign of the largest real part of eigenvalues of a 3 × 3 matrix, A, see (3.24), regardless of the number of neurons.