^{†}Corresponding author. Email: guhuaguang@tongji.edu.cn
^{*}Project supported by the National Natural Science Foundation of China (Grant Nos. 11372224 and 11402039) and the Fundamental Research Funds for Central Universities designated to Tongji University (Grant No. 1330219127).
We investigated the synchronization dynamics of a coupled neuronal system composed of two identical Chay model neurons. The Chay model showed coexisting period1 and period2 bursting patterns as a parameter and initial values are varied. We simulated multiple periodic and chaotic bursting patterns with non(NS), burst phase (BS), spike phase (SS), complete (CS), and lag synchronization states. When the coexisting behavior is near period2 bursting, the transitions of synchronization states of the coupled system follows very complex transitions that begins with transitions between BS and SS, moves to transitions between CS and SS, and to CS. Most initial values lead to the CS state of period2 bursting while only a few lead to the CS state of period1 bursting. When the coexisting behavior is near period1 bursting, the transitions begin with NS, move to transitions between SS and BS, to transitions between SS and CS, and then to CS. Most initial values lead to the CS state of period1 bursting but a few lead to the CS state of period2 bursting. The BS was identified as chaos synchronization. The patterns for NS and transitions between BS and SS are insensitive to initial values. The patterns for transitions between CS and SS and the CS state are sensitive to them. The number of spikes per burst of nonCS bursting increases with increasing coupling strength. These results not only reveal the initial value and parameterdependent synchronization transitions of coupled systems with coexisting behaviors, but also facilitate interpretation of various bursting patterns and synchronization transitions generated in the nervous system with weak coupling strength.
Synchronization is widespread in nature.^{[1, 2]} Chaos synchronization has been one of the most important issues.^{[1, 3]} In particular, synchronization was experimentally detected in the nervous system of the mammalian visual cortex, the human thalamocortical area, and the stomatogastric ganglion of the spiny lobster.^{[4– 6]} Transitions between various synchronization states were observed in the oddball task experiment and evidenced in the human cardiorespiratory system.^{[7– 9]} It is believed that synchronization and synchronization transition are closely related to nervous information transmission and processing.^{[4, 8– 11]} Physiologically, bursting is an important neural rhythm in the nervous system and is also important for information processing.^{[12– 16]} Recently, there have been various simulations of synchronization and synchronization transition in the nervous system.^{[17– 22]}
In most investigations of synchronization, the unit in the network or coupled system always manifests monostable behaviors. With the exception of monostable behavior, multistable and/or coexisting states are also common behaviors of a nonlinear system. The coexisting behavior of the dynamic system for a given set of parameters is determined by the initial values. Multistable behaviors occur in multiple fields, including biological science.^{[23]} In the nervous system, the coexistence of resting condition and firing, ^{[24– 26]} of spiking and bursting, ^{[27– 29]} and of two different bursting patterns^{[30]} were observed in experiments and simulated in theoretical models.
Recently, a few works have focused on the synchronization of unidirectionally coupled systems, ^{[31– 38]} which manifest, when isolated, multistable or coexisting behaviors. These studies are briefly described as follows. Two coupled Rö sslerlike electronic circuits with coexisting chaotic attractors were investigated experimentally and numerically.^{[31– 33]} In that work the synchronization transitions moves from asynchronous motion to intermittent anticipating phase synchronization and generalized synchronization moves to a completely synchronized state as the coupling strength increases. The behaviors of asynchronous synchronization jump between the coexisting attractors. The synchronization dynamics of two coupled bistable Hé non maps were investigated.^{[34]} When two periodic orbits coexist, the slave map evolves into the master map state. When chaotic and periodic attractors coexist, very complex dynamics, such as the emergence of new attractors are observed. Synchronization of coupled optical bistable systems with coexisting different attractors (fixed point and chaos, fixed point and one periodic orbit, and two periodic orbits) was investigated.^{[35]}
Depending on both the control parameters of an isolated system and the coupling strength, different bifurcations and diverse dynamic regimes including periodicity and chaos occur in the route from asynchronous motion to complete synchronization. There also exists work on the nervous systems with bidirectional coupling.^{[36, 37]} One previous study^{[36]} investigates the synchronization transitions of two coupled neurons with coexisting spiking and bursting. Three different configurations of a coupled system, one neuron spiking and the other bursting, both types of spikings and both types of burstings, are considered. The complete synchronization at stronger coupling strength and synchronization transitions depends on the configurations. For the configuration of one spiking and the other bursting, the synchronization transitions for three parameter values within the coexisting parameter region was investigated.^{[37]} The results showed that the complete synchronization and synchronization transition depend on the parameter values. In both investigations, only one initial value is considered, and the influences of different initial values on the complete synchronization and synchronization transitions have not been addressed.
In this paper, we comprehensively study the synchronization and synchronization transitions of two bidirectionally coupled Chay model neurons with coexisting burstings. We considered different configurations of the network, the initial values and the control parameters. We investigated different synchronization states, various bursting patterns corresponding to different synchronization states, and complex synchronization transitions, which range from nonsynchronization to complete synchronization and is composed of different synchronization states. We identified the dependence of the complete synchronization on the initial values and control parameters. We also identified the synchronization transitions containing transitions between spike and burst phase synchronization (i.e., chaos synchronization), and transitions between spike and complete synchronization. The rest of this paper is organized as follows. Section 2 describes the Chay model and the coexisting burstings. Section 3 introduces the network model, initial values, and criteria for synchronization. The synchronization states and the corresponding rhythms are introduced in Section 4. Section 5 provides the detailed synchronization transitions and the influence of the parameter and initial values. Section 6 is the conclusion and discussion.
The Chay model was initially set up to describe the bursting patterns of β cells^{[39]} and manifested bifurcations similar to those observed in an experimental neural pacemaker.^{[40– 42]} The Chay model is described using the following three differential equations:
Here, V, n, and C are the membrane potential, the probability of opening the voltagesensitive K^{+ } channel, and the intracellular concentration of Ca^{2+ }, respectively. V_{i}, V_{k}, and V_{l} are the reversal potentials for mixed Na^{+ }– Ca^{2+ }, K^{+ }, and leakage ions, respectively. g_{i}, g_{kv}, g_{kc}, and g_{l} represent the maximum conductance of the mixed Na^{+ }– Ca^{2+ }, voltagedependent K^{+ }, calciumdependent potassium, and leakage ions, respectively. m_{∞ } and h_{∞ } are the probabilities of activation and inactivation of the mixed inward current channel, and n_{∞ } is the steady state value of n. m_{∞ }, h_{∞ }, and n_{∞ } could be expressed by y_{∞ } = α _{y}/(α _{y}+ β _{y}) in general, y stands for m, h, or n, and α _{m} = 0.1(25+ V)/(1− e^{− 0.1V− 2.5}), β _{m} = 4e^{− (V+ 50)/18}, α _{h} = 0.07e^{− 0.05V− 2.5}, β _{h} = 1/(1+ e^{− 0.1V− 2)}, α _{n} = 0.01(20+ V)/(1− e^{− 0.1V− 2}), β _{n} = 0.125e^{− (V+ 30)/80}. The relaxation time of the voltagegated K^{+ } channel is τ _{n} = 1/(λ _{n}(α _{n}+ β _{n})). K_{C} is the rate constant for the efflux of the intracellular Ca^{2+ } ions, and ρ is a proportionality constant. The parameter values of the Chay model are as follows: g_{i} = 1800, g_{kv} = 1700, g_{kc} = 10, g_{l} = 7, V_{k} = − 75, V_{l} = − 40, V_{i} = 100, K_{c} = 3.3/18, λ _{n} = 227.5, ρ = 0.27. The unit of potential, conductance, and time is mV, pS, and sec, respectively. Detailed descriptions of the Chay model are given in a previous study.^{[39]}
The parameter V_{C}, the reversal potential of Ca^{2+ }, here serves as the bifurcation parameter. The Chay model is solved using the fourthorder Runge– Kutta integration method with time step of 0.001 sec.
When V_{C} = 308.1, the coexistence of period1 bursting and period2 bursting can be simulated with the Chay model. When the initial values of (V, n, C) are set to (− 47.733563, 0.105692, 0.511172) for period1 bursting and (− 48.365978, 0.100857, 0.518610) for period2 bursting, the Chay model exhibits period1 bursting (Fig. 1(a)) and period2 bursting (Fig. 1(b)), respectively. The interspike interval (ISI) of period1 bursting is a constant value, and the period2 bursting exhibits two different ISIs.
Bifurcation of ISI with respect to V_{C} is shown in Fig. 2. As V_{C} increases, period2 bursting transforms to period1 bursting at V_{C} = 309.78. As V_{C} decreases, period1 bursting transforms to period2 bursting at V_{C} = 307.95. V_{C}∈ [307.95, 309.78] is the region where period1 and period2 bursting coexist. When V_{C} is within this range, the Chay model exhibits either period1 bursting or period2 bursting; the choice depends on the initial values of (V, n, C). When V_{C} > 309.78, the behavior is monostable period1 bursting. When V_{C} < 307.95, the behavior is monostable period2 bursting.
Because the attraction domain in the threedimensional space (V, n, C) is difficult to visualize, we use the attraction domain in a twodimensional crosssection (V, C) at a fixed n value which is easier to view. As shown in Figs. 3(a)– 3(c), when V_{C} = 308.1, the attraction domain of period1 bursting (black) and of period2 bursting (white) in the (V, C) plane, corresponds to n = 0.25, 0.15, 0.3, respectively. The attraction domains in the (V, C) plane corresponding to n = 0.25, 0.15, 0.3, when V_{C} = 309.75 are shown in Fig. 3(d)– 3(e), respectively. The ranges of V and C in Fig. 3(a) and Fig. 3(d) are much larger than in the other four figures. The ranges of V and C of Figs. 3(b)– 3(c) and Figs. 3(e)– 3(f) nearly equal those of the stable period1 and period2 bursting patterns. As can be seen from Fig. 3, as V_{C} increases, the black area increases while the white area decreases, thus showing that the attraction domain of period1 bursting increases while that of period2 bursting decreases.
The coupled system is constructed by two bidirectionally coupled identical Chay neurons with a gap junction, ^{[43]} which is an important coupling manner that encourages synchrony.^{[44]} The equations of the coupled system are as follows:
Here, the subscripts 1 and 2 indicate neurons 1 and 2, respectively. The parameter g_{c} is the coupling strength. It varies from 0 to 35 in steps of 0.01. The values of other parameters are the same as those of a single neuron. The coupling item of neuron 1 is g_{c}(V_{2}− V_{1}). The equations are solved numerically using the fourthorder Runge– Kutta integration method with integration time step of 0.001 sec.
Two values of V_{C} within the coexistence region are chosen, indicated by two arrows in Fig. 2. One is 308.1 mV, which is close to the monostable period2 bursting. The other is 309.75 mV, which is close to the monostable period1 bursting.
When V_{C} = 308.1, eight initial values with equal intervals spanning both period1 bursting and period2 bursting, as indicated by dots in Figs. 1(a) and 1(b), are numbered as a_{i} and b_{i} (i = 1, 2, … , 8), respectively, from left to right. For the configuration when one neuron has period1 bursting and the other period2 bursting, the 64 initial values a_{i}b_{j} (i = 1, 2, … , 8; j = 1, 2, … , 8) are different. For the configuration when both neurons show period1 bursting, one can ignore the case of i = j (i = 1, 2, … , 8) that lead to couplingindependent compete synchronization and the case of symmetry of i and j; then there are (64− 8)/2 = 28 different initial values a_{i}a_{j} (i = 1, 2, … , 7; j = i+ 1, … , 8). Similarly, there are 28 different b_{i}b_{j} (i = 1, 2, … , 7; j = i+ 1, … , 8) values for configuration of when both neurons show period2 bursting. Thus, a total of 120 (64 + 28 + 28) different initial values are used for V_{C} = 308.1.
Similarly, 120 different initial values are also used for V_{C} = 309.75.
We investigate complete synchronization and phase synchronization. Both spike and burst phase synchronizations are considered. The spike phase function (ϕ _{l}(t)) of the firing pattern of neuron l (l = 1, 2) can be defined as follows:
Here,
(i) The complete synchronization (CS) state when max(Δ ϕ (t)) = 0.
(ii) The spike phase synchronization (SS) state when 0 < max(Δ ϕ (t)) ≤ 2π .
(iii) Nonspike phase synchronization state when max(Δ ϕ (t)) > 2π .
When the behavior of the two coupled neurons exhibits nonspike phase synchronization, the burst spike synchronization is considered. The burst phase function (Φ _{l}(t)) of firing pattern of neuron l (l = 1, 2) can be defined as follows:
Here,
(1) Burst phase synchronization (BS) state when 0 < max(Δ Φ (t)) ≤ 2π .
(2) Nonburst phase synchronization (NS) state when max(Δ Φ (t)) > 2π .
In addition, lag synchronization (LS) is also considered for some firing patterns. Spike trains of neuron 1 and neuron 2 overlap for LS of firing patterns if one neuron moves forward or backward along time axis and the other remains fixed.
CS of period1 bursting and period2 bursting can be simulated, as shown in Fig. 4. For example, spike trains of, period1 bursting of CS are shown in the upper two rows of Fig. 4(a) for the conditions when V_{C} = 309.75, initial value b_{5}b_{6}, and g_{c} = 5.0. The spike trains of period2 bursting of CS are shown in the upper two rows of Fig. 4(b) for the conditions when V_{C} = 308.1, a_{1}b_{2}, and g_{c} = 4.0. The spike phase differences of both bursting patterns equal zero, as shown in the bottom row of Figs. 4(a) and 4(b).
Spike phase synchronization (SS) can be simulated and the corresponding firing patterns include various periodic and chaotic bursting patterns.
A period3 bursting (V_{C} = 308.1, initial value a_{1}b_{2}, and g_{c} = 1.08) is shown in Fig. 5(a), containing one burst per period and three spikes per burst. A period(3, 3) bursting (V_{C} = 308.1, initial value a_{1}b_{2}, and g_{c} = 1.3) is illustrated in Fig. 5(b), composed of two bursts per period and three spikes in the both bursts. The ISIs of the three spikes in the first burst are different from those of the second. A period(7, 8) bursting (V_{C} = 308.1, initial value a_{1}b_{2}, and g_{c} = 3.3) is shown in Fig. 5(c), containing two bursts per period. One burst contains seven spikes and the other eight spikes. In addition, periodk (k = 2, 4) bursting, period(k, k) (k = 2, 4, 5, 6) bursting and period(k, k + 1) (k = 2, 3, 4, 5, 6) bursting are simulated.
For the periodk (k = 2, 3, 4) bursting and period(k, k) (k = 2, 3, 4, 5, 6) bursting, spikes in the two firing trains within the same time duration exhibit onetoone correspondence, so max(Δ ϕ (t)) is less than 2π , as shown in the lowest row of Figs. 5(a) and 5(b). For period(k, k + 1) bursting (k = 2, 3, 4, 5, 6), the burst with k spikes of one neuron corresponds to the burst with k + 1 spikes of the other neuron. The max(Δ ϕ (t)) almost approaches but remains below 2π (bottom row of Fig. 5(c)). All of the above results show that the periodic bursting events are in the SS state.
Chaotic bursting patterns include chaos2 bursting and chaos(1, 2) bursting. Chaos2 bursting is composed of bursts with two spikes. Chaos(1, 2) bursting is composed of bursts containing one spike or two spikes. An example of chaos2 bursting (V_{C} = 308.1, initial value a_{1}b_{2}, and g_{c} = 0.45), and an example of chaos(1, 2) (V_{C} = 309.75, initial value a_{5}a_{7}, and g_{c} = 0.01) are shown in Figs. 6(a) and 6(b), respectively. The largest Lyapunov exponents (LEs) of the two firing patterns are 0.58 and 0.26 greater than 0, respectively. The chaos(1, 2) bursting manifests behaviors that jump between the two coexisting attractors, which is similar to the behaviors that jump between two coexisting chaotic attractors in a previous study.
For the chaos2 bursting, the max(Δ ϕ (t)) is bounded within 2π due to onetoone correspondence between spikes from the two neurons. For the chaos(1, 2) bursting, bursts with one spike appear continuously and are interrupted by single bursts with two spikes. Bursts with two spikes do not appear continuously. The max(Δ ϕ (t)) is less than 2π .
Bursting patterns of BS include three types: another kind of chaos(1, 2) bursting, chaos(2, 3) bursting composed of bursts with 2 and 3 spikes, and chaos(2, k, k + 1) (k = 3, 4) bursting composed of bursts with 2, k, and k + 1 spikes. A chaos(1, 2) bursting (V_{C} = 309.75, initial value a_{5}a_{7}, and g_{c} = 0.11), a chaos(2, 3) bursting (V_{C} = 308.1, initial value a_{1}b_{2}, and g_{c} = 0.82), and a chaos(2, 3, 4) bursting (V_{C} = 308.1, initial value a_{1}b_{2}, and g_{c} = 1.70) are shown in Figs. 7(a)– 7(c), respectively.
The three firing patterns show nonspike synchronization, as shown by the third row of Figs. 7(a)– 7(c). As can be seen from the second burst to the fourth burst of Fig. 7(a), three spikes appear in neuron 1 and five spikes appear in neuron 2. A procedure changing from uncoupling to coupling of the behaviors of the two neurons is shown in Fig. 7(b). Neurons 1 and 2, when isolated (before t = 0 s), exhibit period1 bursting and period2 bursting, respectively. The two neurons are coupled together at t = 0 s. Then, two coupled neurons immediately exhibit nonspike phase synchronization but burst phase synchronization of bursting pattern with two or three spikes. Some bursts with two spikes of one neuron may correspond to the bursts with three spikes of the other neuron. For example, there are three spikes in the first and fourth bursts of neuron 1 while there are two spikes in the two bursts of neuron 2. The spike phase difference exceeds 2π . The evolvement of chaos(2, 3, 4) and changes of phase difference with varying time are illustrated in Fig. 7(c). Different from Fig. 7(b), the spike phase synchronization lasts about 180 s before changing to nonspike synchronization but burst synchronization. Within the 180 s, the spike phase difference less than π lasts about 150 s at first, and then changed to a value between π and 2π .
For the three chaotic firing patterns, bursts from two firing trains within the same time span exhibit a onetoone correspondence, so the max(Δ Φ (t)) is less than 2π , as shown in the lowest row of Figs. 7(a)– 7(c).
A chaos(1, 2) bursting when V_{C} = 309.75, initial value a_{3}b_{1}, and g_{c} = 0.08 is shown in Fig. 8. When time falls between 2410 and 2480 s (Fig. 8(a)), the number of spikes of neuron 1 is 18, while that of neuron 2 is 15. When time falls between 2017 and 2120 s (Fig. 8(b)), there are 17 bursts for neuron 1 and 19 bursts for neuron 2. Both Δ ϕ (t) and Δ Φ (t) are larger than 2π , showing that the bursting pattern is nonspike or nonburst phase synchronization.
Period(k, k) (k = 2, 3, 4, 5, 6) and period(k, k + 1) (k = 2, 3, 4, 5, 6, 7) bursting with spike phase synchronization are also identified as lag synchronization. For instance, the spike trains of neuron 1 and neuron 2 can overlap if one of the spike trains is moved along the time axis and the other is fixed, as shown by the dashed boxes and arrows in Figs. 5(b) and 5(c). The time difference between two spike trains is also shown in the figures.
Different synchronization states and the corresponding firing patterns are summarized in Table 1.
As g_{c} increases, the two coupled neurons achieve complete synchronization states via a complex processes composed of different synchronization states. For convenience, the transitions between synchronization states with increasing g_{c} is called a synchronization transition process in the present paper. The dependence of parameter V_{C} and the initial values on the transition processes are investigated. One hundred and twenty initial values are considered for both V_{C} = 308.1 and V_{C} = 309.75.
The transitions with an initial value a_{1}b_{2} is shown in Fig. 9. The ISIs of neurons 1 and 2 are shown in Figs. 9(a) and 9(b), respectively. According to the maximum values of spike and burst phase difference shown in Figs. 9(c) and 9(d), the transitions are divided into three parts, as labeled by vertical bold dashed lines. Part 1 exhibits transitions between SS and BS (g_{c} < 2.19). Part 2 manifests transitions between SS and CS (g_{c} ∈ [2.19, 3.32]). Part 3 is the CS state (g_{c} > 3.32).
Figure 10(a) is an enlargement of part 1 in Fig. 9(a). The SS and BS transit four times. The four SS states exist in ranges of g_{c} ∈ [0.01, 0.81], [0.83, 1.11], [1.22, 1.67], and [1.73, 2.17], respectively, and the ranges of the four BS states areg_{c} = 0.82, [1.12, 1.21], [1.68, 1.72], and 2.18, respectively. According to the largest Lyaponov exponent (LE) shown in Fig. 10(b), the BS is identified as chaos synchronization and most of the SS is period synchronization with exception of the chaos2 bursting at g_{c} ∈ [0.43, 0.49]. The transition between BS and SS is also the transition between chaos synchronization and period synchronization.
As g_{c} increases, the bursting patterns appear in the following order: period2, chaos2, period(2, 2), chaos(2, 3), period(2, 3), period3, chaos(2, 3), period(3, 3), period(3, 4), period4, chaos(2, 3, 4), period(4, 4), period(4, 5), and chaos(2, 4, 5), as shown in Fig. 10(a). The largest number of spikes within a burst of bursting patterns of spike and burst synchronization increases with increasing g_{c}.
The detailed transitions of part 1 in Fig. 9 is not very sensitive to the initial values. The 120 transition processes corresponding to different initial values exhibit almost the same transition processes as shown in Fig. 10(a). For instance, a transition process with initial value a_{7}a_{8}, exhibit almost the same process as the one with initial value a_{1}b_{2}, as shown in Fig. 10(c). The difference between the two processes is described as follows. As shown in Fig. 10(a), the period(3, 3) bursting appears when g_{c} ∈ [1.22, 1.45] and period(3, 4) appears when g_{c}∈ [1.46, 1.66]; figure 10(c) shows that the period(3, 3) bursting appears when g_{c} ∈ [1.22, 1.48] and period(3, 4) appears when g_{c} ∈ [1.49, 1.66]. The processes within the remaining ranges of g_{c} shown in Fig. 10(a) are the same as those shown in Fig. 10(c).
Part 2 of the process shown in Fig. 9(a) is highlighted in Fig. 11(a). The right boundary is at g_{c} = 3.32. The transition between CS state and SS state happens three times. The ranges of the three CS states are g_{c}∈ [2.19, 2.43], [2.61, 2.84], and [2.99, 3.23] and of the three SS states areg_{c}∈ [2.44, 2.60], [2.85, 2.98], and [3.24, 3.32]. Asg_{c} increases, ignoring the period2 bursting of CS, the following types appear in order: period(5, 5), period(5, 6), period(6, 6), period(6, 7), and period(7, 8) bursting of SS. The period(7, 8) bursting is the final bursting. This shows that the largest number of spikes within a bursting pattern of SS increases with increasing g_{c}.
Part 2 of the process with initial values a_{7}b_{8} is shown in Fig. 11(b). The right boundary is at g_{c} = 2.60. The transition between the CS state and SS state happens only once. The difference between these two transition processes with initial values a_{1}b_{2} and a_{7}a_{8} indicates that the detailed process is sensitive to initial values. The right boundary and times of transitions between SS and CS depend on the initial values.
Figure 11(c) illustrates the dependence of 120 initial values on the detailed processes. There are 68 processes that transit three times, and the right boundary of g_{c} is between 3.23 and 3.33. Twentythree processes transit two times, and g_{c} of the right boundary is between 2.83 and 3.03. Twentytwo processes transits one time, and the g_{c} values of the right boundary are between 2.38 and 2.63. The remaining seven processes achieve the CS states of period2 bursting that remains unchanged as g_{c} increases and do not experience transition.
Part 3 of the process is in the CS state and is also sensitive to the initial values. The left boundary of g_{c} for part 3 depends on the right boundary of part 2. Out of the 120 processes 109 exhibit CS of period2 bursting, which remains unchanged as g_{c} increases, as shown in Fig. 9(a); the remaining 11 processes manifest a transition from period2 bursting to period1 bursting of CS, which remains unchanged as g_{c} increases. This is shown in Fig. 12. Table 2 summarizes the details of how the bursting patterns of CS depend on the initial values. The results show that most initial values lead to period2 bursting of CS, while only a few lead to period1 bursting.
A detailed transition processes with initial value b_{5}b_{6} can be divided into four parts labeled by the bold vertical dashed lines in Fig. 13. Part 1 is almost in the nonburst synchronization (NS) state. Part 2 exhibits transitions between SS and BS states. Part 3 manifests transitions between SS and CS states. Part 4 is the CS state.
Part 1 exhibits the NS state of chaos(1, 2) bursting (Fig. 8) when g_{c} < 0.12, with three exceptions: g_{c} = 0.01, 0.03, and 0.11. When g_{c} = 0.01 and 0.03, the chaos(1, 2) bursting is in the SS state (Fig. 6(b)). When g_{c} = 0.11, the chaos(1, 2) bursting exhibit BS (Fig. 7(a)). This part is insensitive to the initial value.
Figures 14(a) and 14(b) are enlargement of part 2 from Fig. 13(a) as well as the largest Lyapunov exponent (LE) of neuron 1. The transition between SS and BS states happen three times when g_{c}∈ [0.12, 1.70]. These transitions also represent the transitions between chaos and period synchronizations. The ranges of the three SS states are g_{c}∈ [0.12, 0.82], [0.84, 1.12], and [1.22, 1.69], and those of the three BS states areg_{c} = 0.83, [1.13, 1.21], and 1.70. As g_{c} increases, period2, chaos2, period(2, 2), chaos(2, 3), period(2, 3), period3, chaos(2, 3), period(3, 3), period(3, 4), period4, and chaos(2, 3, 4) bursting appear in order.
The process with initial value a_{3}b_{1} is shown Fig. 14(c). This exhibits almost the same process as that shown in Fig. 14(a), showing that the detailed process of part 2 is insensitive to the initial values. The small difference between the processes shown in Figs. 14(a) and 14(c) is not addressed in this paper. All 120 transition processes exhibit almost the same process.
Part 3 of the process shown in Fig. 13(a) and that of neuron 1 with initial value a_{3}b_{1} are shown in Figs. 15(a) and 15(b), respectively. The g_{c} values of right boundary and the times of transition between SS and CS are different. For example, the right boundary of Figs. 15(a) and 15(b) are g_{c} = 2.96 and g_{c} = 3.30, respectively.
Figure 15(c) shows how part 3 of the 120 processes depends on initial values. The g_{c} values of the right boundary of 70 processes are between 3.26 and 3.34, and transit 4 times between SS and CS. Seventeen processes transit three times, and g_{c} of the right boundary is between g_{c} = 2.86 and 3.01. Twentyfive processes transit two times, and the g_{c} value corresponding to the right boundary is between 2.26 and 2.63. Six processes transit one time, and g_{c} of the right boundary is between 1.89 and 2.19. The remaining two processes achieve the CS of period2 bursting that remains unchanged as g_{c} increases and do not transit.
Part 4 of the process is the CS state. This part is also sensitive to the initial value. Ninetyfive out of the 120 process transit from period2 bursting to period1 bursting of CS, which remains unchanged as g_{c} increases. One example of this case is shown in Fig. 13(a). Two processes manifest the CS of period1 bursting, which remains unchanged as g_{c} increases. The remaining 23 initial values lead to CS of period2 bursting, which also remains unchanged as g_{c} increases, as shown in Fig. 16.
The number of initial values that lead to unchanged period1 bursting and period2 bursting of CS is shown in Table 3. When V_{C} = 309.75, most initial values lead to period1 bursting of CS, while only a few lead to period2 bursting, which is just the opposite of the result when V_{C} = 308.1. This is consistent with attraction domains to a certain extend. When V_{C} = 308.1, the attraction domain of period2 bursting is much bigger than that of the period1 bursting. When V_{C} = 309.75, the attraction domain of period1 bursting increases.
Three differences are extracted by comparing the transition processes of V_{C} = 308.1 and 309.75. First, a NS state exists in the transition process of V_{C} = 309.75; while no NS state exist in the transition process of V_{C} = 308.1. Second, the minimum coupling strength for CS of network with V_{C} = 309.75 (g_{c} = 1.71) is lower than that of network with V_{C} = 308.1 (g_{c} = 2.19). Third, when V_{C} = 308.1, the network is more likely to achieve the CS of period2 bursting that corresponds to large g_{c}. When V_{C} = 309.75, the network is more likely to achieve CS of period1 bursting.
The increase in the maximum number of spikes within a burst with g_{c} increasing can be interpreted by the coupling item. For example, the number of spikes within a burst of neuron 1 when V_{C} = 308.1 and initial value a_{1}b_{2} is shown in Fig. 17(a). Correspondingly, the maximum values (Fig. 17(b)) of the coupling item of neuron 1 increase and the minimum values (Fig. 17(c)) decrease as g_{c} increases. This is consistent with the fact that the stronger the stimulus (i.e., the coupling in this condition) to the neuron, the more spikes there are within the burst.
Recently, the spike and burst synchronization and the transition from the burst synchronization to the spike synchronization have been widely investigated.^{[17– 22]} In this paper, the Chay model shows that period1 bursting and period2 bursting coexist within a parameter range, and the coexisting behavior lies between monostable period1 and period2 bursting. As V_{C} increases, the attraction domain of period1 bursting increases while that of the coexisting period2 bursting decreases. Various synchronization states, including complete synchronization, spike phase synchronization, burst phase synchronization, lag synchronization, and nonsynchronization, are simulated in a network constructed by two bidirectionally coupled Chay model neurons with coexisting period1 and period2 bursting. The burst phase synchronization is identified as chaos synchronization and spike phase synchronization as synchronization of periodicity.
With increasing coupling strength g_{c}, the coupled neurons can achieve complete synchronization states via complex processes. The process depends on the parameter values within the coexisting parameter region and on the initial values. When the coexisting behavior is near period1 bursting, the process can be divided into four parts. Part 1 is the nonsynchronization. Part 2 represents the transitions between burst phase and spike phase synchronization. The processes of part 1 and part 2 are insensitive to initial values. Part 3 represents the transitions between spike phase and complete synchronization. The transit times and range of g_{c} depend on initial values. Part 4 represents complete synchronization. Most initial values lead to complete synchronization of period1 bursting and fewer initial values lead to complete synchronization of period2 bursting when g_{c} is large enough. When the coexisting behavior is near the period2 bursting, three parts similar to parts 2, 3, and 4 are simulated. Most initial values lead to complete synchronization of period2 bursting and fewer initial values lead to complete synchronization of period1 bursting when g_{c} is large enough. The increase of the probability that CS of period1 bursting will be achieve is, to a certain extent, consistent with the increase of attraction domain of period1 bursting. The dependence of the transitions on the parameter values is consistent with the results of a previous study.^{[37]} The synchronization transitions of this paper show more complex characteristics than those in previous works.^{[31– 38]} The most important finding of this paper is the sensitivity of the process to the initial values, especially for those transitions between spike phase and complete synchronization and for the complete synchronization. The dependence of the process on initial values was given little attention in previous works.^{[31– 38]} The synchronization processes that are complex, show a dependence on the initial values and parameters, contain transitions between spike and burst synchronization and between phase and complete synchronization, and chaos synchronization are all helpful to understand the complex spatiotemporal behaviors of networks^{[45– 51]} and synchronization transitions in the nervous system.^{[31– 38]}
The burst synchronization was found to be a chaos synchronization. This shows that chaos synchronization can be caused by cooperation between neurons with periodic behaviors. The transition between spike and burst phase synchronization states also illustrate the transition between chaos synchronization and synchronization of periodic behaviors. When the coexisting behavior is near monostable period1 bursting, we simulate the nonsynchronization of bursting pattern at small coupling strength. The behavior of the nonsynchronization state jumps between the coexisting attractors, which is similar to the asynchronous synchronization reported in a previous work.^{[31]}
The number of spikes per burst of noncomplete synchronous bursting increases with increasing coupling strength. This shows that bursting with multiple spikes per burst can be caused by the cooperation between neurons with one or two spikes per burst. Physiologically, bursting patterns with multiple spikes per burst can often be recorded by neurons of the central nervous system.^{[12, – 16]} The results in this paper indicate that the physiologically recorded bursting with multiple spikes per burst may originate from cooperation between neurons.
From the results of the related investigations^{[31– 38]} and this paper, we conclude that the transitions from nonsynchronization to complete synchronization with increasing coupling strength, which depends on parameter values and initial values, can be simulated in the coupled systems. These systems show multistable or coexisting behaviors when isolated. Various detailed characteristics are shown for specific networks with coexisting behaviors.^{[31– 38]} In the future, more investigations should be performed to reveal more general dynamics of network with coexisting behaviors. The structure of the attraction domain may be helpful to understand the complex transitions. In Ref. [37], the attraction domain is continuous and no chaos synchronization appears. In the present paper, the attraction domain shows a discrete structure or possible fractal structure. In addition, the transitions in Ref. [37] looks simpler than those in the present paper. Combining the trajectory of the coupled system and the attraction domain may be helpful to further identify the dynamics of the complex transitions.
1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 
