Coupling conditions for globally stable and robust synchrony of chaotic systems

TL;DR

This paper introduces general coupling conditions for chaotic systems that ensure global stability and robustness of synchronization, applicable to various network motifs and models, including Lorenz and R"ossler systems.

Contribution

It proposes a novel set of coupling conditions involving cross-coupling links that enhance synchronization stability and robustness against parameter perturbations.

Findings

Coupling conditions ensure global stability of complete synchronization.

Cross-coupling links improve robustness against parameter perturbations.

Large networks can maintain synchrony despite node perturbations, transitioning to generalized synchronization.

Abstract

We propose a set of general coupling conditions to select a coupling profile (a set of coupling matrices) from the linear flow matrix (LFM) of dynamical systems for realizing global stability of complete synchronization (CS) in identical systems and robustness to parameter perturbation. The coupling matrices define the coupling links between any two oscillators in a network that consists of a conventional diffusive coupling link (self-coupling link) as well as a cross-coupling link. The addition of a selective cross-coupling link in particular plays constructive roles that ensure the global stability of synchrony and furthermore enables robustness of synchrony against small to non-small parameter perturbation. We elaborate the general conditions for the selection of coupling profiles for two coupled systems, three- and four-node network motifs analytically as well as numerically using benchmark models, the Lorenz system, the Hindmarsh-Rose neuron model, the Shimizu-Morioka laser model, the R\"ossler system and a Sprott system. The role of the cross-coupling link is, particularly, exemplified with an example of a larger network where it saves the network from a breakdown of synchrony against large parameter perturbation in any node. The perturbed node in the network transits from CS to generalized synchronization (GS) when all the other nodes remain in CS. The GS is manifested by an amplified response of the perturbed node in a coherent state.

Equations

$$\begin{array}[]{l}\dot{x}_{1,2}=\sigma({y}_{1,2}-{x}_{1,2})+\epsilon_{1}(x_{2,1}-x_{1,2})\\ \dot{y}_{1,2}=r_{1,2}x_{1,2}-y_{1,2}-x_{1,2}z_{1,2}\\ \dot{z}_{1,2}=-bz_{1,2}+x_{1,2}y_{1,2}\end{array}$$

$$\begin{array}[]{l}\dot{x}_{1}=\sigma({y}_{1}-{x}_{1})+\epsilon_{1}(x_{2}-x_{1}),\\ \dot{x}_{2}=\sigma({y}_{2}-{x}_{2})+\epsilon_{1}(x_{1}-x_{2})+\epsilon_{2}(y_{1}-y_{2})\\ \dot{y}_{1,2}=r_{1,2}x_{1,2}-y_{1,2}-x_{1,2}z_{1,2}\\ \dot{z}_{1,2}=-bz_{1,2}+x_{1,2}y_{1,2}\end{array}$$

$$\dot{e_x}=\sigma (e_y-e_x)-2{ϵ}_1{e}_x-{ϵ}_2{e}_y,$$

$$\dot{V}(e_x)=-(\sigma +2{ϵ}_1)e_x^2+(\sigma -{ϵ}_2)e_x{e}_y$$

$$\dot{V^{\prime }}(e_y,e_z)=-e_y^2-b{e}_z^2<0$$

$$\begin{array}[]{l}\dot{V}^{*}(e^{*})=-\sigma e_{x}^{2}-e_{y}^{*2}-be_{z}^{*2}<0.\end{array}$$

$$\begin{array}[]{c}\dot{\bf{x}}=g({\bf{x}})=F{\bf{x}}+L({\bf{x}})+{\bf C}\end{array}$$

$$F=\left(\begin{array}{ccc}-\sigma & \sigma & 0\\ r& -1& 0\\ 0& 0& -b\end{array}\right),L(\mathbf{x})=\left(\begin{array}{c}0\\ -xy\\ xz\end{array}\right),\mathbf{C}=0$$

$$F=\left(\begin{array}{ccc}0& 1& -1\\ 0& -1& 0\\ \mu r& 0& -\mu \end{array}\right);\left(\begin{array}{ccc}0& 1& 0\\ r& -b& 0\\ 0& 0& -a\end{array}\right);\left(\begin{array}{ccc}0& -1& -1\\ 1& a& 0\\ 0& 0& -c\end{array}\right);\left(\begin{array}{ccc}0& -a& 0\\ 1& 0& 1\\ r& 0& -1\end{array}\right).$$

$$H_s=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right);\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right);\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right);\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)$$

$$H_c=\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right);\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right);\left(\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ 0& 0& 0\end{array}\right);\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right)$$

$$\begin{array}[]{l}\dot{x}_{1}=ax_{1}^{2}-x_{1}^{3}+{y}_{1}-z_{1}+I+\epsilon_{1}(x_{2}-x_{1})+\epsilon_{2}(y_{2}-y_{1})\\ \dot{x}_{2}=ax_{2}^{2}-x_{2}^{3}+{y}_{2}-z_{2}+I+\epsilon_{1}(x_{1}-x_{2})\\ \dot{y}_{1,2}=b_{1,2}(1-dx_{1,2}^{2})-y_{1,2}\\ \dot{z}_{1,2}=\mu[r(x_{1,2}-c)-z_{1,2}]\end{array}$$

$$\begin{array}[]{l}\dot{e}_{x}=ae_{x}e_{p}-\frac{e_{x}^{3}}{4}-\frac{3}{4}e_{x}e_{p}^{2}+e_{y}-e_{z}-2\epsilon_{1}e_{x}-\epsilon_{2}e_{y},\\ \dot{e}_{y}=b_{1}(1-dx_{1}^{2})-b_{2}(1-dx_{2}^{2})-e_{y},\\ \dot{e}_{z}=\mu re_{x}-\mu e_{z},\end{array}$$

$$\dot{V^{\prime }}(e_x,e_z)=$$

$$-\frac{e_x^4}{4}+(1-{ϵ}_2)e_y{e}_x-\frac{e_z^2}{r}.$$

$$\frac{3}{4}{e}_p^2-a{e}_p+2{ϵ}_1=0$$

$${e_p}_{1,2}=\frac{2a\pm \sqrt{a^2-6{ϵ}_1}}{3}.$$

$$\dot{V}(e_x,e_y,e_z)=-\frac{e_x^4}{4}-e_y^2-\frac{e_z^2}{r}<0$$

$$\begin{array}[]{l}\dot{e}_{y}=b_{1}(1-dx_{1}^{2})-b_{2}(1-dx_{2}^{2})-e_{y}\\ \quad\;=(b_{1}-b_{2})(1-dx_{1}^{2})-e_{y}\;=\frac{b_{1}-b_{2}}{b_{1}}(\dot{y}_{1}+y_{1})-e_{y},\end{array}$$

$$\begin{array}[]{l}\dot{y}_{1}(1-\frac{b_{1}-b_{2}}{b_{1}})-\dot{y}_{2}=-y_{1}(1-\frac{b_{1}-b_{2}}{b_{1}})+y_{2}\\ \dot{y}_{1}\frac{b_{2}}{b_{1}}-\dot{y}_{2}=-(y_{1}\frac{b_{2}}{b_{1}}-y_{2}).\\ \end{array}$$

$${\dot{V}}^{\ast }(e_x,e_y^{\ast },e_z)=-\frac{e_x^4}{4}-e_y^{\ast 2}-\frac{e_z^2}{r}<0$$

$${\dot{\text{x}}}_i=g({\text{x}}_i)+{ϵ}_1\sum _{j=1}^N{A}_{ij}{H}_s({\text{x}}_j-{\text{x}}_i)+{ϵ}_2\sum _{j=1}^N{B}_{ij}{H}_c({\text{x}}_j-{\text{x}}_i),$$

$$A=\left(\begin{array}[]{rrr}0&1&0\\ 1&0&0\\ 1&1&0\\ \end{array}\right),B=\left(\begin{array}[]{rrr}0&0&0\\ 1&0&0\\ 1&0&0\\ \end{array}\right).$$

$$H_s=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right);H_c=\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right).$$

$$\begin{array}[]{l}\dot{x}_{i}=y_{i}+\epsilon_{1}\sum\limits_{j=1}^{N}A_{ij}(x_{j}-x_{i})+\epsilon_{2}\sum\limits_{j=1}^{N}B_{ij}(y_{j}-y_{i})\\ \dot{y}_{i}=r_{i}(x_{i}-x_{i}z_{i})-by_{i}\\ \dot{z}_{i}=-az_{i}+ax_{i}^{2}\end{array}$$

$$A=\left(\begin{array}{ccc}0& 1& 1\\ 1& 0& 1\\ 1& 1& 0\end{array}\right),B=\left(\begin{array}{ccc}0& 0& 0\\ 1& 0& 0\\ 1& 0& 0\end{array}\right).$$

$$H_s=\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right);H_c=\left(\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ 0& 0& 0\end{array}\right).$$

$$\begin{array}[]{l}\dot{x}_{i}=-{y}_{i}-{z}_{i}+\epsilon_{2}\sum\limits_{j=1}^{N}B_{ij}(z_{j}-z_{i})\\ \dot{y}_{i}=x_{i}+ay_{i}+\epsilon_{1}\sum\limits_{j=1}^{N}A_{ij}(y_{j}-y_{i})\\ \dot{z}_{i}=b_{i}+x_{i}z_{i}-cz_{i}\end{array}$$

$$A=\left(\begin{array}[]{cccc}0&1&1&0\\ 1&0&1&0\\ 1&1&0&0\\ 1&1&1&0\\ \end{array}\right),B=\left(\begin{array}[]{cccc}0&0&0&0\\ 1&0&0&0\\ 1&0&0&0\\ 1&0&0&0\\ \end{array}\right)$$

$$H_s=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),H_c=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right).$$

$$\begin{array}[]{l}\dot{x}_{i}=-ay_{i}+\epsilon_{1}\sum\limits_{j=1}^{N}A_{ij}(x_{j}-x_{i})\\ \dot{y}_{i}=x_{i}+z_{i}+\epsilon_{2}\sum\limits_{j=1}^{N}B_{ij}(z_{j}-z_{i})\\ \dot{z}_{i}=r_{i}(x_{i}+y_{i}^{2})-z_{i},\end{array}$$

Full text

††preprint: APS/123-QED

I Introduction

Synchrony is a most desirable state in many real-world dynamical networks such as the human brain Bablo ; Lytt , power grid Macho , ecological network Hirota , which must be stable against small to non-small perturbations. An instability due to a perturbation (intrinsic or extrinsic) may push a synchronous network to the point of failure of a desired performance (schizophrenia, power black-out, desertification etc.). To derive stability of such networks’ synchrony, the master stability function (MSF) is a breakthrough concept Pecora that determines the stability condition of synchronization in networks of identical oscillators. It encompasses the role of a network structure and the dynamics of the nodes as well, however, synchrony is locally stable since it is based on linear stability analysis. The MSF stability condition provided no clue how to prevent a breakdown of synchronization or failure of the networks’ desired function against a parameter perturbation. Recently, a basin stability approach Wiley ; Menck was undertaken to assess the stability region of the basin of attraction of a synchronous state in complex networks and to test how rewiring of coupling links can change the degree of stability to a smaller or a larger subbasin. A global stability can indeed enlarge the synchronous basin of a dynamical network beyond local stability and provide robustness against small or non-small perturbation and, thereby save a synchronous network from drifting away to a disaster. However, imposing global stability of synchronization in a dynamical network is still an enigma.

Realization of global stability of a synchronized state in a relatively small ensemble of oscillators is well known by a design of coupling Pad ; Chen ; Jiang that uses the Lyapunov function stability (LFS) approach. The main weakness of the design of coupling approach Pad ; Chen ; Jiang ; Grosu is the cost of coupling due to a relatively large number of coupling links necessary for global stability of synchronization and the essential presence of nonlinear coupling. In the case of a parameter mismatch, the number of coupling links is even larger. Earlier studies Motter ; Gu ; Hagberg ; Duan on dynamical networks have already shown that many coupling links does not mean better synchrony. By this time, it is also known that synchronization of dynamical networks can be enhanced by rewiring directed links Zeng1 or by deletion of redundant links Zeng2 . Playing Olusola with the number of coupling links for finding an appropriate choice skardal or redundancy of links in a network Schult is a judicious approach to improve the quality of synchronization. A natural question arises if addition or deletion of suitable coupling links can enlarge the basin of a synchronous state and thereby save an ensemble of oscillators from the brink of a breakdown? Is it possible to choose an appropriate set of coupling links for an ensemble of oscillators that may ensure globally stable synchrony? In a real world network, a parameter perturbation (system or coupling parameter) may originate due to internal/external disturbance when synchrony breaks down. Establising robustness of a desired synchronous state against such parameter perturbation is an important issue to resolve although not an easy task. We make an attempt here to address the issues by selecting an appropriate coupling profile consisting of linear coupling only for any two nodes in a network of oscillators.

From our past studies Suman of two coupled systems, we know that addition of cross-coupling links over and above the conventional diffusive coupling can realize global stability of synchrony and adds robustness to induced heterogeneity. However, the choice of coupling links was arbitrarily decided by the LFS condition and there were several available options. This leads us to a search for an appropriate coupling profile for dynamical systems, in general, with fewer diffusive scalar coupling and seletive cross-coupling that could suffice a global stability of synchrony in an ensemble of the dynamical system. Another important target is to ensure robustness of synchrony against heterogeneity due to a parameter mismatch ranging from a small to a large value.

In this paper, we propose a set of general coupling conditions that allows a selection of self-coupling and cross-coupling matrices, which is called as the coupling profile, from the LFM of dynamical systems and ensures global stability of a synchronous state in a network of the dynamical system. A coupling matrix determines the coupling functions that involve specific state variables of two dynamical systems to build up coupling links. Each coupling function is considered as a link between any two dynamical nodes. The conventional diffusive scalar coupling is redefined here as self-coupling because it involves two similar variables of the coupled systems and added to the dynamics of the same state variables. The cross-coupling also involves two similar variables Suman ; Mishra ; Motter1 of the systems, but adds to the dynamics of a different state variable. An appropriate selection of the cross-coupling link as well as the conventional self-coupling link realizes global stability of CS, expands the range of critical coupling for CS beyond the range defined by the MSF and thereby adds robustness to any drifting of the coupling paramater and, furthermore achieves robustness against perturbation in system parameters. As a result, if any one of the coupled systems is manually detuned by a small or a large value, CS simply transits to a type of GS Rulkov . The detuned system’s response is an amplified/attenuated replica of the unperturbed system. The coupled system maintains the dynamics of the isolated system, which we explain later in the text. The salient features of the selection of a coupling profile, especially, the addition of selective cross-coupling links, are exemplified when a larger network is considered. In a large network of dynamical systems, in a CS state, if a parameter of a node is perturbed that particular node transits to a GS state when all the other nodes remain in the CS state. The perturbed node’s attractor is an amplified replica of all the other nodes’ attractors.

We first presented an example of two-coupled Lorenz systems m-1 in section II to elaborate how addition of a cross-coupling link over and above the self-coupling realized global stability of synchrony and adds robustness to synchrony. Then we extended the results to frame a set of general conditions how a coupling profile can be appropriately selected from the LFM of a dynamical system and validate the scheme, in section III, analytically and numerically, using the paradigmatic model systems (Shimizu-Morioka (SM) laser model m-3 , Rössler system m-3 and a Sprott system m-4 ) for 2-coupled model systems. In section IV, we presented analytical details for a second example of a 2-coupled Hindmarsh-Rose (HR) model m-2 . We succesfully extended the results to 3-, 4-node network motifs Alon in section V. Analytical details for 2-coupled systems were presented in an Appendix. For analytical details of the network motifs, see the Supplementary Material (SM). Finally, an example of a larger network was presented in section VI with a summary in section VII.

II Coupled Lorenz systems

Two identical Lorenz oscillators are assumed synchronized under a bidirectional self-coupling,

[TABLE]

where subscripts denote oscillators (1, 2). Two systems are assumed identical, $r_{1}=r_{2}=r$, and a critical coupling $\epsilon_{1}=\epsilon_{MSF}>3.9225$ is determined using the MSF when the CS state, $x_{1}=x_{2}$, $y_{1}=y_{2}$, $z_{1}=z_{2}$, emerges as locally stable. If the coupling strength drifts to a value lower than the critical value, $\epsilon_{1}=3.8<\epsilon_{MSF}$, synchrony is lost (Fig. 1a). To mitigate this breakdown of synchrony, one directed cross-coupling link is added Pad to the coupled system in Eq. (1): a linear diffusive coupling function involving the $y_{1,2}$-variables is added to the dynamics of the $x_{2}$-variable of the second oscillator,

[TABLE]

where $\epsilon_{2}$ is the strength of the cross-coupling link. A directed cross-coupling link is basically added to the second oscillator when the first oscillator’s dynamics remains unchanged. The error dynamics $\dot{e}_{x}=(\dot{x}_{1}-\dot{x}_{2}$) changes,

[TABLE]

while $\dot{e}_{y}=\dot{y}_{1}-\dot{y}_{2}$, $\dot{e}_{z}=\dot{z}_{1}-\dot{z}_{2}$ remain unchanged. For a choice of Lyapunov function, $V(e_{x})=\frac{1}{2}e_{x}^{2}$,

[TABLE]

that establishes an asymptotic stability of $x_{1}=x_{2}$ if $\epsilon_{1}>-\sigma/2$ and $\epsilon_{2}=\sigma$. It leads to a stability relation of a Lyapunov function, $V^{\prime}(e_{y},e_{z})=1/2e_{y}^{2}+1/2e_{z}^{2}$,

[TABLE]

that ensures global stability of the CS manifold, $x_{1}=x_{2}$, $y_{1}=y_{2}$ and $z_{1}=z_{2}$, for identical systems. The addition of the directed cross-coupling link expands the range of critical self-coupling strength largely as defined by the condition, $\epsilon_{1}>-\sigma/2$, which is much lower than $\epsilon_{MSF}$. The global stability of CS thus sustains (Fig. 1b) a large drifting of the self-coupling strength below the critical value $\epsilon_{MSF}$. The synchronous dynamics maintains the isolated dynamics of the coupled systems since both the self- and cross-coupling functions are effectively vanished in the CS state (see Eq.(2)).

Furthermore, the cross-coupling link plays the most important role of realizing robustness of synchrony against a system parameter perturbation. For a demonstration, a parameter $r_{2}=(r_{1}+\Delta r$) is detuned manually in Eq.(2) when the error dynamics is modified (details are given in refs. equ1 ; equ2 ) $\dot{e}^{*}=\left[-(\sigma+2\epsilon_{1})e_{x}+(\sigma-\epsilon_{2})e_{y},\;-e_{y}^{*}-x_{1}e_{z}^{*},\;x_{1}e_{y}^{*}-be_{z}^{*}\right]^{T}$ where $e_{y}^{*}=\frac{r_{2}}{r_{1}}y_{1}-y_{2}$ and $e_{z}^{*}=\frac{r_{2}}{r_{1}}z_{1}-z_{2}$ and this leads to a revised stability condition of a modified Lyapunov function ($V^{*}$),

[TABLE]

provided $\epsilon_{1}>-\sigma/2$ and $\epsilon_{2}=\sigma$. A new globally stable synchronous state, $x_{2}=x_{1}$, $y_{2}=\frac{r_{2}}{r_{1}}y_{1}$ and $z_{2}=\frac{r_{2}}{r_{1}}z_{1}$ emerges which we call as a type of GS state. A linear functional relation emerges between the coupled systems. As an example, for a detuning of $r_{2}$ by twice ($\frac{r_{2}}{r_{1}}$=2), the new emergent coherent state, $x_{2}=x_{1}$, $y_{2}=2y_{1}$ and $z_{2}=2z_{1}$, is globally stable. Effectively, the detuned system’s attractor (red) is amplified along the $y$- and $z$-direction exactly twice that of the unperturbed system’s attractor (black) as shown in Fig. 1(c).

CS simply transits to GS due to the parameter detuning without a loss of overall coherence or synchrony in the coupled system. This is manifested by a rotation of the CS manifold ($y_{1}=y_{2}$ in black line) in Fig. 1(d) to GS manifold ($y_{2}=2y_{1}$, red line). In the GS state, $x_{1}=x_{2}$ condition is still preserved when the self-coupling function vanishes. The first oscillator is thus not affected by the cross-coupling and hence maintains the isolated dynamics and, plays an effective role of a driver. The second oscillator also continues with the original dynamics because it simply works as a slave to the first oscillator in a coherent state, but with an amplified replica (red) as shown in Fig. 1(c). If $r_{1}$ is detuned, instead, the qualitative behavior remains same except that the coupled dynamics will now change since the first oscillator has a change of dynamics. The final effect is a dramatic improvement in the quality of synchronization that is simply implemented by a selection of two linear coupling links, especially, the addition of one selective cross-coupling link.

III Linear Flow Matrix: Coupling profile

The particular selection of a coupling profile (self-coupling and cross-coupling links) for the Lorenz system, can really be made in a systematic manner from the LFM of the system and, we can frame a set of general coupling conditions for the selection of a coupling profile for many dynamical systems,

[TABLE]

$g$:$\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is the flow of the system, $F$ is the LFM (${n\times n}$ matrix, $n$=3 for our example systems), $\textbf{x}=[x,y,z]^{T}$, $T$ denotes a transpose, and $L({\bf x})$ represents the nonlinear functions and C is a constant matrix. For the Lorenz system,

[TABLE]

By an inspection of $F$ of the Lorenz system, we suggest that no self-coupling is necessary, since all the diagonal elements are negative and hence all the error functions due to the self-coupling functions appear as negative definite in $\dot{V}(e)$ of Eq. (4). This is supported by the condition of a critical self-coupling $\epsilon_{1c}>-\sigma/2$ in (6) that includes a condition, $\epsilon_{1}=\epsilon_{1c}=0$, as described above. However, we always apply one scalar self-coupling function to Lorenz systems so as to realize first a locally stable CS state (MSF based stability). On the other hand, a nonzero element exists in the upper triangle of $F$ that is connected to the $\dot{x}_{1,2}$ dynamics. This suggests an addition of one selective directed cross-coupling link defined by a coupling function involving the $y_{1,2}$ variable to the $x_{2}$ dynamics as shown in Eq. (2) since the nonzero element is connected to the $y_{1,2}$ variable in $\dot{x}_{2}$ and it eliminates the additional error functions that appear in $\dot{V}(e)$ due to the presence of $y_{1,2}$ variable in the dynamical equations of the system. As example, $e_{y}$ is removed from Eq. (3) by the addition of a cross-coupling function (representing a coupling link) using the $y_{1,2}$-variable, and it suffices satisfying the condition, $\dot{V}(e)<0$. By a selection of two linear self- and cross-coupling links between two Lorenz systems, the global stability of synchrony and other benefits such as the robustness of synchrony to large parameter perturbation, as described above, are ensured.

This observation from Lorenz systems leads us to define a set of general conditions regarding the selection of a coupling profile, between any two nodes of a network, from the LFM, $F$={$f_{ij}$}, of any dynamical system that represents each node of the network. We always apply a bidirectional self-coupling link first to realize a locally stable (MSF condition) synchronization. The question is which self-coupling link is most appropriate? Then what is the appropriate choice of the cross-coupling link? We frame a general guideline that the diagonal elements of the LFM determine the choice of self-coupling links and the off-diagonal elements, in the upper triangle, determine the choice of cross-coupling links. Accordingly, we propose the coupling conditions, (1) a primary choice of bidirectional self-coupling link is made for each positive diagonal element ($f_{ii}>0$) in the LFM; for any negative diagonal element, self-coupling is redundant, (2) if one or more zeros appear in the diagonal elements, at least one self-coupling link is added to the related dynamical equation, (3) a cross-coupling link is essential for each nonzero off-diagonal element ; however, no cross-coupling is needed if $f_{ij}$ =-$kf_{ji}$ where $k$ is a constant.

Above proposition clearly justifies our selection of self-coupling and directed cross-coupling links in two Lorenz systems as demonstrated above. We present more supporting examples with the paradigmatic HR neuron model m-2 , the SM laser model m-3 , the Rössler system m-4 and a Sprott system m-5 whose LFMs are ordered left to right, respectively,

[TABLE]

For coupling two systems, we now proprose coupling profiles for above four example models. For a HR system, the diagonal element $f_{11}$ of the LFM is 0 and hence we apply one bidirectional self-coupling link, involving the $x_{1,2}$-variables, to the dynamical equations of $x_{1}$; this is an appropriate selection of the self-coupling link for two HR systems since all other diagonal elements ($f_{22}$ and $f_{33}$) are negative. In the upper triangle of the LFM, $f_{12}$ is 1 ($f_{12}\neq-f_{21}=0$) and hence one directed cross-coupling link involving the $y_{1,2}$-variables is to be added to the dynamics of $x_{1}$. No additional cross-coupling is needed since $f_{31}=-kf_{13}~{}(k=\mu r)$. From the LFM of the SM laser system, second in the row, $f_{11}=0$, $f_{22}$ and $f_{33}$ are found negative, and the selection of self-coupling involves the $x_{1,2}$-dynamics and the cross-coupling ($f_{12}=1,f_{21}=r$) involves $y_{1,2}$ variables and both are added to $x_{1}$-dynamics. No other cross-coupling is needed since other off-diagonal elements are $f_{13}=f_{23}=0$. For the Rössler system $f_{11}=0$ and $f_{22}=a$ as seen from the LFM, third in the row, and hence we add one self-coupling involving the $y_{1,2}$-variables to the $y_{1,2}$-dynamics and one cross-coupling involving the $z_{1,2}$-variable to the dynamics of $x_{1}$- variable ($f_{13}=-1;f_{31}=0$); no other cross-coupling is needed since $f_{12}$=-$f_{21}$. The LFM of the Sprott system is shown last in the row. We have two options for bidirectional self-coupling ($f_{11}=f_{22}=0$, $f_{33}=-1$); one can select either of the pairs of state variables $x_{1,2}$ or $y_{1,2}$ witout any effect on the result. One cross-coupling involving the $z_{1,2}$ variable is added to the dynamics of $y_{2}$ since $f_{23}$=1 ($f_{32}=0$). No other cross-coupling is needed since $f_{12}=-af_{21}$.

As a result, guided by the above conditions (1)-(3) and as elaborated above, we obtain the following self-coupling and cross-coupling matrices, $H_{s}$ and $H_{c}$, respectively, for building coupled systems using the HR model, the SM model, the Rössler system and Sprott system, and presented in order from left to right, respectively,

[TABLE]

and

[TABLE]

The elements of the coupling matrices $H_{s}$ and $H_{c}$, in Eqs. (9) and (10), respectively, are chosen accordingly,

(1) HR system: one self-coupling link and one cross-coupling link, all elements of the coupling matrices are 0 except ${H_{s}}_{11}$=1 and ${H_{c}}_{12}$=1,

(2) SM system: one self-coupling link and one cross-coupling link, all elements of the coupling matrices are 0 except ${H_{s}}_{11}$=1 and ${H_{c}}_{12}$=1,

(3) Rössler system: one self-coupling link and one cross-coupling link, all elements of the coupling matrices are 0 except ${H_{s}}_{22}$=1 and ${H_{c}}_{13}$=1,

(4) Sprott system: one self-coupling link and one cross-coupling link. All elements of the coupling matrices are 0 except ${H_{s}}_{11}$=1 and ${H_{c}}_{23}$=1.

The coupling profiles are pictorially demonstrated in Fig. 2(a)-(d) for the above systems that always consist of one self-coupling (black arrow) and one cross-coupling link (red arrow) for 2-node coupled systems. The connectivity matrices (self- and cross-coupling links) are only different for different systems.

IV Two coupled systems

The proposed selection of coupling profile favors global stability of synchrony in two-coupled systems. Numerical results are summarized for all the systems in the immediate right panels of Figs. 2(e)-(h). Analytical details are provided for the HR model in this section and for rest of the models see the Appendix. For each coupled system, the coupling profile is selected, as suggested above, to realize a globally stable CS for identical systems. CS manifolds of identical oscillators are shown in black lines and the GS manifolds of the detuned systems are shown in red lines. CS to GS transition due to a parameter perturbation is reflected by a rotation of the synchronization manifold for a similar reason as explained above for the Lorenz system. The amplified response of the detuned systems’ attractor Suman under a perturbation is explained further using the example of the slow-fast HR systems.

The coupling profile ($H_{s}$ and $H_{c}$) is taken from the Eqs. (9)-(10) for the HR neuron model that suggests one bidirectional self-coupling and one directed cross-coupling for globally stable synchrony. Accordingly, two coupled HR systems appears as,

[TABLE]

The bidirectional self-coupling link $(x_{2,1}-x_{1,2})$ is added to the dynamics of $x_{1,2}$-variables while the cross-coupling link $(y_{2}-y_{1})$ is directed and added only to the dynamics of $x_{1}$ of the first oscillator. This choice can be reversed, i.e., ($y_{1}-y_{2})$ could be added alternatively to the dynamics of $x_{2}$ without any effect on the final results.

The evolution of the error functions $\bf e$=$[e_{x},e_{y},e_{z}]^{T}=[x_{1}-x_{2},y_{1}-y_{2},z_{1}-z_{2}]^{T}$ is given by,

[TABLE]

where, $e_{p}=x_{1}+x_{2}$ so that $x_{1}^{2}-x_{2}^{2}=e_{x}e_{p}$ and $x_{1}^{3}-x_{2}^{3}=\frac{e_{x}}{4}(e_{x}^{2}+3e_{p}^{2})$. For a global stability of $(e_{x}=0,e_{y}=0,e_{z}=0)$, we consider a Lyapunov function, $V(e)=\frac{1}{2}e_{x}^{2}+\frac{1}{2}e_{y}^{2}+\frac{1}{2\mu r}e_{z}^{2}$. We first check the stability of $x_{1}=x_{2}$ and $z_{1}=z_{2}$, separately, by defining a Lyapunov function, $V^{\prime}(e_{x},e_{z})=\frac{1}{2}e_{x}^{2}+\frac{1}{2\mu r}e_{z}^{2}$ when its time derivative is

[TABLE]

$\dot{V}^{\prime}$($e_{x},e_{z})<0$ provided $\epsilon_{2}=1$ and $\frac{3}{4}e_{p}^{2}-ae_{p}+2\epsilon_{1}\geq 0$. To satisfy the condition, we derive the roots of the equation,

[TABLE]

which are given by

[TABLE]

${e_{p}}_{1,2}$ will now be positive real if $a^{2}-6\epsilon_{1}\geq 0$ that implies ${\epsilon_{1}}\leq\frac{a^{2}}{6}$. Thus $\dot{V^{\prime}}(e_{x},e_{z})<0$ is satisfied when ${\epsilon_{2}}=1$ and ${\epsilon_{1}}\leq a^{2}/6$ and this implies $x_{1}=x_{2}$ and $z_{1}=z_{2}$ is asymptotically stable as $t\rightarrow\infty$. Substituting the condition in Eq. (12) and assuming identical systems ($b_{1}=b_{2}$), the error dynamics $\dot{e}_{y}=-e_{y}$ is found when,

[TABLE]

for ${\epsilon_{1}}\leq\frac{a^{2}}{6}$ and ${\epsilon_{2}}=1$. The synchronous state $x_{1}=x_{2},y_{1}=y_{2}$ and $z_{1}=z_{2}$, basically a CS state of the coupled identical HR system, is now globally stable.

Now the effect of heterogeneity on the stability of CS is tested by detuning the parameter, ($b_{1}\neq b_{2}$), other parameters are kept unchanged. The induced heterogeneity has no effect on Eq.(13) and related conditions and hence the stability of $x_{1}=x_{2}$ and $z_{1}=z_{2}$ is still preserved. The stability of $e_{y}$ is only to be checked against a parameter detuning. The equation of $\dot{e}_{y}$ after a perturbation is written from Eq.(12),

[TABLE]

and this leads to

[TABLE]

From (18), the revised error dynamics is $\dot{e}_{y}^{*}=-e_{y}^{*}$, where $e_{y}^{*}=y_{1}\frac{b_{2}}{b_{1}}-y_{2}$ is the modified error function. Accordingly, the Lyapunov function is redefined in terms of the modified error functions whose time derivative is

[TABLE]

and valid for same conditions, ${\epsilon_{1}}\leq\frac{a^{2}}{6}$ and ${\epsilon_{2}}=1$ as given above. It confirms an emergence of a new globally stable synchronous state, $x_{1}=x_{2}$, $y_{1}=\frac{b_{1}}{b_{2}}y_{2}$, $z_{1}=z_{2}$, due to parameter detuning, which we call as a GS state. This is manifested by a rotation of the CS hyperplane ($y_{1}=y_{2}$, black line) to a GS hyperplane ($y_{1}=\frac{b_{1}}{b_{2}}y_{2}$, red line) as shown in Fig. 2(e) where the angle of rotation is decided by the amount of detuning. For any choice of $b_{1}$ and $b_{2}$, the attractor of the detuned system shall be amplified or attenuated along the $y$-direction as determined by the amount of heterogeneity, i.e., by the ratio of (${b_{1}}/{b_{2}}$) as shown in Fig. 3(a). The $y_{2}$-variable (red line) is only amplified twice ($y_{2}=2y_{1}$; $y_{1}$ in black line) in Fig. 3(b) for a choice of $b_{1}=0.5$ and $b_{2}=1.0$. A complete in-phase synchrony is maintained between $y_{1}$ and $y_{2}$.

Another important point is that synchrony is lost under a purely self-coupling ($\epsilon_{2}=0$) if its strength drifts below $\epsilon_{1}<e_{MSF}$, however, restored by the addition of one directed cross-coupling link as shown in $x_{1}$ vs. $x_{2}$ plots in Figs. 3(c)-3(d) since the range of $\epsilon_{1}$ ($\leq\frac{a^{2}}{6}$) is expanded.

V Network Motifs

We focus here on network motifs which are building blocks of many real world networks Alon and check if the proposed general conditions for the selection of coupling profiles work in favor of global stability of synchrony. A few examples are presented that support our proposition and many more are cited in the Supplementary material SM . The dynamics of the $i$-th oscillator in a network of N-oscillators under self-coupling as well as cross-coupling is,

[TABLE]

$g({\bf x_{i}})$ denotes the flow of the $i^{th}$ node in isolation, $i\in\{1,...,N\}$, ${\bf x_{i}}\in\mathbb{R}^{n}$ is the state vector. $\epsilon_{1}$ and $\epsilon_{2}$ are the strength of self- and cross-coupling links, respectively, between any two nodes. $A=\{A_{ij}\}\in\mathbb{R}^{N\times N}$ is the adjacency matrix that defines the network topology via self-coupling links and $B=\{B_{ij}\}\in\mathbb{R}^{N\times N}$ is the connectivity matrix of the cross-coupling links of the network; $A_{ij}$=1 and $B_{ij}$=1, if $i^{th}$ node is connected to the $j^{th}$ ($j\neq i$) and 0 otherwise. $H_{s}$ and $H_{c}$ a usual define the self- and cross-coupling matrices for any pair of nodes.

We make appropriate selection of $H_{s}$ and $H_{c}$ between any two nodes of a network as guided by the LFM of dynamical systems representing each node. We are allowed with $N-1$ directed cross-coupling links from any arbitrarily chosen node to all other nodes that defines the connectivity matrix $B$ of cross-coupling links. A restriction is imposed that the addition of cross-coupling links must not make an original directed link undirected. In other words, $B$ must not change the original topology of a network-motif defined by $A$. The single node from where the outgoing directed cross-coupling links are connected to all other nodes, maintains the isolated dynamics at a globally stable CS state when all coupling functions vanish for identical nodes. This arbitrarily chosen node effectively plays the role of a driving oscillator and, obviously, maintains the isolated dynamics. In the case of a parameter perturbation, the network emerges into a new coherent state of GS type, as explained above, however, continues with the original dynamics so long as the newly emerged driver node’s parameter is not perturbed as discussed above for two Lorenz systems.

Using a few examples of network-motifs, the general applicability of our coupling profile scheme is establised analytically as well as numerically. All the network-motifs are drawn with self-coupling links in black arrows and cross-coupling links in dashed arrows (red). Results are summarized here; see the Supplementary material SM for analytical details and more examples.

V.1 Network motif: Shimizu-Morioka laser model

A particular 3-node network motif is first shown in Fig. 4(a) whose each node is represented by the SM laser model. The network topology $A$ is denoted by the self-coupling links in black arrows and the connectivity matrix of cross-coupling links in dashed arrows (red) is represented by $B$,

[TABLE]

Notice that the cross-coupling links do not change the network topology. $H_{s}$ and $H_{c}$ are obtained from Eqs.(9)-(10), for the SM model representing each node,

[TABLE]

This particular $H_{s}$ suggests that the self-coupling link must involve the $x_{1,2}$-variable and the cross-coupling link must involve the $y_{1,2}$-variables and both be added to the dynamics of $x_{1}$ as suggested by $H_{c}$. Two ($N-1$) directed cross-coupling links (dashed arrows) are added from node 1 to nodes 2 and 3. Replacing $H_{s}$, $H_{c}$ and using the SM model in (20), the governing equations of the $i^{th}$ node is,

[TABLE]

where $i=1,2,3$ and $a=0.375,b=0.826$ and $r_{i}$ is the parameter used for detuning. The global stability of CS in the motif is establised in ref.SM . Now consider $r_{1}=1$ for the oscillator node 1 in Fig. 4 and other two nodes 2 and 3 are perturbed so that $r_{2}=2,r_{3}=4$, respectively. The resultant effect is exactly same as seen for two node coupled system. Attractors of the detuned (positive detuning) nodes are amplified which is reflected in the separation of three synchronization hyperplanes in Fig. 4(b): the attractor of the first oscillator is plotted in a $y_{1}$ vs. $z_{1}$ plane and compared with the perturbed oscillators (2, 3) by taking a 3D plot using the $y_{2,3}$ variables along the $z$-axis. Two of them (blue and red) are amplified with respect to the reference hyperplane (black) and their size scaling depends upon the ratios, $\frac{r_{2}}{r_{1}}$=2 and $\frac{r_{3}}{r_{1}}$=4. The unperturbed oscillator works as a reference here which maintains the original dynamics and accordingly, other two emerge with a GS relation with the first. All three oscillators remain coherent, but the pertubed nodes show amplified versions of the original dynamics. The unperturbed oscillator plays a leadership role which we prove analytically (see Supplementary material SM . For a better clarity of the results, the synchronization manifolds $y_{1}$ vs. $y_{2,3}$ are plotted and the attractors of the 3-nodes are plotted in Fig. 5(a) and 5(b), respectively. Figure 5(a) shows that the synchronization manifolds of all three nodes collapse on the attractor in black when they are identical. In the case two detuned nodes (blue and red), their synchronization manifolds (blue and red lines) are rotated away from the reference (black) where the angle of rotation is decided by $\frac{r_{2}}{r_{1}}$=2 and $\frac{r_{3}}{r_{1}}$=4, respectively. Figure 5(b) shows 2D attractors (black, blue and red) where detuned nodes’ attractors (blue, red) are expanded along the $y$-axis. No other variable is amplified since $x_{1}=x_{2}=x_{3}$ and $z_{1}=z_{2}=z_{3}$ is maintained. The detuned nodes’ attractors expand along the $y$-axis only.

V.2 Network motif: Ring of 3-Rössler Systems

A second example of a network motif is shown in Fig. 6(a), a ring of three Rössler oscillators, basically, a globally coupled network. The adjacency matrix $A$ for self-coupling and the connectivity matrix of cross-coupling $B$ are

[TABLE]

The coupling profile $H_{s}$ and $H_{c}$ from (9) and (10) for the Rössler system is,

[TABLE]

The connectivity matrix of the cross-coupling $B$ involves 2- directed links only as shown in Fig. 6(a).

The dynamics of the network-motif under the chosen coupling profile is obtained from (20),

[TABLE]

$N$=3, $i\in\{1,2,3\}$. $\textbf{x}_{i}=(x_{i},y_{i},z_{i})^{T}$ is the state vector of the $i$-th node. The system parameters, $a=0.2,c=4.8$ and for identical oscillators $b_{i}=0.2$. The coupling profile realizes global stability of CS in the network-motif of identical nodes. The heterogeneity is introduced by changing parameters ($b_{1}<b_{2}<b_{3}$) in two nodes of the network.

All the nodes transits to a coherent GS state and globally stable for induced heterogeneity. In identical case ($b_{1}=b_{2}=b_{3}$=0.2), their dynamics will collapse to the CS hyperplane (black) as shown in Fig. 6(b). For positive detuning of two nodes ($b_{2}=0.4,b_{3}=0.8$), their attractors are amplified by scaling factors $\frac{b_{2}}{b_{1}}$=2 and $\frac{b_{3}}{b_{1}}$=4 and their synchronization hyperplanes (blue, red) are simply rotated away along the $z_{2}$- and the $z_{3}$-direction, respectively, from the CS hyperplane as depicted by their transverse direction (blue, red arrows). The attractor of the unperturbed node (black) act as the driver node ($b_{1}$=0.2). The amplification of the attractors is clarified by the time series plot of $z_{1,2,3}$ of three oscillators in Fig. 7(a), where two of them are the amplified replica of the third one. Other variables of the detuned nodes maintain CS relation with the unperturbed node. A rotational transformation of the synchronization hyperplanes of the detuned nodes compared to the reference CS hyperplane in 7(b) confirms their GS relation. Note that this motif is a symmetric one due to its global structure, $B$ can be chosen arbitrarily from any one of the nodes.

V.3 Network motif: 4-node Sprott system

Finally, a 4-node network-motif of the Sprott system m-5 is considered as shown in Fig. 8(a). $A$ and $B$ of the 4-node network-motif are,

[TABLE]

The coupling profile is,

[TABLE]

Substituting them to (20), the coupled dynamics of the network-motif of Sprott systems,

[TABLE]

where $a$=0.23 and $r_{i}$ is the heterogeneity parameter.

As shown in Fig. 8(b), the projected synchronous hyperplanes for three perturbed nodes are rotated away (blue, green and red, respectively, along the $z_{2}$-, $z_{3}$- and $z_{4}$-axis) from the reference hyperplane (black) and their respective scaling factors are $\frac{r_{i}}{r_{1}}$, $i=2,3,4$.

Further numerical confirmation is given in Fig. 9(a) where it is clearly seen that the time series of $z_{1,2,3,4}$-variables are affected by the induced heterogeneity of three nodes. Actually, three of them are amplified but remain coherent with all each other. This is confirmed by the modified error plot $e^{*}$ in Fig. (Fig. 9(b)) that converges to zero in time when $z_{1}=\frac{z_{2}}{2}=\frac{z_{3}}{3}=\frac{z_{4}}{4}$. Other variables remain unchanged and maintain CS relation in all three nodes.

VI Constructive role of cross-coupling in a network

It is demonstrated here how a single directed cross-coupling link makes a dramatic change in the quality of synchronization in a larger network shown in Fig. 10. The self-coupling links define the network topology (adjacency matrix $A$) in solid lines whose each node is represented by a Rössler system. The self-coupling matrix between any two nodes is chosen from the LFM ($H_{s}$ in Eq. (9) for Rössler system) to realize CS for all identical nodes (strength of each self-coupling link, $\epsilon_{1}$=0.21, $\epsilon_{2}$=0 in Eq.(24)). When node-6 (blue circle) is perturbed, synchrony is lost as shown in Fig. 11(a). One directed cross-coupling link ($\epsilon_{2}=-1$) as defined by $H_{c}$ in Eqn.(10) for a Rossler system is now added from node-12 to node-6 (dashed red arrow). Numerically solving Eq. (24) for N=16, we find that synchrony is restored in the network as shown in Fig. 11(a) after an instant $t=450$ when the cross-coupling link is added. Surprisinly, all the nodes return to CS and follows the isolated dynamics except the node-6 whose attractor is simply amplified/attenuated for a positive/negative detuning as shown in Fig. 11(b). Effectively, node-12 is now playing the role of a driver to node-6 as discussed above for network-motifs. In fact, one cross-coupling link from any one of the neighboring nodes (nodes, 2, 3, 7, 10, 12) to the perturbed node-6 will produce the same result. Synchronization in the network is monitored by defining the error $e_{r}$ as $e_{r}(t)=\sqrt{\sum_{i=1}^{N}\left[(\tilde{x}-x_{i})^{2}+(\tilde{y}-y_{i})^{2}\right]}$, where $N$ is the network size, $\tilde{x}=\frac{1}{N}\sum_{i=1}^{N}x_{i}$ and $\tilde{y}=\frac{1}{N}\sum_{i=1}^{N}y_{i}$. $e_{r}\neq 0$ means no synchronization.

VII Conclusion

We proposed a set of general coupling conditions for a selection of a coupling profile from the LFM of a dynamical system that realized globally stable synchrony in an ensemble of the system and a robustness of synchrony to perturbation in system parameter as well as the coupling parameter. The coupling profile consisted of self-coupling and cross-coupling matrices that defined the coupling links between any two nodes of a network. We explained how the selection of the self-coupling and the cross-coupling matrices were systematically made from the nature of the diagonal and off-diagonal elements, respectively, of the LFM of a dynamical system. Addition of selective cross-coupling links over and above the conventional diffusive coupling (defined here as self-coupling) in particular made a dramatic improvement in the quality of synchronization. Besides establising global stability of CS in the coupled system, it expands the range of critical coupling for CS from the range defined by the conventional MSF conditions and as a result, it saves synchrony from a breakdown in a situation of a drifting of the coupling parameter. Most importantly, the addition of the cross-coupling link added robustness to synchrony against large mismatch of system parameters. We demonstrated the constructive roles of our proposed coupling scheme, especially, the cross-coupling link using many example systems. We exemplified the constuctive role using an example of a larger network in a CS state by a manual detuning of a system parameter when all the nodes of the network remained unpertubed in the CS state, but the perturbd node showed a transition from CS to GS. The manifestation of the GS state is an amplified (attenuated) response of the perturbed node‘s attractor for a positive (negative) detuning of the system parameter. We validated our proposed coupling profile scheme with examples of 2- coupled systems, 3-node and 4-node network motifs using the Lorenz system, the Rössler system, the Hindmarsh Rose model, the Shimizu-Morioka laser model and a Sprott system as dynamical nodes. We supported our results by analytical as well as numerical verifications in all the example systems except the large network. The coupling scheme worked efficiently for many other network-motifs (see Supplementary Material). We plan to work for more rigiorus analysis of our results for large networks.

Acknowledgements.

P.K.R. acknowledge support by the CSIR (India) Emeritus scientist scheme. A.M. is supported by the UGC (India). S.K.D. is supported by the University Grants Commission (India) Emeritus Fellowsip. S.K.B. and S.S. acknowledge support by the DST-SERB (India) File No.-YSS/2014/000687.

Appendix: Two coupled systems

Shimizu-Morioka model: From (9) and (10), bidirectional self-couplings are introduced through the $x_{1,2}$ variables and one cross-coupling link involves $y_{1,2}$ variables and added to the dynamics of $x_{1}$ so that the governing equations are,

[TABLE]

Error dynamics are

[TABLE]

We first examine the synchronization of $x_{1}=x_{2}$ by considering a Lyapunov function $V(e_{x})=1/2e_{x}^{2}$. The time derivative of $V(e_{x}$) is,

[TABLE]

[TABLE]

if ${\epsilon_{2}}$=1 and $\epsilon_{1}>0;$ when $\lim_{t\rightarrow\infty}{e}_{x}(t)\rightarrow 0,$ that implies $x_{1}=x_{2}$ is asymptotically stable. Using this condition, we check the stability of $z_{1}=z_{2}$ when we find,

[TABLE]

$z_{1}=z_{2}$ is now asymptotically stable as $t\rightarrow\infty$. This implies $\dot{e}_{y}$=-$e_{y}$ for identical oscillators, $r_{1}=r_{2}$. The global stability of CS ($x_{1}=x_{2}$, $y_{1}=y_{2}$ and $z_{1}=z_{2}$) is thus established.

For a parameter detuning of $r_{1}$ of oscillator-1 from its identical values ($r_{1}\neq r_{2}$), the error function $e_{y}$ conly changes,

$\dot{e_{y}}=(r_{1}-r_{2})(x_{1}-x_{1}z_{1})-b{e_{y}}=\frac{r_{1}-r_{2}}{r_{1}}(\dot{y}_{1}+by_{1})-b{e_{y}}$;

this implies $\dot{y}_{1}\frac{r_{2}}{r_{1}}-\dot{y}_{2}=-b(\frac{r_{2}}{r_{1}}y_{1}-y_{2})$, or, $\dot{e_{y}^{*}}=-b{e_{y}^{*}}$, where the error variable is redefined by ${e_{y}^{*}}=y_{1}\frac{r_{2}}{r_{1}}-y_{2}$.

One can obtain the time derivative of the Lyapunov function for V($e^{*}_{y}$) as

[TABLE]

Therefore, under the conditions ${\epsilon_{2}}=1$ and $\epsilon_{1}>0$, a new coherent state, ($x_{1}=x_{2}$, $y_{1}=\frac{y_{1}}{y_{2}}$ and $z_{1}=z_{2}$) emerges which we call as a GS state and it is globally stable.

Rössler system: From (9) and (10) the coupled Rössler system, the coupling profile appears as a bidirectional self-coupling that involves $y_{1,2}$ variables. The cross-coupling involves $z_{1,2}$ variables and added the dynamics of $x_{1}$. The dynamical equations of two coupled Rössler systems,

[TABLE]

The error dynamics is

[TABLE]

We determine the stability of $x_{1}=x_{2}$ and $y_{1}=y_{2}$ first. This involves the construction of a Lyapunov function $V(e_{x},e_{y})$ which is positive definite function, $V(e_{x},e_{y})=\frac{1}{2}(e_{x}^{2}+e_{y}^{2}).$ The time derivative of the Lyapunov function,

[TABLE]

For ${\epsilon_{1}}>a/2$ and ${\epsilon_{2}}=-1$, $\dot{V}(e_{x},e_{y})\leq 0$ is negative semidefinite, since we get $\dot{V}(e_{x},e_{y})=0$ if $e_{y}=0$ and for any values of $e_{x}$. However, by using the LaSalle invariance principle LaSalle , the set $S=\{e_{x},e_{y}\}$ does not contain any trajectory except the trivial trajectory $(e_{x},e_{y})$=0, as $\dot{e_{y}}\neq 0$ if $e_{x}\neq 0$. As a result, the trajectory will not stay in the set $S$. So the partial synchronization manifold of $x_{1}=x_{2}$ and $y_{1}=y_{2}$ is asymptotically stable with $t\rightarrow\infty$ above a critical value of $\epsilon_{1}$. We get an unbounded synchronization region for ${\epsilon_{1}}_{c}>a/2$ which is interestingly lower than the critical value found by linear stability analysis (MSF) in Olusola . As a result, the threshold region of self-coupling is expanded by the addition of the cross-coupling link. We consider the error dynamics $\dot{e}_{z}$ separately when considering the effect of heterogeneity.

Assume oscillators have different parameters $b_{1}\neq b_{2}$. The stability in $e_{x}$ and $e_{y}$ is undisturbed by the induced heterogeneity when $\dot{V}(e_{x},e_{y})\leq 0$ remains valid. Now we check the error dynamics ${e_{z}}$,

[TABLE]

From the modified error dynamics of $e_{z}$ one can write $\dot{V^{{}^{\prime}}}(e_{z}^{*})=e_{z}^{*}\dot{e}_{z}^{*}=-ce_{z}^{*2}(1-x_{1}/c)$. Hence, using the LaSalle invariance principle LaSalle the coupled Rössler system (34) is globally synchronized provided $c>\mid x_{1}\mid$. In fact, it shows an error to the limit of $10^{-4}$. Note that this asymptotic stability of $z_{1}=z_{2}$ is valid even for identical systems when $\dot{V}(e_{z})<0$.

Sprott system: We provide analytical evidence for the stability of two coupled Sprott systems. According to (9) and (10), the coupling profile consists of a bidirectional self-coupling involving $x_{1,2}$ variables and a cross-coupling involving $z_{1,2}$ variables is added to $y_{1}$ dynamics. The coupled Sprott system is,

[TABLE]

Error dynamics of the system (37),

[TABLE]

To prove global stability, a Lyapunov function is constructed, $V(e_{x},e_{y},e_{z})=\frac{1}{2a}e_{x}^{2}+\frac{1}{2}e_{y}^{2}+\frac{1}{2}e_{z}^{2}.$ From (38), we examine synchronization in $x$ and $y$ variables of the system (37). We check the stability of $e_{x}$=$e_{y}$=0 first. Now the Lyapunov function is defined as $V(e_{x},e_{y})=\frac{1}{2a}e_{x}^{2}+\frac{1}{2}e_{y}^{2}$. The time derivative of this Lyapunov function is

[TABLE]

The criterion for Eq. (38) to be negative semidefinite, i.e. $\dot{V}(e_{x},e_{y})\leq 0$ are:

[TABLE]

Under the criterion (40) and using the LaSalle invariance principle LaSalle , $x_{1}=x_{2}$ and $y_{1}=y_{2}$ state is stable as $t\rightarrow\infty$. Using the above partial synchronization condition in (38) and after a detuning of $r_{1,2}$ parameter, we obtain modified error dynamics,

[TABLE]

where the modified error variable is,

[TABLE]

and when the time derivative of the Lyapunov function,

[TABLE]

Applying the LaSalle invariance principle LaSalle once again and under the condition (40), the system (37) is globally synchronized in a GS state defined by $x_{1}=x_{2}$, $y_{1}=y_{2}$ and $z_{1}=\frac{r_{1}}{r_{2}}z_{2}$.