3.2 Synchronization in Hindmarsh-Rose Neurons with Chemical and Elec-
3.2.3 Further Discussions and Numerical results
In the subsection, we shall focus on applying Theorem 3.2.1 and Corollaries 3.2.1, 3.2.2 to coupled HR neurons (3.14) to extract more detailed synchronization phenomena.
To this end, we need to know, firstly, the dynamics on synchronous manifold.
(I) Dynamics on synchronous manifold
We begin with the study of the dynamics of synchronous equation (3.16). Its dynamics is to be provided numerically. 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 behavior, bursting. (see e.g., [44]). 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 (3.16) is summarized in Table 3.1. A similar result to Table 3.1 was also reported in Fig. 2 of [16]. On the synchronous manifold, the solution trajectory x(t) of (3.16), depending on initial conditions and kgs, may settle into various stable states. Figure 3.4 provides the maximum Lyapunov exponent (MLE) of synchronous equation (3.16) versus kgs. For 0 ≤ kgs ≤ 0.85, there is a set of initial conditions with positive measure for which their corresponding MLE is positive. However, for 0.809 ≤ kgs ≤ 0.85, 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 conditions with positive measure so that the solution trajectories of (3.16) converge to a stable periodic solution (see Fig. 3.5) and stable regular bursting (see Fig. 3.6), respectively. Specifically, let (xc, yc, zc) be the steady state of (3.16) (see, Fig. 3.7), and let
Cr = {(x, y, z) : |x − xc| < r, |y − yc| < r, and |z − zc| < r}, (3.25a) and
Ir = {(x, y, z) : |x − xc| > r, |y − yc| > r, and |z − zc| > r}. (3.25b) 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 from C0.02 (resp.,
I1), then its trajectory converges to a periodic orbit (resp., a stable regular bursting) (see Figs. 3.5, 3.6). Similarly, for kgs ∈ [0.814, 0.85], synchronous equation (3.16) 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. In fact, a direct calculation shows that a Hopf bifurcation occurs near 0.813.
Furthermore, if one performs the linearized stability at the steady-state (xc, yc, zc), then one sees that (xc, yc, zc) is stable whenever kgs≥ 0.814 (see Fig. 3.8). Such analysis of linearized stability provided some supportive evidence for the validity of Table 3.1.
Table 3.1: The dynamics of synchronous equation (3.16) with various range of kgs. The multi-stability of (3.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
In summary, the numerical results suggest that on the synchronous manifold, for kgs small, the regular bursting behavior of single HR persists. For kgs in an in-termediate range, the multistability of equation (3.16) occurs. Depending on initial conditions, the coexistence of multi stability states including a stable regular burst-ing and a stable periodic solution/a stable fixed point could be observed. When kgs becomes large, equation (3.16) has a globally asymptotically stable fixed point. Such complex dynamical behavior of synchronous equation (3.16) leads to the possibility of stable multi-state synchronization of coupled HR neurons (3.14). If the initial condi-tions and the range of kgs are so chosen that the corresponding synchronous equation leads to a regular bursting solution, then the associated coupled HR neurons (3.14) achieves stable regular bursting synchronization. Likewise, we define stable periodic synchronization and stable steady-state synchronization accordingly. As we can see, via Table 3.1, that for 0.809 ≤ kgs ≤ 0.85, the coexistence of stable multi-state syn-chronization of coupled HR neurons (3.14) could be observed. It should also be mention that the theory of weakly coupled oscillators has often been used to analyze networks of neuron coupled by small chemical synapses gs, see e.g. [60], and the many related
work cited therein. Using this theory enables one to obtain some extensive analytical insight. Furthermore, the ups and downs of synaptic strength can be controlled. For instant, N-methyl-aspartate receptors can both boost and dampen synaptic efficiency in the brain [10]. Such observations give the justification for the consideration of small chemical synapses gs.
0 0.01 0.02 0.03 0.04
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
X:0.85 Y:2e-4
kgs
MLE
Figure 3.4: The maximum Lyapunov exponent (MLE) of synchronous equation (3.16) is computed for various kgs. For 0 ≤ kgs≤ 0.85, MLE> 0 for a set of initial conditions with positive measure. For kgs ∈ [0.809, 0.85], there is also a set of initial conditions for which its corresponding MLE< 0.
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(t), y(t), z(t))
y x
z
Figure 3.5: The solution trajectory with randomly chosen initial conditions from C0.02
converges to a stable periodic orbit. Here kgs = 0.812.
(II) Neurons with only chemical synapse
In [6], a local steady-state synchronization of bursting neurons with no electrical
-3 -2.5-2 -1.5-1 -0.50 0.5 1 1.5 2
Figure 3.6: The solution trajectory with randomly chosen initial conditions from I1 converges to a stable regular bursting. Here kgs = 0.812.
0 10 20 30 40 50 tends to (2, −19, 14.4) as kgs tends to infinity.
coupling is studied without providing mathematical details. Moreover, their approach fails to see if synchronization of neurons can be achieved when the networks are in-termediately 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 chemical signals each neuron receives, and is independent of all other details of the network topology. From (3.21a) and (3.21c), it is clear that the larger the density α of the network is, the greater chance coupled HR neurons (3.14) gets synchronized. In the following, we shall prove that the system of coupled HR neurons (3.14) achieves
steady-state synchronization regardless how sparsely the network is coupled provided that kgs ≥ 0.87.
By Corollary 3.2.1, it suffices to show that the maximum real part λmax(A) of eigenvalue of matrix A defined in (3.24) is negative. From Fig. 3.8, we see that λmax(A) < 0 whenever kgs ≥ 0.814. Upon taking into consideration of the dynam-ics on synchronous manifold as provided in Table 3.1, we have the following conclusion.
Coupled HR neurons (3.14) with d = 0 achieves the steady-state synchronization for a set of initial conditions with positive measure whenever kgs ≥ 0.814 regardless how sparsely the network is coupled, which improves the result obtained in [6]. More-over, the system acquires the steady-state local synchronization for all initial conditions sufficiently close to each other whenever kgs≥ 0.87.
Numerically, the following scenarios are also observed. Coupled two HR neurons achieves synchronization only when kgs ≥ 0.809. Stable regular bursting and steady-state synchronization is found on a set of initial conditions with positive measure, respectively, whenever kgs ≈ 0.85. Such numerical results are illustrated in Figs. 3.9 and 3.10.
0.5 0.6 0.7 0.813 0.9
-0.5 0 1 2 3 3.5
kgs
λmax(
¯ A)
Figure 3.8: The maximum eigenvalue of the linearized operator with respect to the synchronous equation (3.16) is computed for various kgs.
0 200 400 600 800 1000 10-8
10-6 10-4 10-2 100
x1-x2
t
0 200 400 600 800 1000
-4 -2 0 2
x1
t
Figure 3.9: The time series of x1(t) − x2(t) and x1(t). The graphs demon-strate the stable regular bursting synchronization. Here gs = 0.85, initial = [−2, −18, 3, −2.5, −18.5, 2.5].
0 200 400 600 800 1000
10-8 10-6 10-4 10-2
x1-x2
t
0 200 400 600 800 1000
0 0.08 0.16
x1
t
Figure 3.10: The time series of x1(t)−x2(t) and x1(t). The graphs demonstrate the sta-ble steady-state synchronization. Here gs = 0.85, initial = [0.026, 1, 6.5, 0.126, 1.1, 6.6].
(III) Neural synchronization with only electrical synapse
For gs = 0, by Theorem 3.2.1, we obtain stable regular bursting synchronization whenever d > d0+b+d−λ2 1 =: dmin. Consider, for instance, a ring of 2K-nearest-neighbor mutually coupled networks, the predicted minimum electrical coupling strength dmin is
computed with the number m of neurons and K being given. Note that in such case
λ2 = −4 XK
i=1
sin2 iπ
m = sin (K + 12)2πm
− sin mπ
sin mπ − 2K. (3.26a)
The results are listed in Table 3.2. For instance, given a number of (104 + 1) neurons, it takes the electrical coupling strength 6.66×107or greater to reach synchrony for a network with the nearest-neighbor coupling. It only takes 2.7 × 10−3 or greater to do so for an all-to-all network. If the coupling matrix G is of high dimension and without fine structure for computers to be able to calculate its second largest eigenvalue effectively, one may use some known estimates to find the upper bound of λ2. For instance, we have that (see e.g., [70])
λ2 ≤ −2
(m − 1)¯ρ(G) − m−22 ,
(3.26b)
where ¯ρ(G) is the mean distance of the graph associated with G. Upper bounds for dmin by using (3.26b) are listed in the Table 3.2, too. For m = 104 + 1, the upper bounds of dmin are 3.29 × 108 and 65640, respectively, for a network with the nearest-neighbor coupling and all-to-all coupling. In a nutshell, connecting each neuron to more neighbors is an effective way for large-size networks to lower the synchronization threshold.
The upper bound for λ2 in (3.26b) is quite good for the sparsely coupled networks.
Indeed, in the case of the nearest-neighbor coupling, the exact value of λ2 and its estimated upper are both O(m12). On the other hand, if the network is densely coupled, the upper bound in (3.26b) gives a poor estimate for λ2. Nevertheless, if one picks other type of upper bound for λ2, better estimates could be expected. For example, it is also known, see e.g. [103], that
λ2 ≤ −m
β, (3.26c)
where β > 0 is such that βG + L is negative semidefinite. Here L is the Laplacian matrix of the complete graph, i.e., L = m(I − eeT), where e is given as in (2.7a).
For the all-to-all coupling, it is readily verified that βG + L with β = 1 is negative
Table 3.2: The first component of the pair in the table gives the predicted minimum electrical coupling strength dmin by using the exact value of λ2. For instance, with m = 102+ 1, K =m−1
4
= 25, the predicted minimum electrical coupling strength is dmin = 1.4. The second component of the pair in the table is the upper bound of dmin
by using (3.26b).
m 21 102+ 1 104+ 1
K = 1 (the nearest-neighbor coupling) (296, 1320) (6786, 32824) (6.66 × 107,3.29 × 108) K = [m−14 ] (6.10, 269) (1.40, 1320) (0.0145, 1.32 × 105) K = m−12 (the all-to-all coupling) (1.26, 138) (0.27, 663) (0.0027, 65640)
semidefinite. Consequently, the equality in (3.26c) can be achieved, which yields the best possible estimate.
(IV) Neural synchronization with both excitatory chemical and electrical synapses
Herein, networks with both excitatory electrical and chemical connections are considered. To extract more information on synchronization of the system, we further assume C to be a node-balancing matrix. We first observe that x(t) < v and p′(x(t)) ≥ 0 for all t. So hα, defined in (3.21b), is decreasing in α. If
(−λ2)d + (−h−1)kgs> 24, (3.27) then (3.21c) is also satisfied. The synchronization region satisfying (3.27) and kgs ≥ 0.87 is demonstrated in Fig. 3.11. That is to say, if (−λ2d, kgs) is chosen from the shaded region in Fig. 3.11, then multi-state or single-state synchronization can be realized depending on the range of kgs. Consider, for instance, coupled two HR neurons.
Let kgs = 0.812 and d = 30. If (xi(0), yi(0), zi(0)) ∈ C0.02 (resp., I1), i = 1, 2, as given in (3.25a) (resp., (3.25b)), and are distinct, then the stable periodic (resp., stable regular bursting) synchronization occurs (see Figs. 3.12, 3.13).
We further observe that there exists a t1 such that h−1 < 0 (resp., h−1 > 0) whenever kgs ∈ [0, t1) =: J1 (resp., kgs ∈ (t1, 0.87] =: J2) (see Fig. 3.11). Here t1 ≈ 0.6044. Hence, both chemical and electrical synapses enforce the synchronization
phenomena whenever kgs ∈ J1. For kgs ∈ J2, the chemical synapses play dragging roles for system to reach synchrony. To synchronize, the electrical synapses have to be strong enough to suppress the dragging force created by chemical synapses. Such t1 is called a turning point for h−1.
We are then led to compute turning points for hα (see Fig. 3.14). 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 synchrony of the system.
To summarize, a synchronization region is obtained in Fig. 3.11. Particularly, multi-state synchronization of coupled HR neurons can be realized whenever kgs ∈ [0.809, 0.85] and (−λ2d, kgs) lies in the synchronization region. 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.
To conclude this section, we will elaborate more on some crucial points.
(i) As evidence in Table 3.1, the multi-stable state, which depends on the choice of the initial conditions exist in excitatory HR neuron. In fact, other choice of parameters, such as a = 3, µ = 0.006, q = 3, x0 = 1.56, b = 4, would result the neurons burst irregularly (chaotically). Under such circumstance, the presence of both stable chaotic attractor and stable periodic state can be detected. And our approach can be applied to the above described scenario as well.
(ii) From inequality (3.21c), we see that the denser the coupling network is, or equivalently, the larger the density α is, the easier the system gets synchronized.
(iii) We mention that free packages SLEPc developed by V. Hern´andez, J. E. Rom´an, A. Tom´as, and V. Vidal [38] can be used to compute λ2 efficiently.
(iv) In vivo experiments, the strength of an excitatory mono-synaptic connection has a biologically realistic value gs = 0.66 × 10−3 (see e.g., [9]). Using such gs and our theory, we may conclude that if the number of presynaptic neurons that connect to a single cortical neuron is greater than 1319 and in between 1226 and 1288, then the system may reach steady-state synchronization and multistate synchronization,
respectively. It should be mentioned that the quantitative description of the circuit of cat area 17 is given in [9]. Depending on the cell type and its layer location, the presynaptic synapses could range from 0 to 3500. For instance, it seems that excitatory cell types ss4(L4), p5(L2/3), p6(L4) and p6(L5/6) are the most likely candidates to generate multistate synchronization.
10 20 30 40
0 0.5 0.6044 0.87 1 1.5
26.253 -λ2d Sync. Region
kgs
Figure 3.11: The shaded area is the synchronization region satisfied by (3.27) and kgs ≥ 0.87.
0 200 400 600 800 1000
10-10 10-8 10-6 10-4 10-2
x1-x2
t
0 200 400 600 800 1000
0 0.02 0.04 0.06
x1
t
Figure 3.12: The time series of x1(t) − x2(t) and x1(t). The graphs demonstrate the stable periodic synchronization. Here d = 30, gs = 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 200 400 600 800 1000 10-8
10-6 10-4 10-2 100
x1-x2
t
0 200 400 600 800 1000
-4 -2 0 2
x1
t
Figure 3.13: The time series of x1(t) − x2(t) and x1(t). The graphs demonstrate the stable regular bursting synchronization. Here d = 30, gs = 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.
−1 −0.8 −0.67−0.6 −0.4 −0.2 0
0.5 0.604 0.87 1
α kgs
Figure 3.14: Turning points of hα.
Chapter 4
Synchronization in Model II
In this chapter, we consider the synchronization in CMLs (1.2). Throughout the chap-ter, we, as usual, assume coupling matrix G satisfies
(i) λ = 0 is a simple eigenvalue of G and
e= √1m(1, 1, . . . , 1)T1×m is its corresponding eigenvector. (4.1)
4.1 Local Synchronization Criteria
To study the stability of the synchronous manifold M = {xi = s, ∀i} of CMLs (1.2), we consider the variational equation of (1.2):
ξ(k + 1) = DF (s(k))ξ(k) + d(G ⊗ I)DF (s(k))ξ(k)
= [I ⊗ Df(s(k)) + d (G ⊗ I) (I ⊗ Df (s(k)))] ξ(k), (4.2a) where ξ = (ξ1, . . . , ξm) and each ξi is the perturbation to the ith oscillator. Let J = P−1GP, where J = [0] ⊕ J1⊕ · · · ⊕ Jp is the real Jordan Canonical form of G.
Applying the change of variables η = (P−1⊗ I)ξ, we get η(k + 1) = [(I + dJ ) ⊗ Df(s(k))] η(k), or, equivalently, in block diagonal form,
ηi(k + 1) = [(I + dJi)k⊗ Dfk(s(1))] ηi(1)
=: Ai(k)ηi(1). (4.2b)
Let σ(A) denotes the spectrum of A. Then σ(Ai(k)A∗i(k)) equals to σ([(I + dJi)k⊗ Dfk(s(1))] [(I + dJi∗)k⊗ (Dfk(s(1)))∗])
= σ([(I + dJi)k(I + dJi∗)k] ⊗ [Dfk(s(1)) · (Dfk(s(1)))∗])
= σ((I + dJi)k(I + dJi∗)k) · σ(Dfk(s(1)) · (Dfk(s(1)))∗)
= σ((I + d ¯Ji)k(I + d ¯Ji∗)k) · σ(Dfk(s(1)) · (Dfk(s(1)))∗),
where ¯J = [0] ⊕ ¯J1⊕ · · · ⊕ ¯Jp is the Jordan Canonical form of G. Consequently, the Lyapunov exponents of (1.2) are
hj + lim
k→∞
lnp λw,i
k .
Here hj are the Lyapunov exponents of the individual system f and λw,i are the eigenvalues of (I + dJλi)k(I + dJλ∗i)k where Jλi is a Jordan block of matrix G and λi is an eigenvalue of G. Let the size of matrix Jλi be ki× ki, and let N = Jλi− λiI. It should be noted that for sufficiently large k,
(I + dJλi)k = ((1 + dλi)I + dN )k= (1 + dλi)k(I + αN )k
= (1 + dλi)k(I +
kXi−1 j=1
k j
αjNj)
=: (1 + dλi)kTi,
where α = d/(1 + dλi). Clearly, the order of the magnitude of each entry of TiTi∗ is at most O(k2ki−2). We conclude, via the Gershg¨orn disk theorem, that all eigenvalues of TiTi∗ are of the order O(k2ki−2). Consequently, the Lyapunov exponents of (1.2) are
hj + ln |1 + dλi|. (4.3)
We summarize the above as follows.
Theorem 4.1.1. Let coupling matrix G satisfy (4.1). Then CMLs (1.2) achieves local synchronization provided that
hmax+ ln |1 + dλi| < 0, i = 2, . . . , m, (4.4) where hmax is the largest Lyapunov exponent of the individual map f and λi ∈ σ(G) − {0}, i = 2, . . . , m. Otherwise, CMLs (1.2) will not acquire local synchronization.
Remark 4.1.1. The decoupling form (4.2b) of variational equation (4.2a) was first observed and proposed by Pecora and Carroll [75].
We shall assume from here on that the real parts of the eigenvalues of G are nonpositive. To find the range of the coupling d so that (4.4) is fulfilled, we need to solve the following min max problem.
mind∈R max
Here m1 is the number of eigenvalues lying in upper complex plane or on the real axis.
The curves ri(d) are termed the ith-mode of the transverse Lyapunov exponent curves.
The equalities above are due to the facts that |1 + dλi| = |1 + d¯λi|, the real parts of the eigenvalues of G are nonpositive and (4.4) is violated whenever d ≤ 0. Without lose of generality, we may assume those distinct nonzero eigenvalues are λi, i = 2, . . . , m1, with 0 < |λ2| ≤ · · · ≤ |λm1|. The coupling value d := dc solving the min max problem (4.5) is the optimal choice of the coupling in the sense that it gives the fastest convergent rate of the initial values toward the synchronous state. To understand how r(d) is formed, we need to know the ordering of ri(d). For d > 0, direct computations yield
ri(d) =
Let Ai = {j : 2 ≤ j ≤ m1 and |λi| = |λj|}. Then max
j∈Ai |1+dλj| = |1+dλw|, where w is so chosen that Re(λw) ≥ Re(λj) for all j ∈ Ai. This gives that within each of the index set Ai, their corresponding quantities |1+dλi| are well ordered for any d > 0. Consequently, to solve (4.5), we may assume, without loss of generality, that 0 < |λ2| < · · · < |λm1| from here on. Using the terminology in [37], we see that the numbers 2 and m1
correspond to the longest and shortest wavelength modes, respectively. The numbers in between 2 and m1 are to be called intermediate wavelength modes. Since dij = dji, for any i and j, we consider only dij with i > j. Our reduction process is now complete.
The following procedures are proposed to determine the “actual” node points of r(d) from the candidate set {dij : i > j}.
(A) Set k0 = 0, and k1 = max{l| Re(λi) ≤ Re(λl), ∀i = 2 . . . , m1}. Let k2 be the largest index so that 0 < dk2k1 ≤ dwk1 for all k1 < w ≤ m1.
(B) Let k3 be the largest index so that dk2k1 < dk3k2 ≤ dwk2 for all k2 < w. The process can be continued until kp = m1 for some p ≤ m1.
The next result shows that {ki}pi=1 is the set of “actual” node points of r(d).
Theorem 4.1.2. Let coupling matrix G satisfy (4.1). Assume that the real parts of the eigenvalues of G are nonpositive. Then r(d) = rki(d) whenever dkiki−1 ≤ d ≤ dki+1ki, i = 1, . . . , p. Here dk1,k0 = 0 and dkp+1kp = ∞.
Proof. Denote by Ij = [dkjkj−1, dkj+1kj]. It then follows from (4.6c) that if i > j and dij > 0, then ri(d) > rj(d) whenever d > dij and ri(d) < rj(d) whenever 0 < d < dij. We then conclude that
(i) the ordering of ri(d) and rj(d) remains the same until both curves meet; (4.7a) (ii) if ri(d∗) > rj(d∗) for some d∗ > 0 with i > j, then ri(d) > rj(d) for all d ≥ d∗.
(4.7b) Using the first inequality in (4.6a), we have that r(d) = rk1(d) for ǫ1 > d ≥ 0. Here ǫ1 is sufficiently small. It then follows from (4.7a), (4.7b) and procedure (A) that r(d) = rk1(d) on I1. Upon using (4.7a), we conclude that r(d) = rk2(d) for d ∈ (dk2k1, dk2k1+ǫ2).
Here ǫ2 is sufficiently small. Similarly, r(d) = rk2(d) on I2. We omit the proof of the remaining assertions of the theorem due to the similarity.
Note that not all cki given in (4.6a) could be critical points of r(d). In fact, the critical points of r(d) may not even come from the set {cki}. We next identify the
“actual” critical points of r(d). Our next main result shows that r(d) has exactly one critical point.
Theorem 4.1.3. The curve r(d) has a unique critical point dc that solves the min max problem (4.5). Moreover, the optimal range of coupling d to sustain stably synchronous chaos of (1.2b) is (dl, dr). Here dl and dr, dl< dr, are the intersection points (if exist) of the straight line y = e−hmax and the curve y = r(d). Consequently, CMLs (1.2b) acquires local synchronization if and only if d ∈ (dl, dr).
Proof. We break up the proof of the theorem into the following three steps.
(Step I) We first claim that the number of cki lying in the interiorI◦i of Iiis at most one. Indeed, suppose there exist cka ∈I◦aand ckb ∈I◦bwith cka < ckb. Then the following hold true. (i) rka(ckb) > rka(cka). (ii) rka(cka) > rkb(cka). (iii) rkb(cka) > rkb(ckb).
Inequalities (i) and (iii) hold true since cka and ckb are, respectively, the minimum points of rka(d) and rkb(d). The fact that rka(d) lies above all other curves on Ia leads to inequality (ii). Combining these inequalities, we have that rka(ckb) > rkb(ckb), a contradiction to the fact that rkb is the maximum curve on Ikb.
(Step II) We next show that if cki ∈ I◦i, then r(d) is decreasing on (0, cki) and increasing on (cki, ∞). Indeed, for d ∈ Ii+1, r(d) = rki+1(d) > rki(d) > rki(dki+1ki) = rki+1(dki+1ki). Using the conclusion in Step I and the fact that rk2i(d) is parabolic, we conclude that rki+1(d) must be increasing on Ii+1. On the other hand, rki−1(d) must be decreasing on Ii−1 since rki−1(d) > rk(d) > rki(dkiki−1) = rki−1(dkiki−1). The monotonicity of r(d) on each interval Ij, 1 ≤ j ≤ m1, can be similarly determined.
(Step III) Since r(d) is decreasing initially on I1 and increasing eventually on Ip, there must have at least one critical point. If such points do not lie in the set of node points, then r(d) has a unique critical point. Suppose cki ∈/ I◦i for all i = 1, . . . , p. Then r(d) is monotonic on each interval Ii. Suppose r(d) first changes its monotonicity at
dkl+1kl for some l. Then a similar argument as given in the Step II shows that once r(d) becomes increasing on Il+1 it will stay increasing the rest of the way. We have just completed the proof of the theorem.
Now, we define some terminologies occurring in Theorem 4.1.3. The curve r(d) given in is to be called the synchronization curve of CMLs (1.2), and the interval Nm,f := (dl, dr) is termed the synchronization interval. (1.2), if exists. The value dc is called the optimal coupling strength of CMLs. The value r(dc) is called the synchronization index of CMLs since local synchronization occurs iff it is less than e−hmax. The value hmax + ln r(dc) is called the Lyapunov index of CMLs.
Remark 4.1.2. (i) The optimal coupling strength dc ∈ (dl, dr) and depends only on the connectivity topology. (ii) If the straight line y = e−hmax and the curve y = r(d) do not intersect, then CMLs (1.2) will not achieve synchronization for any coupling strength.
Suppose dr and dl exist. Then as soon as d exceeds dr, a certain wavelength mode is excited, which, in turn, causes the instability of the synchronous state. The illustration in Examples 2 and 3 shows that the excited wavelength mode could be either the shortest wavelength mode, the intermediate wavelength mode or the longest wavelength mode.
In any event, dr is the exact critical value where wavelength bifurcation (WB), as the terminology using in [37], occurs. On the other hand, dlis the exact critical value where all wavelength modes become deexcited enough to induce the stability of the synchronous state.
Theorem 4.1.4. Suppose the coupling matrix G ∈ Rm×m has nonpositive real eigen-values. Denote by {λi}mi=21 the distinct nonzero eigenvalues of G. Then
(i) The synchronization curve is r(d) =
λ2(d), d ∈ [0, dm12] = I1
λm1(d), d ∈ (dm12, ∞) = I2 . (4.8a) (ii) The optimal coupling strength is
dc = dm12 = −2 λ2+ λm1
. (4.8b)
(iii) The synchronization interval is, if exists, Nm,f :=
1 − e−hmax
−λ2 ,1 + e−hmax
−λm1
. (4.8c)
(iv) The synchronization index is
t2,m :=
λ2− λm1
λ2+ λm1
. (4.8d)
(v) The Lyapunov index is
δm,f := hmax+ ln |t2,m|. (4.8e) Consequently, depending on the quantity of hmax, either CMLs (1.2) achieves no syn-chronization or short wavelength bifurcation occurs as d varies.
δm,f := hmax+ ln |t2,m|. (4.8e) Consequently, depending on the quantity of hmax, either CMLs (1.2) achieves no syn-chronization or short wavelength bifurcation occurs as d varies.