# Synchronization optimized networks for coupled nearly identical oscillators and their structural analysis

###### Abstract

The extension of the master stability function (MSF) to analyze stability of generalized synchronization for coupled nearly identical oscillators is discussed. The nearly identical nature of the coupled oscillators comes from some parameter mismatch while the dynamical equations are the same for all the oscillators. From the stability criteria of the MSF, we construct optimal networks with better synchronization property, i. e. the synchronization is stable for widest possible range of coupling parameter. In the optimized networks the nodes with parameter value at one extreme are selected as hubs. The pair of nodes with larger parameter difference are preferred to create links in the optimized networks. And the optimized networks are found to be disassortative in nature, i. e. the nodes with high degree tend to connect with nodes with low degree.

###### pacs:

05.45.Xt,02.60.Pn## I Introduction

Synchronization processes of locally interacting dynamical systems has been the focus of intense research in physical, biological, chemical, technological and social sciences sync_ref (); PecoraPRL1990 (); Heagy1994 (); PecoraPRL1998 (); Arenas2008 (). The simplest and the most studied form of synchronization is the complete synchronization (CS) which is observed when identical dynamical systems are coupled and their state variables become equal as time goes to infinity. As a result of this equality of the state variables the motion of the coupled systems collapses onto a subspace of the overall phase space. This subspace is known as synchronization manifold and the remaining directions in the phase space define the transverse manifold. The complete synchronization is stable when all perturbations in the transverse direction decay with time.

One very important tool for the study of stability of complete synchronization of coupled identical oscillators is the Master Stability Function (MSF), introduced by Pecora and Carroll PecoraPRL1998 (). The MSF is defined as the largest Lyapunove exponent, calculated from a set of equations known as Master Stability Equation (MSE). The MSF simplifies the study of stability for complete synchronization by separating the effect of the network structure from that of the dynamics of individual systems. The MSF also facilitates the stability analysis of complete synchronization for coupled identical oscillators by constructing a single function which can be used to compare synchronizability of different networks.

In practical world, most of the interacting dynamical systems are nonidentical in nature and being nonidentical they cannot exhibit complete synchronization; instead they undergo generalized synchronization Abarbanel1995 (); Abarbanel1996 (). Thus, to better understand synchronization processes of natural systems it is important to construct a master stability function for generalized synchronization. Motivated by this problem, we have extended the formalism of the master stability function to the generalized synchronization of coupled nearly identical oscillators Acharyya2012 (); Acharyya2013 (); Acharyya2014 ().

During the last two decades, the theory of complex network has evolved tremendously compnet_ref (). Many real world complex systems can be modelled as complex networks of interacting dynamical oscillators. Recently, the effect of structure of complex network on the synchronization dynamics of coupled oscillators has been investigated FernandezPRL2000 (); Barahona2002 (); Nishikawa2003 (); Donetti2005 (). In an earlier work, it is shown that the small world scheme enhances synchronizability for a network of coupled identical systems Barahona2002 (). A general argument underlying this phenomenon is that the communication between the coupled systems is more efficient because of the smaller average network distance. In Ref Nishikawa2003 () it is shown that having smaller network distance is not sufficient for performing best synchronizability properties, it is also required to have homogeneous degree distribution among the coupled dynamical systems. In Ref. Donetti2005 (), the authors introduced a new family of graphs, namely the entangled networks, which optimizes synchronizability for many dynamical systems. These entangled networks are interwoven and has an extremely homogeneous structures, i.e. the degree distributions are very narrow. All of the above mentioned results are for networks with identical dynamical oscillators.

The study of finding optimal topology of networks which exhibit better synchronizability has been a subject of paramount interest. Recently, the edge rewiring method has been used vastly for finding optimal topologies of networks for better synchronizability Donetti2005 (); JaliliChaos2008 (); JaliliChaos2010 (). The optimal networks obtained in this way are mostly homogeneous networks, i.e. the degree distributions of these networks are very narrow.

We use the stability criteria provided by the master stability function to construct optimal networks which shows better synchronizability for nearly identical systems. In the optimized networks the nodes with parameter value at one extreme are selected as hubs. The pair of nodes with larger parameter difference are preferred to create links in the optimized networks. And the optimized networks are found to be disassortative in nature, i. e. the nodes with high degree tend to connect with nodes with low degree.

## Ii Master stability function for nearly identical oscillators

We consider a network of coupled nearly identical chaotic oscillators. The dynamics of -th oscillator is given by

(1) |

where, is dimensional state variable and provides the dynamics of the isolated oscillator, is scalar coupling strength. is coupling matrix; if systems and interact then , otherwise and the diagonal elements of are . Thus the elements of satisfy . is a linear coupling function. is some parameter of the dynamics that depends on oscillator . Let , where is a typical value of the parameter and is a small parameter mismatch. We call the parameter as a Node Dependent Parameter (NDP).

When , the coupled oscillators will be identical and for suitable coupling function and coupling parameter the coupled oscillators will undergo complete synchronization, i.e. .

Now, we consider the case where and in this case the synchronization between the coupled oscillators will be of generalized type, i.e. there will be a functional relationship between the variables, .

To determine the stability of generalized synchronization, we do linear stability analysis. Throughout this paper we consider to be small, i.e. . Due to this condition the attractors of the coupled oscillators are not very different from each other. This enable us to expand Eq. (1) in Taylor series about the solution of an isolated oscillator with NDP .

The exponential nature of solution of a linear differential equation is dominated by the homogeneous term of that differential equation. The effect of the NDP appears first in the homogeneous part from the quadratic terms in Taylor series expansion of the function . We retains terms up-to second order in and . The dynamics of deviation is given by

(2) | |||||

As an equation for , the RHS of Eq. (2) contains both homogeneous and inhomogeneous terms. To a first approximation, the inhomogeneity of a linear ordinary differential equations does affect the Lyapunov exponents or the exponential rate of convergence or divergence of the solutions though it can shift the solutions Acharyya2012 (). In nonlinear systems, in addition to the shift the attractor may also deform which can lead to a change in the exponent. This is the case in the desynchronized state. However, in the synchronized state attractors of the coupled oscillators are in generalized synchrony and are related to each other, i.e. . Hence it is reasonable to conjecture that in the synchronized state the shifted solution preserves the nature of the attractor so that the average expansion and contraction rates are not significantly affected Acharyya2012 ().

Hence, to calculate Lyapunov exponent from Eq. (2) we consider the homogeneous equation obtained from Eq. (2)

(3) |

In matrix for Eq. (3) can be written as

(4) |

where is the transpose of the coupling matrix and is is a diagonal matrix with the diagonal entries as the mismatch in NDP. We can see from Eq. (4) that it is necessary to include the quadratic terms in and in the Taylor series expansion as the effect of the NDP is not present in the linear terms.

Let be the -th eigenvalue and be the corresponding right eigenvector of . Define an dimension vector . The dynamics of the vector is given by

(5) |

In general, are not eigenvectors of and hence Eq. (5) is not easy to solve. To solve Eq. (5) we use first order perturbation theory matrix-computation-book () and write Eq. (5) as

(6) |

where is the first order correction and and are the right and left eigenvectors of corresponding to the eigenvalue .

Since both and can be complex, treating them as complex parameters and respectively, we can construct the master stability equation as

(7) |

We call as network parameter and as mismatch parameter.

The master stability function (MSF) is defined as the largest Lyapunov exponent calculated from Eq. (7). The stability of the synchronization is given by the negativity of the MSF. For coupled identical systems, the above equation reduces to the master stability equation given by Pecora and Carroll PecoraPRL1998 () by setting the mismatch parameter .

Here we note that the eigenvalue of the coupling matrix corresponds to the synchronization manifold and the remaining eigenvalues correspond to the transverse manifold. For a given network the synchronization is stable when all Lyapunov exponents corresponding to the eigenvalues of are negative, i.e. they fall in the region where the MSF is negative.

For many chaotic oscillators it is observed that the MSF is negative in a finite interval of the network parameter . Let, the interval be for identical oscillators and for nearly identical oscillators. We can write the condition for stable synchronization of a given network of coupled nearly identical oscillators as

(8) |

The above condition can also be written as

(9) |

### ii.1 Stable interval in coupling parameter

When the variations in the NDP are small the master stability function can be approximated as a linear function near the bifurcation points and and thus one can write

where, and are the mismatch parameters corresponding to the eigenvalues and of the coupling matrix and and are the slopes of master stability function near the points and respectively.

The interval of the coupling parameter where the synchronization is stable then can be written as

(10) |

where, is the stable interval for coupled identical oscillators. We choose as the order parameter to construct optimized networks with better synchronizability.

Finally, we consider an example of component coupled nearly identical Rössler oscillators. The dynamical equations are give as

where, the frequency parameter is the NDP and and are the other Rössler parameters. The MSF for nearly identical Rössler oscillators is shown in Fig. 1. From Fig. 1, we can see that the negative region of the MSF increases as the mismatch parameter increases.

## Iii Synchronization optimized networks

In this section we discuss the construction of optimized networks with better synchronizability for coupled nearly identical oscillators. By better synchronizability we mean that the synchronization is stable for the widest possible interval of the coupling parameter , i.e. is maximum for the optimal networks. We construct the optimized networks with two constraints. The number of vertices and the number of edges of the network are fixed and there are no multiple edges and self loops in the networks. By rewiring the network using Metropolis algorithm we obtain the optimal network. Now, we briefly discuss the Metropolis algorithm for construction of optimized networks with better synchronizability.

Let us start with a connected network of coupled nearly identical oscillators and edges and the coupling matrix be . Let the stable interval of the coupling parameter for this initial network be where the value of is determined using Eq. (10). Now we randomly delete one existing edge and create one new edge at an edge vacancy. Thus, we avoid creating multiple edges and self loops. We reject the resultant network if it is disconnected. Otherwise the stable interval for the resultant network is determined. We accept the resultant network if , otherwise we accept the resultant network with a probability , where and is a temperature-like parameter. This rewiring procedure which defines a Monte Carlo step, is repeated several times. We start with a high value of (). is kept fixed for 1000 Monte Carlo steps or 10 accepted ones, whichever occurs first. Then is reduced by a certain factor () so that stimulated annealing or slow cooling occurs Binder_book (); Binder_book_2 (); numerical_recipes (). We keep on repeating this process until there are no more changes during five successive temperature steps.

## Iv Structural analysis of optimized networks

For our numerical simulation we consider an undirected and unweighted network of coupled nearly identical Rössler oscillators with and the total number of links . The non-identity nature of the coupled Rössler oscillators is introduced by considering the frequency parameter as NDP. The frequency parameter is chosen randomly from an interval . In Fig. 2(a), the order parameter is plotted as a function of the Monte Carlo iterations. The stable interval increases and saturates to a higher value.

In Figs. 3(a), a sample of initial network of 32 coupled nearly identical Rössler oscillators with 174 edges is shown. The vertex size is proportional to frequency parameter , i.e. the vertex with larger has bigger size. In Fig. 3(b), the optimal network obtained from the network of Fig. 3(a) is shown. In Fig. 3(b) we can see that in the optimal network the nodes with higher values of frequencies (), have more connections.

In the optimized network we investigate the structural properties of the network. First, we study the vertices which have more connections than other nodes, i.e. the nodes which are selected as hubs. In Fig. 4, we plot the degrees of the vertices as a function of the NDP for the initial network (blue squares) and the optimal network (red circles). For the initial network all vertices have almost similar degree, but in the optimal network the vertices with larger has higher degree than other nodes.

To quantify this effect, we define the Pearson correlation coefficient between the parameter and degree of a node as,

(11) |

where, is the degree of -the node.

Fig. 2(b) shows (solid line) as a function of Monte Carlo steps. For the random network . We find that increases and saturates to a positive value. Thus, in the synchronized optimized network the nodes which have larger frequencies have more connections and are preferred as hubs. The reason for this is the “V” shape of the stability region in Fig. 1, i.e. the stability range increases as increases.

The degree distribution gives the probability that a randomly chosen node will have degree . In Fig. 4(b), the degree distribution of the initial network (blue dotted line) and the optimal network (red solid line) are shown. The degree distribution of the initial network is Gaussian and has one peak while the degree distribution of the optimal network has a smaller peak at higher degree. The reason for this is the presence of some hubs in the optimal network.

To investigate the question of which edges are preferred, we define the correlation coefficient between the absolute parameter differences between two nodes and the edges as,

(12) |

where, if nodes and are connected and 0 otherwise.

Fig. 2(c) shows as a function of Monte Carlo iterations. We find that increases from 0 (the value for the random network) and saturates. Thus, in the synchronized optimized network the pair of nodes which have a larger relative frequency mismatch are preferred as edges for the optimized network. Again, the reason for this preference of edges is probably the conical shape of the stability region in Fig. 1. The edges are to be chosen so that the parameter increases and the stability region increases.

The clustering coefficient is another important parameter which quantifies the possibility that two neighbors of a common node are also neighbors. The clustering coefficient of vertex is defined as

(13) |

where, is the number of edges that exist among the neighbors of vertex and is the degree of vertex . The clustering coefficient of the entire network is defined as

(14) |

In Fig. 5(a) the clustering coefficient of the network is plotted as a function of the Monte Carlo iterations. From Fig. 5 we can see that the clustering coefficient of the network increases and saturates to a higher positive value. Thus, the optimized network has more local structure than the random network, i.e. there are more triangles than the random network. The result is intuitively easy to understand. Forming a loop will enhance the stability of synchronization due to a faster feedback and smaller the size of the loop better will be the result. We note that the behavior is similar to that for coupled identical oscillators where it has been noticed that networks with larger value of clustering coefficient have better stability of synchronization Donetti2005 ().

Assortative mixing in networks NewmanPRL2002 (); NewmanPRE2003 () gives the tendency of vertices to be connected with vertices of comparable degrees. Let the degrees of verices at the ends of the th edge connecting vertices and be and . Following Ref NewmanPRL2002 () the degree mixing coefficient can be calculated as

(15) |

where, is the total number of edges in the network and the sums are over the edges. When comparable degree nodes get connected the correlation coefficient is positive and the network is called assortative network. The network is called disassortative network when the coefficient is negative. This happen when high degree vertices get connected with low degree vertices. For networks which show no assortative mixing the correlation coefficient is zero. The random networks of Erdős and Rényi and the scale free network model of Barabási and Albert shows no assortative mixing. It has been observed that many naturally evolving networks, such as internet, WWW, protein interaction, neural networks, etc. shows disarrortative mixing of degree NewmanPRL2002 (); ChavezPRE2006 (); Bernardo2005 (); Sorrentino2007 ().

In Fig. 5(b) the degree mixing correlation coefficient is shown as a function of Monte Carlo iterations. The degree mixing coefficient decreases and becomes negative. Thus the optimized network is disassortative in nature. This behavior is consistent with that of coupled identical oscillators..

## V Conclusion

To conclude we have extended the MSF formalism to analyze stability of generalized synchronization for nearly identical oscillators. Using the stability criteria given by the MSF we construct optimal networks with better synchronizability by using Metropolis algorithm. We find the hubs of the optimal networks are nodes with larger frequency. The optimal network is disassortative in nature, i.e. nodes with higher degree tend to connect with node with lower degree.

## References

- (1) A. Pikovsky, M. Rosenblum, J. Kurths, Synchronization, Cambridge University Press, Cambridge, UK, 2001; S. Boccaletti, J. Kurths, G. Osipov, D.L. Valladares, C.S. Zhou, Phys. Rep. 366, 1 (2002); G.V. Osipov, J. Kurths, C. Zhou, Springer, Berlin, Germany, 2007.
- (2) L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. 64, 821 (1990).
- (3) J. F. Heagy, T. L. Carroll, and L. M. Pecora, Phys. Rev. E, 50, 1874 (1994).
- (4) L. M. Pecora and T. L. Carroll, Phys. Rev. Lett. 80, 2109 (1998).
- (5) A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, C. Zhou, Phys. Rep. 469, 93 (2008)
- (6) N. F. Rulkov, M. M. Sushchik, L. S. Tsimring, and H. D. I. Abarbanel, Phys. Rev. E 51, 980 (1995).
- (7) H. D. I. Abarbanel, N. F. Rulkov, and M. M. Sushchik, Phys. Rev. E 53, 4528 (1996).
- (8) S. Acharyya and R. E. Amritkar, EPL 99, 40005 (2012);
- (9) S. Acharyya and R. E. Amritkar, Eur. Phys. J. Special Topics, 222, 939 (2013).
- (10) S. Acharyya and R. E. Amritkar, Synchronization of nearly-identical dynamical systems. I Size instability (to be published)[arXiv:1403.7725].
- (11) Réka Albert and Albert-László Barabási, Rev. Mod. Phys. 74, 47 (2002); M. E. J. Newman, SIAM Rev. 45, 167 (2003); S. Boccaletti, V. Latora, Y. Moreno, M. Chavez , D. -U. Hwang, Phys. Rep. 424, 175 (2006), A.-L. Barabási, and R. Albert, Science 286, 509 (1999),D. J. Watts, and S. H. Strogatz, , Nature (London) 393, 440 (1998).
- (12) L. F. L-Fernández, R. Huerta, and F. Corbacho, and J. A. Sigüenza, Phys. Rev. Lett. 84, 2758 (2000)
- (13) M. Barahona and L. M. Pecora, Phys. Rev. Lett. 89, 054101 (2002).
- (14) T. Nishikawa, A. E. Motter, Y. -C. Lai,and F. C. Hoppensteadt, Phys. Rev. Lett. 91, 014101 (2003).
- (15) L. Donetti, P. I. Hurtado, and M. A. Muñoz, Phys. Rev. Lett. 95, 188701 (2005).
- (16) A. A. Rad, M. Jalili and M. Hasler, Chaos 18, 037104 (2008)
- (17) M. Dadashi, I. Barjasteh and M. Jalili, Chaos 20, 043119 (2010)
- (18) Gene Howard Golub and Charles F. Van Loan, Matrix computations, Johns Hopkins University Press (1983).
- (19) K. Binder, and D. W. Heermann, Monte Carlo Simulation in Statistical Physics: An Introduction, Springer - Verlag, 2010.
- (20) D. P. Landau, and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, Cambridge University Press, 2009.
- (21) W.H.Press and S. A. Teukolsky and W. T. Vetterling and B. P. Flannery,Numerical Recipes in C, Cambridge University Press (1993).
- (22) M. E. J. Newman, Phys. Rev. Lett. 89, 208701 (2002).
- (23) M. E. J. Newman, Phys. Rev. E 67, 026126 (2003).
- (24) M. Chavez, D.-U. Hwang, J. Martinerie, and S. Boccaletti, Phys. Rev. E 74, 066107 (2006).
- (25) M. di Bernardo, F. Garofalo, and F. Sorrentino, arXiv:cond-mat/0506236 (2005).
- (26) F. Sorrentino, M. di Bernardo, and F. Garofalo, Int. J. Bif. and Chaos 17, 2419 (2007).