Transitions from chimeras to coherence: An analytical approach by means of the coherent stability function

TL;DR

This paper develops an analytical method using the coherent stability function to study the transition from chimera states to coherence in coupled oscillators, validated through numerical simulations with leech neuron models.

Contribution

It introduces the coherent stability function approach for analyzing transitions from chimeras to coherence, extending stability analysis tools to complex dynamical states.

Findings

Analytical conditions for chimera to coherence transition derived.

Numerical validation using leech neuron models confirms the approach.

Observation of complete synchronization at high coupling strengths.

Abstract

Chimera states have been a vibrant subject of research in the recent past, but the analytical treatment of transitions from chimeras to coherent states remains a challenge. Here we analytically derive the necessary conditions for this transition by means of the coherent stability function approach, which is akin to the master stability function approach that is traditionally used to study the stability of synchronization in coupled oscillators. In chimera states, there is typically at least one group of oscillators that evolves in a drifting, random manner, while other groups of oscillators follow a smoother, more coherent profile. In the coherent state, there thus exists a smooth functional relationship between the oscillators. This lays the foundation for the coherent stability function approach, where we determine the stability of the coherent state. We subsequently test the analytical prediction numerically by calculating the strength of incoherence during the transition point. We use leech neurons, which exhibit a coexistence of chaotic and periodic tonic spiking depending on initial conditions, coupled via non-local electrical synapses, to demonstrate our approach. We systematically explore various dynamical states with the focus on the transitions between chimeras and coherence, fully confirming the validity of the coherent stability function. We also observe complete synchronization for higher values of the coupling strength, which we verify by the master stability function approach.

Tables

Table 1. Table 1: Parameter values for the Leech neuron model ( 3 ) r7 .

Parameter	Value	Parameter	value
${\overline{g}}_{K2}$	$30ns$	$E_1$	$-0.046V$
${\overline{g}}_{Na}$	$200ns$	$C$	$0.5nF$
$g_1$	$8ns$	${\tau }_{K2}$	$0.25sec$
$E_K$	$-0.07V$	${\tau }_{Na}$	$0.0405sec$
$E_{Na}$	$0.045V$	$V_{K2}^{shift}$	$-0.025361V$
$A_1$	$-150$	$B_1$	$0.0305$
$A_2$	$-83$	$B_2$	$0.018$
$A_3$	$500$	$B_3$	$0.0333$

Equations

$$\begin{array}[]{l}\dot{\bf x}_{i}=F({\bf x}_{i})+\dfrac{\epsilon}{d(i)}\sum\limits_{j=1}^{N}\mathscr{A}_{ij}\Gamma({\bf x}_{j}-{\bf x}_{i}),\end{array}$$

$$\begin{array}[]{lll}\mathscr{L}_{ij}&=-1\hskip 10.0pt{\text{for}\leavevmode\nobreak\ }0<|i-j|\leq P,\\ &=2P\hskip 10.0pt{\text{for}\leavevmode\nobreak\ }i=j,\\ &=0\hskip 18.0pt{\text{otherwise}}.\end{array}$$

$$\begin{array}[]{l}\dot{V}=-\frac{1}{C}\big{[}\bar{g}_{K2}m^{2}_{K2}(V-E_{K})+g_{1}(V-E_{1})+\\[2.0pt] \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \bar{g}_{Na}\big{(}f(A_{1},B_{1},V)\big{)}^{3}h_{Na}(V-E_{Na})\big{]},\\[5.0pt] \dot{m}_{K2}=\frac{f(A_{2},B_{2}+V^{shift}_{K2},V)-m_{K2}}{\tau_{K2}},\\[5.0pt] \dot{h}_{Na}=\frac{f(A_{3},B_{3},V)-h_{Na}}{\tau_{Na}},\end{array}$$

$$\begin{array}[]{l}f(A,B,V)=\dfrac{1}{1+e^{A(B+V)}}.\\ \end{array}$$

$$\begin{array}[]{l}\dot{V}_{i}=-\frac{1}{C}\big{[}\bar{g}_{K2}m^{2}_{K2_{i}}(V_{i}-E_{K})+g_{1}(V_{i}-E_{1})+\\ \bar{g}_{Na}\big{(}f(A_{1},B_{1},V_{i})\big{)}^{3}h_{Na_{i}}(V_{i}-E_{Na})\big{]}+\frac{\epsilon}{2P}\sum\limits_{\begin{subarray}{c}j=i-P\\ j\neq i\end{subarray}}^{i+P}(V_{j}-V_{i}),\\ \dot{m}_{K2_{i}}=\dfrac{f(A_{2},B_{2}+V^{shift}_{K2},V_{i})-m_{K2_{i}}}{\tau_{K2}},\\[7.0pt] \dot{h}_{Na_{i}}=\dfrac{f(A_{3},B_{3},V_{i})-h_{Na_{i}}}{\tau_{Na}},\end{array}$$

$$\begin{array}[]{l}L_{i}(t)=\bigg{|}\dfrac{1}{2\delta}\sum\limits_{|i-k|\leq\delta}e^{\phi_{k}(t)\sqrt{-1}}\bigg{|},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ i=1,2,\cdots,N,\end{array}$$

$$\begin{array}[]{l}\phi_{k}(t)=\arctan\Big{(}\frac{y_{k}(t)}{x_{k}(t)}\Big{)}.\end{array}$$

$$\begin{array}[]{l}{\bf x}_{i}=\Phi({\bf x}_{1},{\bf x}_{2},\dots,{\bf x}_{i-1},{\bf x}_{i+1},\dots,{\bf x}_{N}).\end{array}$$

$$\begin{array}[]{l}\begin{array}[]{lll}\dot{\bf y}_{i}=F({\bf y}_{i})+\dfrac{\epsilon}{d(i)}\sum\limits_{j=1}^{N}\mathscr{A}_{ij}\Gamma({\bf x}_{j}-{\bf y}_{i}).\end{array}\end{array}$$

$$\begin{array}[]{l}\dot{\bf z}_{i}=F({\bf x}_{i}+{\bf z}_{i})-F({\bf x}_{i})-\epsilon\Gamma{\bf z}_{i},\hskip 10.0pti=1,2,\dots,N.\end{array}$$

$$\begin{array}[]{l}\dot{\bf z}_{i}=[JF({\bf x}_{i})-\epsilon\Gamma]{\bf z}_{i},\hskip 20.0pti=1,2,\dots,N.\end{array}$$

$$\begin{array}[]{l}\left\lVert{\bf z}_{i}\right\rVert\rightarrow 0\Leftrightarrow\leavevmode\nobreak\ \left\lVert{\bf x}_{i}-{\bf y}_{i}\right\rVert\rightarrow 0\mbox{ as }t\rightarrow\infty.\end{array}$$

$$\begin{array}[]{l}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \dot{z}_{i}^{V}=-\frac{1}{C}\Big{[}\bar{g}_{K2}m^{2}_{K2_{i}}+g_{1}+\bar{g}_{Na}\big{(}f(A_{1},B_{1},V_{i})\big{)}^{3}h_{Na_{i}}+3\bar{g}_{Na}\big{(}f(A_{1},B_{1},V_{i})\big{)}^{2}f_{V}(A_{1},B_{1},V_{i})(V_{i}-E_{1})h_{Na_{i}}\Big{]}z_{i}^{V}\\[5.0pt] \hskip 36.0pt-\frac{1}{C}\Big{[}2\bar{g}_{K2}m_{K2_{i}}(V_{i}-E_{K})z_{i}^{m_{K2}}+\bar{g}_{Na}\big{(}f(A_{1},B_{1},V_{i})\big{)}^{3}(V_{i}-E_{Na})z_{i}^{h_{Na}}\Big{]}-\epsilon z_{i}^{V},\\[8.0pt] \dot{z}_{i}^{m_{K2}}=\dfrac{f_{V}(A_{2},B_{2}+V^{shift}_{K2},V_{i})z_{i}^{V}-z_{i}^{m_{K2}}}{\tau_{K2}},\\[8.0pt] \leavevmode\nobreak\ \dot{z}_{i}^{h_{Na}}=\dfrac{f_{V}(A_{3},B_{3},V_{i})z_{i}^{V}-z_{i}^{h_{Na}}}{\tau_{Na}},\hskip 50.0pti=1,2,\dots,N.\end{array}$$

$$\begin{array}[]{l}\sigma(m)=\left\langle\sqrt{\frac{1}{n}\sum\limits_{j=n(m-1)+1}^{nm}\big{[}z_{j}-\langle z\rangle\big{]}^{2}}\right\rangle_{t},\end{array}$$

$$\begin{array}[]{l}V_{\mbox{SI}}=1-\frac{1}{M}\sum\limits_{m=1}^{M}\Theta(\delta-\sigma(m)),\end{array}$$

$$\begin{array}[]{l}\left\lVert F({\bf x})-F({\bf y})\right\rVert\leq M\left\lVert{\bf x}-{\bf y}\right\rVert.\end{array}$$

$$\begin{array}[]{lll}G(t)=\frac{1}{2}\sum\limits_{i=1}^{N}{\bf z}_{i}^{tr}{\bf z}_{i},\end{array}$$

$$\begin{array}[]{lll}\dot{G}(t)=\sum\limits_{i=1}^{N}{\bf z}_{i}^{tr}\dot{{\bf z}}_{i}\\ \hskip 22.0pt=\sum\limits_{i=1}^{N}{\bf z}_{i}^{tr}[F({\bf x}_{i}+{\bf z}_{i})-F({\bf x}_{i})]-\epsilon\sum\limits_{i=1}^{N}{\bf z}_{i}^{tr}\Gamma{\bf z}_{i}\end{array}$$

$$\begin{array}[]{lll}\dot{G}(t)\leq\sum\limits_{i=1}^{N}M{\bf z}_{i}^{tr}\dot{{\bf z}}_{i}-\epsilon\sum\limits_{i=1}^{N}{\bf z}_{i}^{tr}\Gamma{\bf z}_{i}=\big{[}M-\epsilon\lambda(I_{N}\otimes\Gamma)\big{]}{\bf Z}^{tr}{\bf Z},\end{array}$$

$$\begin{array}[]{lll}\epsilon>\frac{M}{\lambda[\Gamma]}.\end{array}$$

$$\begin{array}[]{l}\delta\dot{\bf x}_{i}=JF({\bf x}_{0})\delta{\bf x}_{i}-\dfrac{\epsilon}{2P}\sum\limits_{j=1}^{N}\mathscr{L}_{ij}\Gamma\delta{\bf x}_{j}.\end{array}$$

$$\begin{array}[]{l}\dot{\xi}_{i}=[JF({\bf x}_{0})-\frac{\epsilon}{2P}\gamma_{i}\Gamma]\xi_{i},\hskip 10.0pti=1,2,\dots,N.\end{array}$$

$$\begin{array}[]{l}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \dot{\xi}_{i}^{V}=-\frac{1}{C}\Big{[}\bar{g}_{K2}m^{2}_{K2}+g_{1}+\bar{g}_{Na}\big{(}f(A_{1},B_{1},V)\big{)}^{3}h_{Na}+3\bar{g}_{Na}\big{(}f(A_{1},B_{1},V)\big{)}^{2}f_{V}(A_{1},B_{1},V)(V-E_{1})h_{Na}\Big{]}\xi_{i}^{V}\\[5.0pt] \hskip 35.0pt-\frac{1}{C}\Big{[}2\bar{g}_{K2}m_{K2}(V-E_{K})\xi_{i}^{m_{K2}}+\bar{g}_{Na}\big{(}f(A_{1},B_{1},V)\big{)}^{3}(V-E_{Na})\xi_{i}^{h_{Na}}\Big{]}-\dfrac{2\epsilon}{P}\sum\limits_{l=1}^{P}\sin^{2}\big{(}\frac{\pi l(i-1)}{N}\big{)}\leavevmode\nobreak\ \xi_{i}^{V},\\[8.0pt] \dot{\xi}_{i}^{m_{K2}}=\dfrac{f_{V}(A_{2},B_{2}+V^{shift}_{K2},V)\xi_{i}^{V}-\xi_{i}^{m_{K2}}}{\tau_{K2}},\\[8.0pt] \leavevmode\nobreak\ \dot{\xi}_{i}^{h_{Na}}=\dfrac{f_{V}(A_{3},B_{3},V)\xi_{i}^{V}-\xi_{i}^{h_{Na}}}{\tau_{Na}},\hskip 50.0pti=2,3,\dots,N.\end{array}$$

$$\begin{array}[]{l}E_{\mbox{CS}}=\Bigg{\langle}\sqrt{\dfrac{1}{N(N-1)}\sum\limits_{\begin{subarray}{c}i,j=1\\ j\neq i\end{subarray}}^{N}\Big{[}(V_{i}-V_{j})^{2}+(m_{K2_{i}}-m_{K2_{j}})^{2}+(h_{Na_{i}}-h_{Na_{j}})^{2}\Big{]}}\Bigg{\rangle}_{t}\end{array}$$

Full text

I Introduction

Chimera states are emergent spatiotemporal patterns consisting of the coexistence of coherent and drifting domains in a network of identical oscillators. They were first observed in non-locally coupled identical phase oscillators by Kuramoto and Battogtokh r12 . Later these patterns was named chimera states by Abrams and Strogatz r13 . With the invention of these states, research on the emergent behaviors of a network of coupled oscillators has been reanimated. Chimera states have been studied extensively in recent years r15a ; chimera_r1 ; r14a ; chimera_r2 . They have been investigated in different types of oscillators, such as phase oscillators r18 , chaotic and periodic maps r19 , pendulum oscillators r20 , neuronal models r15a , mechanical systems r21 , and chemical systems r22 . Since their discovery, several types of chimera states have been identified based on different schemes of spatio-temporal patterns. An amplitude chimera state amplt_cmra is one in which in the incoherent subpopulation, the amplitudes of the oscillators differ, but the phase velocity of the entire network remains the same. A spiral wave chimera r31 ; epjst is a two-dimensional chimera consisting of a randomizing desynchronized core surrounded by phase-locked oscillators. Chimera death is a partially coherent inhomogeneous pattern consisting of spatially coherent and incoherent oscillation death r25 . In a stationary chimera state, the position of coherent and incoherent clusters is fixed. However, in a non-stationary chimera state, the position of coherent and incoherent clusters changes over time r40 . There are other types of non-stationary chimera states, such as breathing chimera r18 , imperfect chimera r20 , traveling chimera r46 , multi-headed chimera r29 and alternating chimera r30 . Recently, a new type of nonstationary chimera state called a spike chimera state spike has also been discovered, where the chimera and fully coherent states appear alternately with respect to time in a neuronal hypernetwork. In the context of chimera states, the underlying network architecture and the coupling function also play an important role. With linear diffusive coupling, the chimera states emerge only for non-local r13 or global chimera_global network architecture, while for non-linear chemical synaptic coupling, chimera states can emerge irrespective of local, non-local, global r39 , and even unidirectional local and non-local coupling spike . Chimera states have also has been observed in many real biological and physical systems, including ecology cmra_eco1 ; r29 , SQUID metamaterial cmra_SQUID , and quantum systems cmra_quantum .

The nervous system is a complex biological process that consists of billions of neurons. Neurons as systems with strong non-linear dynamical properties can generate various complex activities from a rest state to tonic spiking and bursting states r7 ; r9 . Among these non-linear complex activities, a ubiquitous observed pattern is tonic spiking, which correlates to different functions of the nervous systems, such as the coding of sensory information, information processing, memory formation, attention, and motor control r8 . There are many mathematical models that effort to illustrate and investigate different behaviors of neurons r1 . One effective way to gain a better understanding of the behaviors of the entire nervous system is to investigate them in the network of coupled neurons r2 . Such investigations, which provide the advantage of considering the interaction among all particles of the network, are significant and useful tools to perceive the potential mechanism of neuronal diseases r6 . Several mathematical neuronal models can exhibit the coexistence of tonic spiking, bursting, and rest states, which has been studied extensively r2 . In a dynamical system, the coexistence of several attractors in phase space is referred to as multi-stability, where each attractor is determined by initial conditions. Lencher et al. r10 exhibited bistability in a Leech neuron via experiments. The isolated Leech neuronal model exhibits bistable states, consisting of one chaotic and one periodic states, depending on the initial conditions r7 .

Previous studies have shown that a chimera state is closely related to various types of brain phenomena, such as hemisphere sleep in animals, epilepsy, Parkinson’s, Alzheimer’s, and schizophrenia r32 ; r33 ; r34 . Chimera states also play an important role in neuroscience and cardiac fibrillation r35 . Andrzejak et al. r33 exhibited an abrupt conversion of the coexistence of coherent and incoherent oscillators into completely coherent oscillators in an epileptic seizure in humans. It has been shown that spiral wave chimera states appear in the contractions of heart cells during ventricular fibrillation, which is formed on the surface of a sphere chimera_r1 . Chimera states have also been observed in unihemispheric slow-wave sleep of some aquatic mammals and a few migratory birds cmra_unih1 ; cmra_unih2 . During slow-wave sleep, one-half of the brain is awake and the neurons in that portion are desynchronized, while the other half is at rest and the corresponding neurons are synchronized. Chimera states have also been investigated in neuronal networks, such as the Hindmarsh-Rose r37 ; r37a ; r39 and FitzHugh-Nagumo network r38 ; r38a . Vüllings et al. r37 considered a two-dimensional network of Hindmarsh-Rose neurons by a rotating coupling matrix, representatives of excitability type I, and they obtained different multi-headed chimera patterns. Omelchenko et al. r38 constructed a network of identical non-locally coupled Fitzhugh-Nagumo oscillators and found new types of chimera states.

Motivated by the above, in this paper, we investigate chimera states in an identical Leech neuronal network that interact through non-local electrical synapses. The non-local coupling emerges in several real-world applications ranging from neuronal networks neuronal_nl , more particularly in ocular dominance strips ocular_nl , to chemical oscillators r12 and Josephson function arrays jj_nl . In particular, the appearance of non-stationary chimera patterns is observed, which have an important role in neural networks and need more attention r41 . The most interesting part of our work is the analysis of the instability of chimera states towards coherent states.

In our observed non-stationary chimera states, the entire network splits into two groups, one with an incoherent domain and the other with a coherent domain. Contrary to these chimera states, all the oscillators in the coherent state follow a smooth coherent profile, i.e., they follow an amplitude-locked oscillation. From this, we claim that there may exist a smooth functional relationship among the oscillators. That is, the network of oscillators follows the generalized synchronization (GS) state. This is a typical coordination of the rhythms among the coupled oscillators. This type of synchronization was first noticed by Rulkov et al. ausys_3 between two unidirectionally coupled systems. The auxiliary system approach ausys_1 has been widely used to detect and analyze the stability of GS states. At first, this approach was proposed in a drive-response system ausys_1 . The necessary and sufficient conditions for GS based on the auxiliary system approach were later illustrated in Ref. ausys_4 .

Later, this seminal approach was extended to bidirectionally coupled systems and complex networks by introducing a system consisting of an auxiliary node for each original node ausys_2 . In Ref. hung , the authors uncovered GS in a bidirectional scale-free network of identical chaotic oscillators using an auxiliary-system approach. They noticed that the GS starts from the hubs of the network and then spreads throughout the whole network upon increasing the coupling strength. The occurrence and development of GS on various mutually coupled complex networks, including scale-free networks, small-world networks, random networks, and modular networks are also investigated guan . It has been shown that, this phenomenon generally takes place in such networks for both coupled identical and nonidentical oscillators. For coupled identical oscillators, they found a typical path from non-synchronization to complete synchronization (CS) through GS, i.e. nonsynchronization $\rightarrow$ GS $\rightarrow$ CS. In Ref. chen , the authors investigated GS in two typical complex dynamical networks, namely small-world networks and scale-free networks, in terms of an impulsive control strategy, and they analyzed how the development of GS depends on impulsive control gains. By applying the auxiliary-system approach, they demonstrated their theoretical analysis. GS was also investigated via the adaptive bidirectional couplings strategy liu . Adopting the auxiliary-system approach and the Lyapunov function method, the authors proved that for any given initial coupling strengths, GS can take place in complex networks consisting of nonidentical dynamical systems. Subsequently, this state was further characterized by analytical and numerical experiments in mutually coupled networks gs_njp . Recently, Moskalenko et al.olga compared the applicability of the auxiliary-system approach in unidirectional and bidirectional coupled systems.

Hence, we analyze the stability of the coherent state in terms of GS by the auxiliary system approach. The instability of chimera states indicates the stability of coherent states, while the stability of coherent states implies the instability of chimera states. Thus, we could identify the critical transition point of the chimera to the coherent state, analytically obtained by the coherent stability function approach. The effect of initial conditions is also studied by the calculation of basin stability measurement. For a higher value of the coupling strength, complete synchronization occurs in the non-locally coupled neuronal network, and it is verified by the master stability function approach.

II Mathematical model

Generally, a network of $N$ identical coupled nodes can be written as

[TABLE]

where $i=1,2,\dots,N$; ${\bf x}_{i}$ is the state variable of node $i$. Here each isolated node modeled by a $d-$dimensional dynamical system $\dot{\bf x}=F({\bf x})$, which is a continuous vector valued function. $\epsilon$ is the coupling strength, which modulates the amount of the information propagation among the nodes, $d(i)$ is the degree of node $i$, and $\Gamma$ is the inner coupling matrix that determines which node’s variable is used in the coupling. $\mathscr{A}$ is the adjacency matrix of the network. It represents the topology of the connectivity in the network, where $\mathscr{A}_{ij}=1$ if the $i^{th}$ node and the $j^{th}$ node are connected and zero otherwise. If $\mathscr{L}$ is the corresponding zero row sum Laplacian matrix of the network, then $\mathscr{L}_{ij}=-\mathscr{A}_{ij}$ if $i\neq j$ and $\mathscr{L}_{ii}=\sum_{j=1}^{N}\mathscr{A}_{ij}$. The chimera state is the broken of symmetry among the symmetrically coupled identical oscillators; only their initial conditions are different. Therefore, the underlying network is considered to be a regular network that is either local or non-local, or global. For a non-locally coupled network, the Laplacian matrix can be written as

[TABLE]

Therefore, the degree of the $i^{th}$ node is $d(i)=2P$. For $P=1$, it is a locally coupled network and for $1<P<\frac{N}{2}$, it is non-local. The Laplacian matrix $\mathscr{L}$ is a real symmetric matrix, so all it’s eigenvalues are real, and it is orthogonally diagonalizable by the basis of eigenvectors $U$. Moreover, $\mathscr{L}$ is a circulant matrix. It’s $N$ eigenvalues are $\sum_{l=1}^{N}\mathscr{L}_{1l}\omega_{j}^{l-1}$, where \omega_{j}=\exp\Big{(}\frac{2\pi j\sqrt{-1}}{N}\Big{)} and $j=0,1,\dots,N-1$, which yields all $N$ eigenvalues as 4\sum_{l=1}^{P}\sin^{2}\big{(}\frac{\pi lj}{N}\big{)}, where $j=0,1,\dots,N-1$. For $j=0$, we get it’s unique zero eigenvalues, and others are non-zeroes.

III Leech neuronal network

A Leech neuron exhibits bi-stability of chaotic and periodic attractors. It can exhibit different neuronal activities, such as tonic spiking and bursting. Cymbalyuk et al. r7 simplified the Leech neuron in such a way that it depends on the dynamics of the potassium current $I_{K2}$ and the sodium current $I_{Na}$. This model is described by a three-dimensional differential equation

[TABLE]

where the variables $V$, $m_{K2}$, and $h_{Na}$ describe the membrane potential, the activation of $I_{K2}$ and the inactivation of $I_{Na}$ respectively. The parameters: $\bar{g}_{K2}$, $\bar{g}_{Na}$, and $g_{1}$ demonstrate the maximum conductance of $I_{K2}$, $I_{Na}$ and the conductance of the leak current, respectively. $E_{k}$, $E_{Na}$, and $E_{1}$ refer to the reversal potential of $K^{+}$, $Na^{+}$, and the leak current, respectively. $C$ is the membrane capacitance, $\tau_{K2}$ and $\tau_{Na}$ illustrate the time constants of activation of $I_{K2}$ and inactivation of $I_{Na}$, respectively and $V^{shift}_{K2}$ characterizes the shift of the membrane potential of half-activation of $I_{K2}$ from its canonical value. $f(A,B,C)$ is the Boltzmann function r7 ,

[TABLE]

The values of the constant parameters that are used in this paper are represented in Table 1.

Initially for $V^{shift}_{K2}=-0.026\leavevmode\nobreak\ V$ in system (3), both tonic spiking modes are represented by two periodic attractors. As the bifurcation parameter $V^{shift}_{K2}$ increases, through a cascade of period-doubling bifurcation, one periodic attractor becomes chaotic. Finally, system (3) exhibits the coexistence of a chaotic attractor, and a limit cycle attractor for $V^{shift}_{K2}=-0.025361\leavevmode\nobreak\ V$. These two bistable attractors of a single Leech neuronal oscillator are depicted in Fig. 1(a). One is a periodic orbit, which is plotted by a red (gray) curve, and the another is a chaotic attractor, shown by a blue (black) curve. Figure 1(b) demonstrates their basins of attraction, the blue (black) region is the basin of attraction of the chaotic state, while the red (gray) region depicts the basin of attraction of the period-1 limit cycle state. This shows that their basins of attraction are almost well-separated. To get a clearer view of the well-separated basin of attraction, we plot it on the $h_{Na}=0.5$ plane in Fig. 1(c).

To construct the neuronal network, a one-dimensional ring of $N$ identical oscillators of the system (3) is considered, in which each oscillator is coupled to its $2P$ nearest neighbors. Then the evolution equations of the dynamical network are given as follows:

[TABLE]

where $i=1,2,\dots,N$ is the index of the $i^{th}$ oscillator, and the parameter $\epsilon$ is the strength of coupling. For numerical simulation, the number of oscillators $N$ and the number of neighbors in each side $P$ are considered to be fixed at $N=200$ and $P=20$. Throughout this work, we integrate Eq. (5) using the fourth-order Runge-Kutta-Fehlberg method with a fixed time step $dt=0.001$ for $3\times 10^{5}$ time iterations. We check that our observations are not specific to either the chosen integrator or the numerical accuracy. In our work, we consider three types of ‘V’ shaped initial conditions (with small random fluctuations) as follows.

PIs: ‘V’ shaped initial conditions for each oscillator are taken from the basin of attraction of the periodic state, i.e., the dotted line in Fig. 1(c) as $V_{i}(0)=0.1-0.2\times\frac{i-1}{\frac{N}{2}-1},\leavevmode\nobreak\ m_{K2_{i}}(0)=-1.25V_{i}+0.125,\leavevmode\nobreak\ h_{Na_{i}}(0)=0.5$ for $1\leq i\leq\frac{N}{2}$; and $V_{i}(0)=-0.1+0.2\times\frac{i-\frac{N}{2}-1}{\frac{N}{2}-1},\leavevmode\nobreak\ m_{K2_{i}}(0)=1.25V_{i}+0.375,\leavevmode\nobreak\ h_{Na_{i}}(0)=0.5$ for $\frac{N}{2}+1\leq i\leq N$. 2. 2.

CIs: ‘V’ shaped initial conditions for each oscillator are taken from the basin of attraction of the chaotic state, i.e., the dashed line in Fig. 1(c) as $V_{i}(0)=0.1-0.135\times\frac{i-1}{\frac{N}{2}-1},\leavevmode\nobreak\ m_{K2_{i}}(0)=-0.4444V_{i}+0.2344,\leavevmode\nobreak\ h_{Na_{i}}(0)=0.5$ for $1\leq i\leq\frac{N}{2}$; and $V_{i}(0)=-0.035+0.135\times\frac{i-\frac{N}{2}-1}{\frac{N}{2}-1},\leavevmode\nobreak\ m_{K2_{i}}(0)=0.0741V_{i}+0.2576,\leavevmode\nobreak\ h_{Na_{i}}(0)=0.5$ for $\frac{N}{2}+1\leq i\leq N$. 3. 3.

MIs: ‘V’ shaped mixed type initial conditions are taken from the entire phase space, i.e., the solid line in Fig. 1(c) as $V_{i}(0)=0.1-0.135\times\frac{i-1}{\frac{N}{2}-1},\leavevmode\nobreak\ m_{K2_{i}}(0)=-0.8519V_{i}+0.27019,\leavevmode\nobreak\ h_{Na_{i}}(0)=0.5$ for $1\leq i\leq\frac{N}{2}$; and $V_{i}(0)=-0.035+0.135\times\frac{i-\frac{N}{2}-1}{\frac{N}{2}-1},\leavevmode\nobreak\ m_{K2_{i}}(0)=0.8519V_{i}+0.3298,\leavevmode\nobreak\ h_{Na_{i}}(0)=0.5$ for $\frac{N}{2}+1\leq i\leq N$.

In the next section, our main emphasis will be to investigate the emergence of chimera states and the effects of these three types of initial conditions on them.

IV Emergence of chimera states in a Leech neuronal network

We first investigate the emergence of chimera states in the non-locally coupled Leech neuronal network, Eq. (5), by varying the interaction strength $\epsilon$. The first, second, and third columns of Fig. 2 respectively show the snapshot, two-dimensional phase space in the $(m_{k2_{i}},V_{i})$ plane, and local order parameter $L_{i}$ (which is defined later). In the top panel of Fig. 2, initial conditions are taken from PIs, with coupling strength $\epsilon=0.2$. The scatter population in the right side of Fig. 2(a) indicates that the corresponding neurons are drifting relative to the amplitude locked oscillators. Here, the neurons from $i=1$ to $i=116$ spike coherently, while the remaining neurons spike with different frequencies. From the intrinsic periodic state, all the neurons exhibit the quasi-periodic state depicted in Fig. 2(b), which is lying near that periodic state. Here, the oscillators corresponding to the yellow curves are in the incoherent state, and all the oscillators from the coherent region follow the blue trajectories. The yellow trajectories of the incoherent oscillators are more spread out in the phase space compared to the those of coherent oscillators. The middle panel of Fig. 2 is drawn for CIs initial conditions, where $\epsilon$ is fixed at $0.4$. Figure 2(d) represents the snapshot, where the oscillators from $i=103$ to $i=171$ maintain the coherent state. Corresponding two-dimensional phase space is plotted in Fig. 2(e), where the blue trajectories are for the coherent group of oscillators, while the yellow curves represent the trajectories of the incoherent group of oscillators. Here, CIs-type initial conditions are chosen. But due to the interaction, the trajectories of a few neurons spike quasi-periodically, while the rest of the neurons remain in the chaotic state, and one neuron is oscillating from the chaotic to the quasi-periodic state (deep yellow). For MIs initial conditions, the snapshot of the chimera states is drawn in Fig. 2(g) for $\epsilon=0.3$. In one domain, the oscillators are desynchronized, and the another two domains, they are in coherent motion. As in the previous case, all the coherent oscillators follow a quasi-periodic trajectory (blue), but the incoherent oscillators are divided into two groups, where one of them is in the chaotic state while the other follows the quasi-periodic state.

By using the notion of a local order parameter, we characterize the spatial coexistence of coherent-incoherent patterns of the chimera states. To determine the coherent-incoherent clusters and chimera states, the concept of a local order parameter is often used r46 . The local order parameter demonstrates the neighbor ordering of the neuronal network and quantifies the coherent and incoherent clusters. This parameter is defined as

[TABLE]

where $\delta$ denotes the window size of the $i^{th}$ neuron used for the spatial average; here we consider $\delta=3$. The geometric phase of the $k^{th}$ neuron is determined by the formula

[TABLE]

The color bars in the third column of Fig. (2) represent the variation of $L_{i}(t)$ for three different types of initial conditions, which corroborate with the corresponding snapshot plots in the first column. White region corresponds to the values of $L_{i}(t)$ equal to $1$, which indicate the maximum coherency of the $i^{th}$ neuron, and it belongs to the coherent population. However, for the $i^{th}$ neuron that belongs to the incoherent domain of the chimera state, $L_{i}$ is obtained at less than $1$, which is shown as other colors in the color bars.

When the coupling strength gets higher values than $\epsilon=8.0$, the chimera states disappear, and the network of neurons tends to oscillate coherently. Figure 3(a) shows the snapshot of the coherent state, taking the ‘V’ shaped periodic initial conditions for the coupling strength $\epsilon=10.0$. This depicts a smooth coherent manifold among the oscillators. In this coherent state, all the oscillators follow a chaotic motion near its intrinsic periodic orbit, drawn in Fig. 3(b). At the chimera states, the neuronal network (5) possesses different types of spiking resting behavior as we observe from the phase space plots of Fig. (2) (second column). But at the coherent state, the neuronal network exhibits the same chaotic trajectory regardless of the other two types of initial conditions (checked separately, but the results are not presented here). The corresponding local order parameter is plotted in Fig. 3(c). Throughout the time interval $t\in[200,220]$, $L_{i}(t)\simeq 1$ for all $i=1,2,\dots,N$. This describes the local coherency of all the oscillators. From these results, it is clear that the oscillators in the coherent states are amplitude-locked, and all of them oscillate at an identical instantaneous frequency and smooth amplitude profile.

V Analytical study of the chimera to coherent transition

In chimera states, symmetrically coupled identical oscillators are split into two regions; one is a coherent region where the phases of the oscillators are locked, and the other is an incoherent region where the oscillators are in a drifting motion. So, in the chimera states, there exists a subpopulation of oscillators that oscillate in a drifting manner, which is random. Any oscillator from this incoherent domain, could not obey a smooth transformation of the remaining oscillators. Therefore, in the whole network, there could not exist a smooth transformation function of a particular oscillator in terms of all other oscillators. But for coherent states, the amplitudes and phases are locked for all the oscillators. At every time, their snapshots show a smooth profile. From this fact, we claim that for all $i=1,2,\dots,N$, there exists a smooth transformation $\Phi:\mathbb{R}^{d(N-1)}\rightarrow\mathbb{R}^{d}$ such that

[TABLE]

In other words, we assert that in the coherent state, the network of oscillators follows generalized synchronization, which is a general type of correlation among the coupled oscillators.

If there exists such a smooth mapping $\Phi$, and a manifold \mathcal{M}=\big{\{}({\bf x}_{1},{\bf x}_{2},\dots,{\bf x}_{N})\leavevmode\nobreak\ \big{|}\leavevmode\nobreak\ {\bf x}_{i}=\Phi({\bf x}_{1},{\bf x}_{2},\dots,{\bf x}_{i-1},{\bf x}_{i+1},\dots,{\bf x}_{N})\big{\}}\subset\mathbb{R}^{dN}, such that the trajectories of network (1) with random initial conditions from its basin of attractions approach $\mathcal{M}$ as $t\rightarrow\infty$, then the network of oscillators is said to realize generalized synchronization. Here $\mathcal{M}$ is the generalized synchronization manifold, which is an invariant manifold, and it is attracting if the generalized synchronization state is stable. The existence of the smooth mapping $\Phi$ guarantees the functional relationship of ${\bf x}_{i}$ with ${\bf x}_{1},{\bf x}_{2},\dots,{\bf x}_{i-1},{\bf x}_{i+1},\dots,{\bf x}_{N}$.

To mathematically detect the occurrence of the generalized synchronization states, we use the auxiliary-system approach ausys_1 . It is a powerful mathematical tool to detect the above-mentioned functional relationship among coupled oscillators. This approach is applied to bidirectionally coupled dynamical networks by introducing an auxiliary node system corresponding to each original node ausys_2 . It can be described for the $i^{th}$ node as

[TABLE]

If ${\bf y}_{i}$ identically synchronizes with ${\bf x}_{i}$, then the GS of the network (1) is assumed to occur. Obviously, all the nodes show identical behavior with their auxiliary counterpart. The system (1) is said to achieve generalized synchronization if $\left\lVert{\bf x}_{i}(t)-{\bf y}_{i}(t)\right\rVert\rightarrow 0$ as $t\rightarrow\infty$, where $\left\lVert\cdot\right\rVert$ is the Euclidean norm. Considering the variation of the $i^{th}$ node ${\bf z}_{i}(t)={\bf y}_{i}(t)-{\bf x}_{i}(t)$, we can get the equation of motion of the error systems as

[TABLE]

For the local stability of the generalized synchronization manifold, we consider a small difference ${\bf z}_{i}(t)$ of the state variables of the original system and its auxiliary partner. Then we get the coherent stability equation as

[TABLE]

Here $J$ denotes the Jacobian of the function with respect to the arguments, and $({\bf x}_{1},{\bf x}_{2},\dots,{\bf x}_{N})$ is the state variable of the generalized synchronization manifold obeying Eq. (1). Now,

[TABLE]

Therefore, generalized synchronization of the network (1) is achieved if and only if the trivial solution of the above error system (11) is asymptotically stable. The linearized Eq. (11) will be solved in parallel to the original non-linear Eq. (1) for ${\bf x}_{i}$, which then allows us to compute all its Lyapunov exponents. The coherent stability function $\Lambda_{\mbox{GS}}$, the maximum of these Lyapunov exponents as a function of $P$ and $\epsilon$, actually gives the necessary condition for the stability of the chimera and coherent states. On the basis of the linear stability theory, the negative Lyapunov exponents can bring about the stability of the trivial solution of the error system. So, the maximum Lyapunov exponent, $\Lambda_{\mbox{GS}}$, of the system (11) should be negative for the local stability of the coherent states. Since the chimera and the coherent states are mutually exclusive and exhaustive, $\Lambda_{\mbox{GS}}>0$ implies the stability of the chimera states, while $\Lambda_{\mbox{GS}}<0$ indicates the stability of the coherent states. The coherent stability function allows one to quickly establish whether the coherent state is stable or not. In this way, it reveals the instability of the chimera states. Here, $\Lambda_{\mbox{GS}}>0$ is the necessary condition for the chimera states, but it is not sufficient, as it may also exhibit an incoherent state.

Finally, for the non-locally coupled Leech neuronal network, our required coherent stability equation can be written as

[TABLE]

Here, f_{V}(A,B,V)=\frac{\partial}{\partial V}f(A,B,V)=\frac{-A\exp[A(B+V)]}{\big{\{}1+\exp[A(B+V)]\big{\}}^{2}}, and the state variables $(V_{i},m_{K2_{i}},h_{Na_{i}})$ of the generalized synchronization are manifold dominated by Eq. (5). When the chimera states of the non-locally coupled system become unstable, the coherent states turn out to be stable. The maximum Lyapunov exponent, $\Lambda_{\mbox{GS}}$, of the above error system (13) becomes negative, and with the stable coherent state, all the oscillators evolve generalized synchronously.

To distinguish the chimera states from the other states, such as incoherent and coherent states, the time-averaged statistical measure of the strength of incoherence $(V_{\mbox{SI}})$ is widely used r34 . To compute the value of $V_{\mbox{SI}}$, we divide the total number of oscillators into an even number of $M$ bins, each of which is of equal length $n=\frac{N}{M}$. Then the local standard deviation $\sigma(m)$ is defined as

[TABLE]

where $m=1,2,\dots,M$; the difference of the potential state variables $V_{i}$ is $z_{i}$, which is defined as $z_{i}=V_{i}-V_{i+1}$ for $i=1,2,\dots,N$ and $V_{N+1}=V_{1}$. This new state variable better illustrates the occurrence of the different correlation state in the coupled systems. When the two neighbours, namely the $i^{th}$ and $(i+1)^{th}$ oscillators, oscillate coherently, the value of $z_{i}$ will be near zero. If the two neighboring oscillators, $i^{th}$ and $(i+1)^{th}$, evolve incoherently, $z_{i}$ will take a value much greater than zero. Here $\langle\cdot\rangle_{t}$ denotes the long time average, and $\langle z\rangle=\frac{1}{N}\sum\limits_{i=1}^{N}z_{i}$. Then the strength of incoherence is defined as

[TABLE]

where $\delta$ is the pre-defined threshold. In our work we choose $\delta=0.0025$, and $\Theta(\cdot)$ is the Heaviside step function. Here, $V_{\mbox{SI}}=1$ and $V_{\mbox{SI}}=0$ respectively represent the incoherent and coherent states, while the chimera states are identified by the values of $V_{\mbox{SI}}$ belonging in the interval $(0,1)$.

The variation of the maximum Lyapunov exponent $\Lambda_{\mbox{GS}}$ of the error systems (13) is shown in Fig. (4) with respect to $\epsilon$. For the coupling strength $\epsilon\geq 8.0$, the zero solution of the system (13) becomes stable, and accordingly, each node of the original system (5) is synchronized with its counterpart node in the auxiliary system (9). Here the coherent state is achieved together with the generalized synchronization state. To verify our analytical conditions, we plot $V_{\mbox{SI}}$ together with $\Lambda_{\mbox{GS}}$ by varying $\epsilon$. The top, middle, and bottom panels are, respectively, for PIs-, CIs-, and MIs-types initial conditions. The values of $V_{\mbox{SI}}$ are plotted by blue circle curves, while $\Lambda_{\mbox{GS}}$ are plotted by green square curves. For all three types of initial conditions, as the values of $V_{\mbox{SI}}$ abruptly drops from $1$, it is indicated for a slight increment of $\epsilon$ from zero that the incoherent states become unstable and chimera states appear. The values of $V_{\mbox{SI}}$ remain nonzero up to $\epsilon<8.0$; in that range of $\epsilon$ the values of $\Lambda_{\mbox{GS}}$ remains non-negative. From $\epsilon=8.0$, as soon as $V_{\mbox{SI}}$ becomes zero, $\Lambda_{\mbox{GS}}$ turns out to be negative. That is, for $\epsilon\geq 8$, the chimera states become unstable and coherent states become stable. For still higher values of $\epsilon$, $V_{\mbox{SI}}$ remains zero and $\Lambda_{\mbox{GS}}$ remains negative, which indicates the persistence of the coherent states. In conclusion, for all three types of initial conditions, our analytical study exactly matches the numerical simulations. This analysis provides a straightforward and precise characterization of the chimera to the coherent transition in coupled dynamical systems, probably irrespective of the complexity of the underlying network architecture and system dynamics.

To explore the complete scenario of the variation of $V_{\mbox{SI}}$ and $\Lambda_{\mbox{GS}}$, we compute the phase diagram by simultaneously varying $P$ and $\epsilon$ in the range of $P\in[1,99]\cap\mathbb{N}$ and $\epsilon\in[0,15]$. In the $(P,\epsilon)$ parameter space of a non-locally coupled Leech neuronal oscillator, the variations of $V_{\mbox{SI}}$ and $\Lambda_{\mbox{GS}}$ are drawn in Figs. 5(a) and 5(b) respectively. For this result, the PIs’ initial conditions are chosen. For the other two types of initial conditions (CIs and MIs), almost similar results are presented in the Appendix. From Fig. 5(a), it is noted that for lower values of $P$, the incoherent states persist up to $\epsilon\simeq 5.5$. Beyond that, chimera states emerge. But for higher values of $P$, chimera states appear as soon as we switch on the interaction. For the incoherent and chimera regions, the values of $\Lambda_{\mbox{GS}}$ remain non-negative in Fig. 5(b). For larger values of $\epsilon$, all nodes oscillate coherently, irrespective of the values of $P$. It emphasizes that the coupling strength plays an important role in changing the dynamical properties of the network. For higher values of $\epsilon$, as soon as chimera states are annihilated and coherent states emerge, the values of $\Lambda_{\mbox{GS}}$ become negative together with $V_{\mbox{SI}}=0$.

VI Global stability of the chimera and coherent states

From the discovery of the chimera states, the initial conditions play a crucial role. For many coupled systems, chimera states emerge only for a few specific initial conditions cmra_bs1 ; cmra_srep . Especially in multi-stable systems, the emergence of chimera states depends on the initial conditions. Since the isolated Leech neuronal model itself possesses multi-stable behavior, we now investigate the robustness of the initial conditions on the chimera states using basin stability measurement cmra_srep ; bs_delay . This is a non-linear, non-local approach, that can be easily applied to high-dimensional multi-stable systems. It is closely related to the volume of the basin of attraction. The basin stability tells us how stable the incoherent, chimera, and coherent states are against multifarious initial conditions. Recently, this measurement was used to quantify the multi-stable incoherent, chimera, and coherent states in non-locally coupled time-delayed Mackey-Glass oscillators cmra_srep . Here we have characterized different dynamical states by calculating the strength of incoherence. Now our main emphasis will be to identify the variation of basin stability in the parameter range of coupling strength $\epsilon$.

For numerical simulation, we choose random initial conditions from the phase-space volume $[-0.05,0.05]\times[0.0,0.5]\times[0.0,1.0]$ of the Leech oscillator. Then we integrate the entire system (5) for $Q$ (sufficiently large) different initial conditions, for which $T$ initial conditions finally arrive at our desired state. Then the basin stability $(B)$ of that particular state can be estimated as $\frac{T}{Q}$. The value $B\in[0,1]$, where $B=0$ indicates the instability of that state, and it is globally stable only if $B=1$. $B\in(0,1)$ indicates the coexistence of more than one dynamical state. Actually, basin stability corresponds to the probability of getting the state for any random initial conditions in the classical sense.

Now we investigate the probability of the emergence of the chimera state together with the incoherent and coherent states. The variations of $B$ for incoherent, chimera, and coherent states with respect to the coupling strength $\epsilon$ are plotted in Fig. (6). Here blue (black), red (dark gray), and green (gray) bars respectively represent the incoherent, chimera, and coherent states. For the values of $\epsilon$ almost near zero, the incoherent state solely dominates. By increasing the values of $\epsilon$, the $B$ of the chimera states abruptly increases, while the incoherent states are annihilated as well. From $\epsilon=2.0$, the chimera states become mono-stable, and they remain so up to $\epsilon=5.0$. Beyond that value of $\epsilon$, the basin stability of the chimera states gradually decreases and that of the coherent states improves. At $\epsilon=10$, the basin stability of the coherent states reaches $1$ and acquires full space in the basin volume, which is maintained for still higher values of $\epsilon$. Such changes of $B$ give an indication of the annihilation of incoherent and chimera states and the global stability of coherent states by varying the coupling strength.

VII Stability of the complete synchronization states: master stability function approach

Now, we investigate the persistence of the coherent states in terms of the auxiliary-system approach. For this purpose, we assume that the isolated evolution function $F({\bf x})$ satisfies the Lipschitz condition, i.e. there exists a positive constant $M$ such that

[TABLE]

The tenacity of the trivial equilibrium point of systems (10) implies the persistence of the coherent states. For this, we choose a Lyapunov function

[TABLE]

where $tr$ denotes the transpose of a matrix. Taking the derivative of $G(t)$ with respect to time $t$ along the trajectory of the systems (11), we get

[TABLE]

Using the Lipschitz condition (16), it is easy to obtain the relation $[{\bf x}-{\bf y}]^{tr}[F({\bf x})-F({\bf y})]\leq M[{\bf x}-{\bf y}]^{tr}[{\bf x}-{\bf y}]$. Then considering ${\bf Z}=[{\bf z}_{1},{\bf z}_{2},\cdots,{\bf z}_{N}]^{tr}$, we get

[TABLE]

where $\otimes$ is the Kronecker matrix product and $\lambda(I_{N}\otimes\Gamma)$ denotes the smallest eigenvalue of the matrix $I_{N}\otimes\Gamma$. Now for the two square matrices $A$ and $B$ of order $N_{1}$ and $N_{2}$ respectively, if $\{\lambda_{i}:1\leq i\leq N_{1}\}$ and $\{\mu_{j}:1\leq j\leq N_{2}\}$ are the respective sets of eigenvalues, then the eigenvalues of $A\otimes B$ are $\lambda_{i}\mu_{j}$ for $i=1,2,\dots,N_{1}$ and $j=1,2,\dots,N_{2}$. This yields $\lambda(I_{N}\otimes\Gamma)=\lambda(I_{N})\lambda(\Gamma)=\lambda(\Gamma)$. If all the eigenvalues of the inner coupling matrix $\Lambda$ are positive, then we get the persistence condition as

[TABLE]

By the Lasalle’s invariance principle LaSalle , the system (10) is asymptotically stable for these values of $\epsilon$, hence our original system oscillates coherently. The persistence condition (20) provides that once the coherent states are achieved, they will remain for any higher coupling strength if all the eigenvalues of the matrix $\Gamma$ are non-zero.

In our case, we observe that for higher values of $\epsilon$, qualitatively two different types of snapshots are exhibited: one is a smooth curve and the other is a straight line. The straight line snapshot indicates the complete synchronization among the coupled neurons. For this, we analytically derive the stability condition of the complete synchronization state, a special type of coherent state. For this type of coherent state, the snapshot of the oscillators exhibits a straight line. Each neuronal oscillator is identical and coupled by bidirectional electrical synapses. The individual evolution and interaction dynamics are both differentiable. Our coupling scheme is linear diffusive coupling; the complete synchronization manifold will be an invariant manifold. So the stability criterion of the complete synchronization states can be analyzed by the master stability function approach msf , which gives the local stability criterion of the complete synchronization state. When complete synchronization occurs, let each node of the network evolve with the synchronization solution ${\bf x}_{0}(t)$, which satisfies the evolution equation $\dot{\bf x}=F({\bf x})$ from Eq. (1). To gauge the local stability of the complete synchronization state, we perturb the $i^{th}$ node from the complete synchronization state by an amount $\delta{\bf x}_{i}(t)$. Then the current state of the $i^{th}$ node is ${\bf x}_{i}={\bf x}_{0}+\delta{\bf x}_{i}$. Considering small perturbations and expanding around the complete synchronization solution up to first order, we get the linearized equation of the error systems for the complete synchronization state as

[TABLE]

To determine the transverse error components, we project the error vector $(\delta{\bf x}_{1},\delta{\bf x}_{2},\dots,\delta{\bf x}_{N})$ to the space spanned by the Laplacian eigenvectors $U$. If $(\xi_{1},\xi_{2},\dots,\xi_{N})^{tr}$ is the projected error vector, then (\xi_{1},\xi_{2},\dots,\xi_{N})^{tr}=\big{(}U\otimes I_{d}\big{)}^{-1}(\delta{\bf x}_{1},\delta{\bf x}_{2},\dots,\delta{\bf x}_{N})^{tr}, which yields the projected error dynamics as

[TABLE]

The linearized equation for $i=1$, corresponds to the variational equation of the synchronization solution ${\bf x}_{0}$, which we have succeeded in separating from the other transverse directions. In our considered problem, the transverse error components can be written as

[TABLE]

Here $(V,m_{K2},h_{Na})$ satisfies the evolution equation of the synchronization solution (3). Next we investigate the effect of periodic initial conditions near and far away from the complete synchronization manifold by varying the coupling strength $\epsilon$.

Figure 7(a) shows the variation of $V_{\mbox{SI}}$ (blue circle line), $\Lambda_{\mbox{GS}}$ (green square line), and $\Lambda_{\mbox{CS}}$ (red diamond line) for the range of $\epsilon\in[0,40]$ with PIs near the complete synchronization manifold. Initially, for a very small value of $\epsilon$, the network shows the incoherent state. After increasing the value of $\epsilon$, the transition from chimera to coherent states emerges at $\epsilon=8.0$. From here, the value of $V_{\mbox{SI}}$ becomes zero, and it remains unchanged for still higher values of $\epsilon$. This is in good agreement with our stability analysis of the coherent state, as $\Lambda_{\mbox{GS}}<0$ for $\epsilon\geq 8.0$, while the values of $\Lambda_{\mbox{CS}}$ become negative at $\epsilon=25.0$. Hereafter, $\Lambda_{\mbox{CS}}$ maintains its negative values. So the complete synchronization state is locally stable for $\epsilon\geq 25.0$. To understand the qualitatively different coherent states, we plot the snapshots (middle panel) and corresponding spatio-temporal plots (lower panel) for two different values of $\epsilon$. At $\epsilon=15.0$, where $\Lambda_{\mbox{GS}}<0$ and $\Lambda_{\mbox{CS}}>0$, Fig. 7(b) shows a smooth coherent profile, but it is not a straight line parallel to the $i-$axis. Therefore, at this coupling strength, the coupled systems are in coherent states, but not a complete synchronous state. The inset figure shows the phase-space in the $(m_{K2_{i}},V_{i})$ plane. Here the network of coupled oscillators evolves with a quasi-periodic state, and Fig. 7(c) shows the corresponding spatio-temporal plot. At another value of $\epsilon=27.0$, where $\Lambda_{\mbox{GS}}<0$ and $\Lambda_{\mbox{CS}}<0$, the snapshot in Fig. 7(d) shows that the entire network is in a complete synchronization state. The complete synchronization solution is a period-1 limit cycle (inset figure), and the spatio-temporal plot of Fig. 7(e) depicts its amplitude variations. The variation of $\Lambda_{\mbox{CS}}$ in the $(P,\epsilon)$ plane for a complete synchronization state is shown in the Appendix. For higher values of $\epsilon$, the complete synchronization state persists when the initial conditions are taken near the complete synchronization manifold.

However, this complete synchronization state is not globally stable for all $\epsilon$ in Fig. 7(f) (left axis, blue circle curve). Particularly for the V-shaped PIs (far away from the complete synchronization manifold), the complete synchronization error

[TABLE]

first becomes zero at $\epsilon=25.0$. But it again shows a non-zero value for $\epsilon\in[29.0,33.0]$. In this range of $\epsilon$ with V-shaped initial conditions, $\Lambda_{\mbox{GS}}$ remains negative. Therefore, here the coherent state in the network is maintained. So here the complete synchronization state transforms to the generalized synchronization state. The complete synchronization state again returns at $\epsilon=33.0$ for these initial conditions.

Due to this multi-stable behavior of the complete synchronization state, we plot the basin stability $(B_{CS})$ of this state in Fig. 7(f) (right axis, green square curve). The values of $B_{CS}$ remains zero up to $\epsilon=22.67$. By increasing $\epsilon$, the values of $B_{CS}$ gradually increases up to $\epsilon=26.67$, where it reaches the value $0.97$. Upon further increasing the values of $\epsilon$, it decreases up to $\epsilon=30.67$, where $B_{CS}=0.7578$. For higher values of $\epsilon$, the $B_{CS}$ again increases and finally reaches its maximum value $1.0$ at $\epsilon=34.67$. At $\epsilon=30.0$, however, $\Lambda_{\mbox{CS}}$ is negative in Fig. 7(a) for small perturbations of their initial conditions. At $\epsilon=30.0$, the complete synchronization may not be stable (cf. Fig. 7(f)), while the generalized synchronization state is stable (cf. Fig. 6). Figures 7(g) and 7(h) show the snapshot and corresponding spatio-temporal plot for the generalized synchronization state. Again at $\epsilon=36.0$, the complete synchronization state is globally stable as $B_{CS}\backsimeq 1.0.$ Figure 7(i) shows the snapshot for the complete synchronization state, where each neuron oscillates with the periodic solutions (inset figure), and their amplitude variations are shown in Fig. 7(j).

VIII Conclusions

In this paper, a neuronal network composed of Leech neurons assumed to be coupled non-locally through gap junctions has been studied. Leech neuronal oscillator exhibits the coexistence of periodic and chaotic attractors. These two stable states possess a well-separated basin of attraction in phase space. For numerical simulations, three types of ‘V’ shaped initial conditions were chosen: a periodic basin, a chaotic basin, and a mixed basin. For a weaker synaptic strength, non-stationary chimera states emerge in this neuronal network. Three types of initial conditions possess different types of phase space of the entire network. Beyond a critical interaction strength, the chimera states become unstable, and coherent states appear, where the snapshot of the state variables follows a smooth profile. From this, we assert that in coherent states, the oscillators are in a generalized synchronization state. Using the auxiliary-system approach, we analytically derive the necessary condition for the chimera to coherent transitions, which are in excellent agreement with the value of the strength of incoherence. In a specified coupling range, the observed states are dependent on the initial conditions. By basin stability measurement, we have investigated the robustness of the chimera states with respect to the initial conditions. We also derived analytically the persistence conditions of the coherent states. This work contributes to a better understanding of the emergence of chimera states and their transition to coherent states in neuronal networks. Our analysis of the stability of coherent states will be useful in shedding light on the remarkable domain of chimera states.

Acknowledgements.

DG was supported by DST-SERB (Department of Science and Technology), Government of India (Project no. EMR/2016/001039). MP was supported by the Slovenian Research Agency (Grants J4-9302, J1-9112, and P1-0403).

IX Appendix

In this appendix, we study the transition from chimera to coherent states in the $(P,\epsilon)$ plane using measurements of the strength of incoherence $V_{\mbox{SI}}$ and the maximum Lyapunov exponent for generalized synchronization $\Lambda_{\mbox{GS}}$. For this, we choose the initial conditions from the chaotic and mixed basins of attraction. Then we plot the $(P,\epsilon)$ plane for the complete synchronization state using the master stability function $\Lambda_{\mbox{CS}}$.

In the $(P,\epsilon)$ parameter space, the variations of $V_{\mbox{SI}}$ and $\Lambda_{\mbox{GS}}$ are drawn in Figs. 8(a) and 8(b), respectively, by taking CIs initial conditions. For lower values of $P$, the incoherent states persist up to $\epsilon\simeq 5.5$, while for higher values of $P\geq 10$, they persist up to $\epsilon=0.25$. Beyond that, chimera states emerge, that persist up to a certain value of the synaptic strength, after which coherent states emerge. The chimera to coherent transition point behaves non-monotonically up to $P=60$. After that, this critical value decreases monotonically as $P$ increases. For the incoherent and chimera regions, the values of $\Lambda_{\mbox{GS}}$ remain non-negative in Fig. 5(b). For larger values of $\epsilon$, when coherent states appear the values of $\Lambda_{\mbox{GS}}$ becomes negative along with the $V_{\mbox{SI}}=0$.

Figures 9(a) and 9(b) respectively represent the values of $V_{\mbox{SI}}$ and $\Lambda_{\mbox{GS}}$ in the $(P,\epsilon)$ plane by taking MIs-type initial conditions. Both subplots are in excellent agreement with each other. For lower values of $P$, the incoherent states persist up to $\epsilon\simeq 5.5$, but for higher values of $P$, chimera states appear once we switch on the interaction. When Fig. 9(a) shows incoherent or chimera states, the value of $\Lambda_{\mbox{GS}}$ is positive there. For higher values of $\epsilon$, as soon as chimera states annihilate and coherent states emerge, the values of $\Lambda_{\mbox{GS}}$ become negative.

Figure 10 represents the parameter region in the $(P,\epsilon)$ plane for the complete synchronization states, by taking PIs initial conditions. The color bar shows the variation of the master stability function $\Lambda_{\mbox{CS}}$. The complete synchronization states are enhanced monotonically with respect to both the varying parameter $p$ and $\epsilon$.