Partial synchronization phenomena in networksof identical oscillators with non-linear coupling

Partial synchronization phenomena in networks
of identical oscillators with non-linear coupling

Celso Freitas    Elbert Macau Associate Laboratory for Computing and Applied Mathematics - LAC, Brazilian National Institute for Space Research - INPE, Brazil    Arkady Pikovsky Department of Physics and Astronomy, University of Potsdam, Germany Department of Control Theory, Nizhni Novgorod State University, Gagarin Av. 23, 606950, Nizhni Novgorod, Russia
July 5, 2019

We study a Kuramoto-like model of coupled identical phase oscillators on a network, where attractive and repulsive couplings are balanced dynamically due to nonlinearity in interaction. Under a week force, an oscillator tends to follow the phase of its neighbors, but if an oscillator is compelled to follow its peers by a sufficient large number of cohesive neighbors, then it actually starts to act in the opposite manner, i.e. in anti-phase with the majority. Analytic results yield that if the repulsion parameter is small enough in comparison with the degree of the maximum hub, then the full synchronization state is locally stable. Numerical experiments are performed to explore the model beyond this threshold, where the overall cohesion is lost. We report in detail partially synchronous dynamical regimes, like stationary phase-locking, multistability, periodic and chaotic states. Via statistical analysis of different network organizations like tree, scale-free, and random ones, and different network sizes, we found a measure allowing one to predict relative abundance of partially synchronous stationary states in comparison to time-dependent ones.

Synchronization, complex networks, Kuramoto model

Lead paragraph: Large population of coupled oscillators synchronizes if the coupling is attractive and desynchronizes if it repulsive. However, the nature of coupling may depend on the strength of the compelling field: if the force on the oscillator becomes strong, it can switch its behavior from a “conformist” to a “contrarian”. We study such a population on a network. Here the oscillators connected to many others become contrarians first, so that synchrony breaks. We show that the regimes of appearing partial synchrony can be rather complex, with a large degree of multistability, and suggest a network measure which allows to predict relative abundance of static and dynamic regimes.

I Introduction

In a seminal work Kuramoto (1975), aiming to understand synchronization phenomena, Kuramoto proposed a mathematical model of non-identical, nonlinear phase-oscillators, mutually coupled via a common mean field. Studying this system, he identified a synchronization transition to an oscillating global mode when the coupling strength is larger than a critical value, which is proportional to the range of the distribution of the natural frequencies. Over the time, subsequent outcomes based on Kuramoto propositions have shown that his approach can be used as a framework to several natural and technological systems where an ordered behavior (synchronization) emerges from the interactions of a large number of dynamical agents Acebrón et al. (2005); Strogatz (2000). Furthermore, works have shown that the Kuramoto model can be exploited as a building block to develop highly efficient strategies to process information Follmann et al. (2014); Vassilieva et al. (2011).

Recently, generalizations of the Kuramoto model toward interconnections between the elements more complex than the mean field one, have received considerable attention. Indeed, in many real-world problems, each dynamical agent interacts with a subset of the whole ensemble Dörfler et al. (2013); Leonard (2014); Sadilek and Thurner (2014), which can better described using a network topology. A variety of studies have analyzed the onset of the synchronization regime in this context. For a general class of linearly coupled identical oscillators, the Master Stability Function, originally proposed by Pecora and Carrol Pecora and Carroll (1998), allows one to determine an interval of coupling strength values that yields complete synchronization, as a function of the eigenvalues of Laplacian matrix of the coupling graph. For networks of oscillators with non-identical natural frequencies, Jadbabaie et al. Jadbabaie et al. (2004) were able to give similar bounds for the coupling strength of the Kuramoto model without the assumption of infinitely many phase-oscillators. Among related works, Ref. Franci et al. (2010) deals with a model whose natural frequency oscillators change with time, even when they are isolated. Ref. Papachristodoulou et al. (2010) explores the effects of delay in the communication between oscillators. Besides, Ref. Pecora et al. (2013) builds a bridge between graph symmetry and cluster synchronization.

Taking into consideration all of these previous results, one can roughly state that the Kuramoto transition to synchronization always happens if the coupling between oscillators is attractive; while, when it changes to be repulsive, this synchronization state is absent  Tomov and Zaks (2014); Tsimring et al. (2005). However, the structure of the coupling can nontrivially depend on the level of synchrony itself. Such a dependence, called nonlinear coupling scheme, has been explored in recent theoretical Filatrella et al. (2007); Rosenblum and Pikovsky (2007); Pikovsky and Rosenblum (2009); Baibolatov et al. (2009) and experimental Temirbayev et al. (2012, 2013) studies dealing with setup of global coupling. The main effect here is the partial synchrony, which establishes at moderate coupling strengths, where the coupling is balanced between the attractive and repulsive one.

Here, we consider the effects of the non-linear coupling on a network: a set of identical oscillators, which communicate via a connected simple coupling graph. each element is forced by a (local) mean field, which encompasses the oscillators that are connected to it. The coupling function is tailored so that its influence is attractive, if the local acting field is small, or repulsive, otherwise. This coupling strategy implies that only nodes with a large enough number of connections may become repulsive. Thus, the hubs play a key role for the ensemble dynamics. A non-linear coupling parameter in the system tunes the critical quantity of connections and how cohesive this mean field must be in order to allow this transition. So, our approach can be considered as a dynamical generalization of the inhomogeneous populations of oscillators consisting of conformists and contrarians Hong and Strogatz (2011). However, the type of oscillator’s depends on the force acting on it.

Overall dynamics in the model can be qualitatively understood as follows: Let us assume initially that all the mean fields are small. Then, there are only attractive interactions (conformists) in the system. So, in a first moment, they start to mutually adjust their phases. Above a threshold, the most connected oscillators start to feel a repulsive effect that drives them away from the synchronous state. In other words, if an oscillator has a sufficiently large number of neighbors and if it suffers enough cohesive pressure from them, instead of attractiveness, it becomes a contrarian, wishing to be in anti-phase with the force. Then, due to the repulsiveness of some nodes, other mean fields may also become smaller. Finally, this tendency can shift nodes to attractiveness again. As a consequence, an intermediate configuration may emerge due to the balance these conflicting tendencies in the system.

Depending on the non-linear coupling parameter, we report a variety of qualitative dynamic behaviors. In general, for small values of the non-linear coupling parameter, we observed full synchronization and phase-locked states. When this parameter is increased, multistability, periodic and chaotic dynamics take place.

The paper is organized as follows. Initially, we discuss the basic details of the model in section II. In section III, the analytical result about the stability of full synchronization is presented. Numerical experiments in section IV illustrate different possible regimes that the present model can display. Finally, in section IV.3, we perform a numerical exploration to address the correlation of stationary phase locking states with partial synchronization with the network parameters, by exploring different network topologies and sizes.

Ii Model of oscillator network with nonlinear coupling

Let us consider a system of identical phase-oscillators represented by coupled through a simple and connected undirected graph . The dynamics for the -th oscillator is given by the following ordinary differential equation


where denotes the set of neighbors of in the coupling graph . Equations (1) are formulated in the reference frame rotating with the common frequency of the oscillators, so that the latter one does not appear in the equations. The time is normalized by the linear coupling strength. The main feature of our model (1) is the nonlinear coupling parameter, , which modifies the coupling at each node, with a standard setup of the Kuramoto model on a network is recovered. This parameter enters as a coefficient at the square of the norm of the local order parameter:


which measures the magnitude of the force acting on oscillator with index . On the other hand, we represent the (global) order parameter by


where is its norm and is its phase.

We stress that is not normalized (in the sense that there is no division by the number of the terms in the summation, like in ), as it measures the total action of the neighbors on the -th oscillator, which is called local mean field. Simple calculations show that , where denotes the degree of the -th vertex, that is, the number of incoming or outgoing connection, since the graph is undirected. Thus, a necessary condition for a node to suffer repulsive coupling, i.e. , is that .

The introduced order parameters and are maximal in the case of full synchronization , whiles they decrease when oscillators begin to move apart from each other.

If , where , then all oscillator will attract each other so that the full synchronization is established. Next, if becomes larger than , the corresponding oscillator begins to be repulsive related to its local mean field, and the full synchronization breaks. As a result, may decrease and switch again the node to be attractive. Depending on the coupling graph , on the initial conditions , and on the intensity of the nonlinear coupling parameter , numerical simulation reveals that model (1) can exhibit different qualitative behaviors.

If the largest degree in the coupling graph satisfies , with , then no node can be repulsive. We demonstrate in the next section via the Lyapunov analyses, that in fact this condition guarantees that the full synchronized state is stable.

Iii Stability of full synchronization

The basic steps to obtain the results in this section follow the main ideas from Jadbabaie et al. (2004). We begin presenting some preliminary concepts, including elements of the graph theory needed, and a generalized norm of the order parameter to define our Lyapunov function.

Let be the directed incidence matrix of a graph . Thus, is a matrix with rows and columns, where is the number of directed edges of the matrix. The number of undirected edges, \ieignoring the direction, equals is . The columns of represent the edges of the graph: if the -th arrow (directed edge) of the graph goes from to , then the -th column of is zero, except at positions and , where and . Regarding the dynamics of the system, an arrow from node to node in the graph means that node influences node . Although the directed incidence matrix is generally defined for directed graphs, it must be emphasized that only undirected graphs are considered here. We abuse terminology and identify a graph with its adjacency matrix, which is an matrix where ; , if there is and edge between nodes ; and , otherwise. So, . Another common characterization of a graph is its Laplacian matrix, . One can check that . A simple illustration of these concepts is given at Fig. 1.

Figure 1: Example of graph with and , its directed incidence matrix and its Laplacian matrix .

The usage of the directed incidence matrix allows us to rewrite model (1) in a vector form as follows:


where , and stands for the matrix with the elements of a vector on the leading diagonal, and elsewhere.

The square of the global order parameter can be expressed as

However, to build our Lyapunov function, we define a generalized norm of order as


Note that requires the sum of all with (for ), but its generalization takes the sum () only through the edges of the graph. In the case of full coupling graph, direct substitution yields that both global and generalized norm of the order parameter have the same expression .

For any connected symmetrical coupling graph, one can check that the maximum of is the unit, and that if and only if this value is achieved 111On the other hand, the minimum value of does not necessarily correspond to ..



be a candidate Lyapunov function. It is clear that the minimum value of corresponds to the maximum value of , which is equivalent to the fully synchronized state.

In fact, algebraic manipulations reveal that


and that the differential of is given by


As a result, we synthesize in the next theorem the previous intuitively stated idea that if is small enough then full synchronization is a robust phenomenon related to small perturbations over initial conditions.

Theorem 1.

In Model (1), if is smaller than a critical value , then the synchronized stated () is Lyapunov stable.


Consider the potential field defined in (6). So, using the vector form of the model (4) and the expression of the differential from (8), we have that equals to

If we set , then we have that is larger or equal than . Moreover, we can also define a lower bound for , since ; where is the algebraic connectivity of the graph. In the last inequality we used that and that both matrices and have the same non-trivial eigenvalues , where is strictly larger than zero because the coupling graph is connected Godsil and Royle (2001). Therefore,

As a result, implies that , then the fully synchronized state is stable. ∎

Iv Dynamics of partially synchronous states

In this section numerical simulations are performed to illustrate the rich repertoire of behaviors that model (1) may exhibit, specially beyond the threshold , where Theorem 1 cannot be applied.

iv.1 Quantification of dynamical regimes

The numerical integration scheme applied is a forth order Adams-Bashforth-Moulton Method (see Burden and Faires (2010)) with discretization time step . We calculate the partial synchronization metric s from Ref. Gómez-Gardeñes et al. (2007), which, for every two oscillators in the network, takes values , indicating how much the mean phase difference between and varies in the time interval , with . This metric is defined as

One can check that if for some constant , then the exponent in the previous integral is constant and . Nevertheless, if assumes every possible value over the unit circumference with not clear trend, then is close to zero. Now, we average contributions of all neighbor oscillators under a graph with nodes to write

where is the quantity of undirected edges in the graph.

To exclude transients and to detect the statistically stationary state, we adopted the following procedure. For all experiments the time interval is always considered as transient time. Then, the numerical integration is performed in the subsequent intervals , with , until the first such that , or is achieved. Only such a time interval is regarded as non-transient. For the subsequent analysis, we use values of the phases in the stationary time interval regime (whose beginning is shifted to without loss of generality) at points .

iv.2 Examples of complex behaviors

As it was claimed before, in dependence on the network structure, very different types of the dynamics are possible. In order to give impression on it, we present simulations of model (1) with two different coupling graphs displayed as inserts in Fig. 2. Both networks have nodes and they differ only by the rewiring of a single edge. We performed simulations for random initial conditions chosen with uniform distribution over for each experiment. For all these initial conditions , the norm of the order parameter , according to (3), is computed from the time series. As explained in the previous section, in these calculations a transient time is eliminated and that a statistically stationary regime of units of time and points is considered. Then, also for each distinct initial condition, the maximum, average and minimum values of the associated norm of the order parameter are computed, respectively denoted by ; ; and . Of course, converges to a constant if and only if . Now, having different simulations for a fixed coupling graph, we evaluated the maximum, average and minimum value of the average value of the norm of the order parameters over this ensemble, respectively denoted by ; ; and . So, if the norm of the order parameter converges to the same value for all initial conditions simulated, then . For the cases where the norm of the order parameter does not converge over all initial conditions, it will be useful to examine the overall maximum and overall minimum values of the norm of the order parameter, respectively denoted by ; and . Thus, if there is no fixed phase synchronization for all the initial conditions simulated, but the norm of the order parameter presents only small deviations around a common value, then the gap between and is also small. Also notice that , since for all initial conditions. Finally, the maximum Lyapunov exponent for each initial condition is also computed, according to the algorithm in Alligood et al. (1997). The maximum Lyapunov exponent over all the chosen initial conditions is also analyzed.

Figure 2: Numerical results for Model (1) as a function of , for the coupling graphs despited as insect, including 10 random initial conditions. A black line corresponds to , while the interval between and is shown as a gray strip. The gap between and is shown as an orange strip. Since the orange strip is by construction larger or equal than the gray one, the first one is not displayed in the figure when they coincide. Left vertical axes show values related to norm of the order parameter, while the right ones represents the maximum Lyapunov exponent , shown as a red dashed line. Letters in green vertical lines from the upper experiment correspond to subfigures in Fig. 3.
Figure 3: Evolution of for different values of indicated in green at the upper experiment from Fig. 2. Every color represents a different initial condition, while pairs of solid/dashed lines with the same color correspond to solutions whose initial conditions differ not more than at each coordinate. (a) : full synchronization; (b) : fixed phase synchronization; (c,d,e) respect.: examples of multi stability; (f) : example with .

We now describe different regimes observed in the networks, using also Fig. 3, where we depict time series of for some particular choices of , indicated as green letters in the upper panel from Fig. 2 (this is the case we choose for illustrating different regimes). Notice that in both cases, so Theorem 1 guarantees that for the full synchronization state, , is locally stable as illustrated in Fig. 3 (a) (with ).

Panel (a) in Fig. 3 illustrates full synchronization in the network for . For slightly bigger than , simulations suggest that a stationary regime of partial phase synchronization, where , is locally stable as shown in Fig. 3(b) (). Details of this state are clear from Fig. 4. There we show the that the synchronization between the individual oscillators is complete if measured by quantity , and all the oscillators have the same frequency. However, the oscillators are split into two groups with a constant phase shift between them; this division originates in the edge which connects the two largest hubs in the network (vertexes ).

Figure 4: Example of group formation. Details of one of the trajectories from Fig. 3 (b) . On the left side, the coupling graph with in its edges is shown. On the right side, a histogram of is presented with color code representing the normalized frequency.

For larger values of , the regimes are still static but with multistability. For instance, at (see Fig. 3 (c)) two stable configurations emerge with , with (black) or (blue), depending on the initial condition. Fig. 5, which is analogous to Fig. 4, shows the existence of three subgroups, whose members may vary according to the initial condition.

Figure 5: Example multi-stability with group formation. Details of two trajectories from Fig. 3 (c) , the left picture corresponds to the solid black line and the right one to the solid blue. Histograms of are presented with color code representing the normalized frequency.

Other types of multistabilities appear for instance at and , as illustrate in Figs. 3 (d,e). For (panel d) some initial conditions do no converge to a fixed phase synchronization, but to a regime where the order parameter is periodic in time. For (panel c), the norm of the order parameter of all trajectories simulated becomes periodic. Fig. 6 provides an illustration of this regime.

Figure 6: Example of periodic norm of the order parameter. Details of one of the trajectories from Fig. 3 (d) . On the left side, the coupling graph with in its edges is shown. On the right side, a histogram of is presented with color code representing the normalized frequency. We denote by the argument of the order parameter. The bottom picture shows that the curve is closed.

Finally, for (Fig. 3 (f)), one observes a chaotic state with , the distribution of phases and frequencies is illustrated in Fig. 7.

If , we also obtained multistability, with the coexistence of solutions converging to phase-lock and irregular order parameter after the transient, similar to Fig. 3(d).

Figure 7: Example of trajectory with . Details of one of the trajectories from Fig. 3 (f) . On the left side, the coupling graph with in its edges is shown. On the right side, a histogram of is presented with color code representing the normalized frequency.

Now, we compare the results for two slightly different networks depicted in panels (a) and (b) in Fig. 2. The interval of values of with fixed phase synchronization for all initial conditions simulated is very similar for both networks, namely ; also multistability of static partial synchronous regimes have been observed in both cases.

When , contrary to case (a), we observed that the solution for all initial conditions converged to the same phase-lock regime, similar to Fig. 3 (b).

In conclusion of this section, in Fig. 8 we show simulation results for two other networks. Panel (a) shows a random network with nodes and undirected edges. Here predominantly static regimes are observed, only in small ranges of coupling constant chaos with a positive Lyapunov exponent appears. Static regimes, however, demonstrate a large degree of multistability. In panel (b) we show a scale-free network with nodes and undirected edges. Here static states are rare, typically irregular regimes with low values of the order parameter are observed.

Figure 8: Numerical results for Model (1) as a function of , for the coupling graphs depicted as insect, including 10 random initial conditions. The legend of the pictures is the same as in Fig. 2


iv.3 Dependence of partial synchronization regimes on network structure

We have seen that partially synchronous states can be rather different even for similar networks. It is therefore difficult to make general predictions for a relation between the network properties and the dynamical behaviors. Here we attempt such a description, focusing on the property of abundance of static regimes in comparison to time-dependent ones. For this purpose we define the convergence index as the ratio of values of such that converges to a constant value, considering all the random initial conditions. So, both networks Fig. 2 have similar values of the index: in case (a) while in case (b). In contradistinction, network shown in Fig. 8(a) has very large value of the index , while that in Fig. 8(b) a rather low value .

In order to explore which features of the coupling graph are related with , we performed numerical experiments with three sets of graphs, with nodes. Each set consists in three common types of networks, each one with members, generated as: (i) random (Erdös-Rényi) graphs with edges; (ii) scale-free graphs, also with ; and (iii) tree graphs ( edges). The Barabási-Albert algorithm is applied for the last two types of networks (ii),(iii), with an initial clique of nodes and with other nodes been connected to existing ones. For the -edges scale-free graphs, we fixed and ; while for the tree graphs ( edges scale-free graphs), . We point out that all graphs created are connected and symmetrical. Additionally, three sets of initial conditions , with uniform distribution over and , have been explored. So, for each of the coupling graphs we computed its correspondent values by numerical integration of model (1) for .

In Table 1 we report the mean value and the standard deviation of for each topology and size of coupling graph. From these data we see that the mean value of increases if we go from tree to scale-free and to random graphs, respectively. However, these difference becomes less noticeable for larger values of . Both the mean value and the standard deviation of decrease with larger networks.

Tree 0.421 (0.260) 0.016 (0.006) 0.008 (0.004)
Scale-free 0.857 (0.029) 0.050 (0.015) 0.013 (0.003)
Random 0.872 (0.090) 0.183 (0.063) 0.077 (0.022)
Table 1: Mean value of and its standard deviation (in brackets) for each network type and size simulated.

We have explored different networks metrics, searching for one mostly correlated with the convergence index . Let denote the Laplacian eigenvalues of the coupling graph Godsil et al. (2001). Recall that this graph is assumed to be simple and connected. We stress that these eigenvalues express fundamental characteristics of the graph. For instance, is related with graph diameter and with its largest degree size.

We found that the quantity , defined as the ratio between the maximum eigenvalue and the average of the non-trivial eigenvalues of the Laplacian matrix of the graph, is rather suitable for this purpose. Formally, it is defined as

Figure 9: Convergence index versus . Networks with nodes are represented as circles, squares and triangles, respectively. The two experiments from Fig. 2 are shown as disks. We show in red an exponential fit for the data.

In Fig. 9 a correlation plot between and measure for the correspondent graph is presented. From there, we observe a clear trend indicating that the greater the value of is, the smaller is the value of . Independently of the network type and size, static regimes of partial synchronization, full synchronization and phase-lock, are typical for values of , like in the experiments from Fig. 2. On the other hand, graphs with larger values of this measure yields more irregular dynamics, like time-dependent periodic and chaotic regimes, as the ones from Fig. 8.

V Conclusion

In this work we introduced and studied a Kuramoto-like model of identical oscillators with non-linear coupling. Our main parameter was , which governs the coupling nonlinearity strength. It is clear that the most influence of nonlinearity in the coupling is on the hubs which experience strong forcing from many connected oscillators, while less connected nodes may still operate in a linear-coupling regime.

We proved that if this parameter is smaller than the inverse of the square of the maximum vertex degree in the network, then the full synchronized state is stable. Via numerical experiments, we showed that our model can display a variety of other qualitative behaviors of partial synchronization, like stationary phase locking, multistability, periodic order parameter variations, and chaotic regimes. We explored the relative abundance of stationary phase locking regimes under different network topologies. Statistical analysis performed suggests that tree graphs are much less likely to exhibit stationary phase locking in comparison with scale-free or random networks. In addition, this type of behavior becomes more rare if we increase network sizes, irrespective to the network topology. Finally, we also found a good correlation between the ration between the maximum eigenvalue and the average of the non-trivial eigenvalues of the Laplacian matrix of the graph, and the proportion of the repulsion parameter values which yield stationary phase locking. The greater this measure is, the smaller tend to be presence of stationary phase locking states in the system.

As a future research, we want to investigate analytical conditions and correlations involving other graph measures related to other forms of synchronization in the model.

We would like to thank the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - CAPES (Process: BEX 10571/13-2) for financial support. AP V.V. thanks the IRTG 1740/TRP 2011/50151-0, funded by the DFG/FAPESP, CNPq, and the grant/agreement 02.B.49.21.0003 of August 27, 2013 between the Russian Ministry of Education and Science and Lobachevsky State University of Nizhni Novgorod.


  • Kuramoto (1975) Y. Kuramoto, in International Symposium on Mathematical Problems in Theoretical Physics, Lecture Notes in Physics, Vol. 39, edited by H. Araki (Springer Berlin Heidelberg, 1975) pp. 420–422.
  • Acebrón et al. (2005) J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort,  and R. Spigler, Rev. Mod. Phys. 77, 137 (2005).
  • Strogatz (2000) S. H. Strogatz, Physica D: Nonlinear Phenomena 143, 1 (2000).
  • Follmann et al. (2014) R. Follmann, E. Macau, E. Rosa,  and J. Piqueira, Neural Networks and Learning Systems, IEEE Transactions on PP, 1 (2014).
  • Vassilieva et al. (2011) E. Vassilieva, G. Pinto, J. Acacio de Barros,  and P. Suppes, Neural Networks, IEEE Transactions on 22, 84 (2011).
  • Dörfler et al. (2013) F. Dörfler, M. Chertkov,  and F. Bullo, Proceedings of the National Academy of Sciences 110, 2005 (2013).
  • Leonard (2014) N. E. Leonard, Annual Reviews in Control 38, 171 (2014).
  • Sadilek and Thurner (2014) M. Sadilek and S. Thurner, arXiv preprint arXiv:1409.5352  (2014).
  • Pecora and Carroll (1998) L. M. Pecora and T. L. Carroll, Physical Review Letters 80, 2109 (1998).
  • Jadbabaie et al. (2004) A. Jadbabaie, N. Motee,  and M. Barahona, in American Control Conference, 2004. Proceedings of the 2004, Vol. 5 (2004) pp. 4296–4301 vol.5.
  • Franci et al. (2010) A. Franci, A. Chaillet,  and W. Pasillas-Lépine, in Decision and Control (CDC), 2010 49th IEEE Conference on (IEEE, 2010) pp. 1587–1592.
  • Papachristodoulou et al. (2010) A. Papachristodoulou, A. Jadbabaie,  and U. Munz, Automatic Control, IEEE Transactions on 55, 1471 (2010).
  • Pecora et al. (2013) L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy,  and R. Roy, arXiv preprint arXiv:1309.6605  (2013).
  • Tomov and Zaks (2014) P. Tomov and M. Zaks, in Nonlinear Dynamics of Electronic Systems (Springer, 2014) pp. 246–253.
  • Tsimring et al. (2005) L. S. Tsimring, N. F. Rulkov, M. L. Larsen,  and M. Gabbay, Phys. Rev. Lett. 95, 014101 (2005).
  • Filatrella et al. (2007) G. Filatrella, N. F. Pedersen,  and K. Wiesenfeld, Phys. Rev. E 75, 017201 (2007).
  • Rosenblum and Pikovsky (2007) M. Rosenblum and A. Pikovsky, Phys. Rev. Lett. 98, 064101 (2007).
  • Pikovsky and Rosenblum (2009) A. Pikovsky and M. Rosenblum, Physica D: Nonlinear Phenomena 238, 27 (2009).
  • Baibolatov et al. (2009) Y. Baibolatov, M. Rosenblum, Z. Z. Zhanabaev, M. Kyzgarina,  and A. Pikovsky, Phys. Rev. E 80, 046211 (2009).
  • Temirbayev et al. (2012) A. A. Temirbayev, Z. Z. Zhanabaev, S. B. Tarasov, V. I. Ponomarenko,  and M. Rosenblum, Phys. Rev. E 85, 015204 (2012).
  • Temirbayev et al. (2013) A. A. Temirbayev, Y. D. Nalibayev, Z. Z. Zhanabaev, V. I. Ponomarenko,  and M. Rosenblum, Phys. Rev. E 87, 062917 (2013).
  • Hong and Strogatz (2011) H. Hong and S. H. Strogatz, Phys. Rev. Lett. 106, 054102 (2011).
  • (23) On the other hand, the minimum value of does not necessarily correspond to .
  • Godsil and Royle (2001) C. Godsil and G. Royle, Algebraic Graph Theory, Graduate Texts in Mathematics., Vol. 207 (volume 207 of Graduate Texts in Mathematics. Springer, 2001).
  • Burden and Faires (2010) R. L. Burden and J. D. Faires, Numerical Analysis (Brooks Cole, 2010).
  • Gómez-Gardeñes et al. (2007) J. Gómez-Gardeñes, Y. Moreno,  and A. Arenas, Phys. Rev. Lett. 98, 034101 (2007).
  • Alligood et al. (1997) K. Alligood, T. Sauer,  and J. Yorke, Chaos: An Introduction to Dynamical Systems, Chaos: An Introduction to Dynamical Systems (New York, NY, 1997).
  • Godsil et al. (2001) C. D. Godsil, G. Royle,  and C. Godsil, Algebraic graph theory, Vol. 207 (Springer New York, 2001).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description