Using synchronism of chaos for adaptive learning of network topology
In this paper we consider networks of dynamical systems that evolve in synchrony and investigate how dynamical information from the synchronization dynamics can be effectively used to learn the network topology, i.e., identify the time evolution of the couplings between the network nodes. To this aim, we present an adaptive strategy that, based on a potential that the network systems seek to minimize in order to maintain synchronization, can be successfully applied to identify the time evolution of the network from limited information. This strategy takes advantage of the properties of synchronism of chaos and of the presence of different communication delays over the network links. As a motivating example we consider a network of sensors surveying an area, in which information regarding the time evolution of the network connections can be used, e.g., to detect changes taking place within the area. We propose two different setups for our strategy. In the first one, synchronization has to be achieved at each node (as well as the identification of the couplings over the network links), based solely on a single scalar signal representing a superposition of signals from the other nodes in the network. In the second one, we incorporate an additional node, termed the maestro, having the function of maintaining network synchronization. We will see that when such an arrangement is realized, it will become possible to effectively identify the time evolution of networks that are much larger than would be possible in the absence of a maestro.
In recent years, much effort has been devoted to the study of synchronization of networks of coupled dynamical systems Winfree (1980); Kuramoto (1984); Ermentrout and Kopell (1984); Strogatz (2001); Strogatz and Stewart (1993); Pecora and Carroll (1998); Barahona and Pecora (2002); Boccaletti et al. (2006); Nishikawa et al. (2003); Motter et al. (2005); Moreno and Pacheco (2004); Oh et al. (2005); Lee (2005); Yook and Meyer-Ortmanns (2005); Zhou and Kurths (2006); Wang (2002); Wang and Chen (2002); Parlitz (1996); Maybhate and Amritkar (2000); Hegazi et al. (2002); Chen et al. (2004); Lu and Cao (2005); Atay (2003); Atay et al. (2004). Important results have been obtained under the assumption that communication between the connected systems is immediate; i.e., the signal from system at time is received by systems at the same time . Most of these studies have also focused on the case where the network topology is static, i.e., the coupling strengths between the connected systems is fixed and does not evolve in time. Yet in real situations networks of dynamical systems are often time-varying and affected by communication delays. In this paper we will address situations incorporating both the time varying structure common in typical real networks and the delays which inherently affect communication across the network. Our main emphasis will be on the situation where the time-dependent connection strength from system to another system in the network is unknown by system , and we investigate how system can deduce these strengths by an adaptive strategy making use of the phenomenon of synchronization of chaos. As we shall show, diversity in the delay time associated with different links, as well as chaos of the exchanged signals, will play a crucial role in achieving this goal.
As a possible motivation, we will consider as our reference application that of a network of sensors that are dispersed over a ground and interact with each other through direct wireless communication. Namely, we will assume that the sensor network is employed to detect the evolution of the configuration within the surveyed area, e.g., due to motions of objects on the ground. In that case, a movement of an object across a communication pathway affects the strength of the signals exchanged by the nodes at its endpoints. Therefore, if we can succeed in identifying the time-evolution of the network couplings, this will provide information potentially enabling us to localize the position of the objects moving across the ground. We consider that each of the dynamical systems (‘nodes’) are identical chaotic systems broadcasting chaotic signals. As the network evolves, each node uses an adaptive strategy (which we specify) to maintain synchronization, and, in so doing, it determines the strengths of the couplings to it.
Moreover, we will also show how our adaptive strategy could greatly benefit from the introduction of a special node, having the specific function of maintaining network synchronization. In what follows, we will refer to this particular node as the maestro, in order to distinguish it from the remaining nodes of the network, which we will refer to as the orchestra. For example, for the case of our reference application, we will make the assumption that this node is positioned in such a way that the strength of the signal sent by this node is not affected by the movements of the objects on the ground (e.g., it might be located atop a tower, hill or tall building). We will see that when such an arrangement is realized, it will become possible to effectively identify the time evolution of much larger networks than would be possible in the absence of a maestro.
The problem of identifying the structure of real complex networks, arises in all those situations where the exact network topological structure is unknown or uncertain, therefore it is of relevance in various fields of science and technologies. Examples are metabolic networks, biological neural networks, electric power grids, and so on Sales-Prado et al. (2007); Jin et al. (1999); Goto et al. (1998). Moreover, in many real world situations, there is often an intrinsic difficulty in extracting information on the network structure. Hence, the problem of learning the network topology, represents a challenge in a number of real applications, that involve both natural and man-made systems Sales-Prado et al. (2007); Jin et al. (1999); Goto et al. (1998); Zhou and an Lu (2007).
In Sec. II we present our proposed adaptive strategy for maintaining synchronization of chaos and determining network topology. We do this both for the case with no maestro and for the case with a maestro. We also show that it is essential for our technique that signals be chaotic and experience time delays along the links of the network. In Sec. III we test our strategy through numerical experiments. We show that the presence of a maestro can greatly benefit our strategy. Conclusions and further discussion are presented in Sec. IV.
Ii Adaptive strategy
In a previous paper Sorrentino and Ott (2008) we showed how a time-varying complex network can be synchronized by using a ‘potential’ that each network node seeks to minimize. In this paper, we will extend this concept with the goal, of not only maintaining synchronization, but of also identifying the time-varying evolution of the network itself. In Refs. Yu and Cao (2006); Wu (2008) it was shown how to identify the static structure of a complex network of synchronizing systems by introducing a response (replica) network that can adjust its couplings to converge to those of the original one. The main focus of our paper will be on identifying the evolution of time-varying couplings of a complex network from limited information. We will show how this may be achieved by solving a small number of auxiliary differential equation, provided that some general conditions on the time delays of the network links are met. We will consider two versions of our problem. In Version I the maestro is absent, while in Version II we have a maestro.
In Version I, the set of dynamical equations describing the dynamical state vector of dimension at each node of the time evolving network is the following:
is the incoming signal at node , which is a linear combination of the different scalar transmitted signals , through the time varying coupling coefficients . Here is a constant -vector specifying the coupling of the scalar received signal to node . (In our numerical experiments in Sec. III, we will take the components of to be all zero except for one component, say , for which .) We assume that the time delays , between each pair of connected nodes and , are known, e.g., the sensors occupy fixed positions and thus the communication delays can be easily obtained from knowledge of the relative distances between the sensors. In contrast, we assume that the strengths of the network interactions, represented by the , are unknown and variable in time. Here, the ’s might represent the spreading, scattering, and attenuation of a wireless output signal broadcast by node to node . In the case where the input to from is accomplished by propagation of an electromagnetic wave from to , the link strength will vary with time as a body moves through the propagation path from to . Although the vary with time, we shall be most interested in the case where the time scale for this variation is slow compared to the natural time scale of the chaotic dynamics of an individual uncoupled system,
The quantity , appearing in Eq. (1) represents an estimate made at node of the true coupling . How this estimate can be made, as well as the accuracy of the estimate, will be the main focus of this paper.
In Version II, we consider an additional dynamical state vector, i.e., , associated with the maestro node, and the following set of dynamical equations:
and the maestro obeys the same dynamical evolution as an uncoupled node, i.e.,
Here is the strength of the coupling of the maestro to each node ; is the communication delay between the maestro and node ; is a parameter that tunes the relative strength of the coupling exerted on the network nodes by the maestro, with respect to that of the other nodes in the network (the orchestra). When , each node is affected only by the forcing by the maestro, and therefore it can be seen as an isolated slave system connected to the master. On the other hand, when , each node will not feel the influence of the maestro and synchronization can only be realized from the mutual coupling with the other nodes in the network (the orchestra); thus the case corresponds to our problem in Version I. Values of between and (the case of most interest) correspond to intermediate situations where each node feels both the influence of the maestro, as well as that of the orchestra.
In what follows we assume that is known a priori at each node . In Sec. III.D we also address the case that the assumed known coupling with the maestro is affected by a small error.
Thus, in a similar way to the case of a network of coupled oscillators without delays or time variability of the coupling Pecora and Carroll (1998), the synchronization manifold,
Once the dynamical functions and are assigned, the stability of synchronized states to perturbations from the synchronization manifold will depend on the dynamical evolution of the network adjacency matrix , on the overall coupling , and on the time delays (plus on and in the case of Version 2). In what follows, we will proceed under the assumption that the synchronization manifold (8a - 8b) is stable for our choices of the time evolution of , , and the time delays and , and we will focus on the problem of identifying the network structure (i.e., obtaining the estimates ).
We assume that at each node direct information on the couplings is unavailable. Thus, the estimates at node must be obtained solely using the aggregate incoming signal . In order to make this estimate, we introduce a ‘potential’ at each node :
If the time scale for variation of is larger than , then can be approximated by
As given by (10a) and (10b), the potential is approximately the time averaged squared synchronization error, and is the temporal extent over which the averaging is performed. Note from (10a) that if . Similarly, from (10b) we note that if and . Due to the chaotic nature of the ’s, provided that the are all distinct, these conditions can only be satisfied if (i.e., our estimates of the time evolutions of the network couplings are correct) and (i.e., the network systems are synchronized in time). Thus we seek to minimize . Therefore, we impose that the evolve according to the following gradient descent relation:
and is a constant larger than zero.
In practice it is useful to avoid doing the above integrations at each time step of our adaptive computation. Rather, we note that , , and given by (12) and (13), and (14) satisfy simple differential equations, and we employ these in our adaptive computation of :
The simplest version of this setup, used in our numerical example in Sec. III, corresponds to letting , in which case Eqs. (11) implies that we can obtain the by solving the system of linear equations,
According to our previous discussion (in particular, the condition for the approximate equality in Eq. (10a)), in order for our identification strategy to work properly we should choose such that it is in the intermediate range,
where is the time scale of the chaotic dynamics of an uncoupled nodal system (), and is the time scale on which the network itself changes (i.e. the time scale on which varies). We will subsequently give quantitative definitions of and for our examples.
We emphasize that we try to identify the time evolutions of the several ’s, by using only one piece of external information at each node , namely, the incoming signal . Our ability to extract estimates for the values of () at each node from each incoming signal, relies on the temporal signal diversity of due to the presence of both chaos and different time delays along the network links. To illustrate the necessity of different time delays, consider the following example. Assume, for instance, that our network consists of three coupled systems in the Version I scenario ( in Eq. (1)), and we want to identify the time evolutions of and when the only available information at node 1 is , the time evolution of the incoming signal at node 1. If we are on the synchronization manifold, then Eqs. (1) and (2) yield
and, for , the terms are equal, yielding for the term in (20) that is in curly brackets
Thus, if , we see that we could obtain synchronization at node for any and such that , and we cannot hope to estimate both and only from . In contrast, if , then on the synchronization manifold and the previous degeneracy is removed. Thus we now have the possibility that it may be feasible to obtain the desired two unknowns and ; in other words, it may become possible to obtain good estimates and , provided that . Hence the presence of different delays is crucial in the identification process.
Moreover our strategy benefits from the chaotic nature of the exchanged signals at the network nodes. In particular, we take advantage of the high information content encoded in chaotic signals. As an illustrative counterexample, assume that, in place of the chaotic signals, we had sine waves of the same frequency broadcast by the network systems. The communication delays would result in phase shifts of the sine waves, and the signal received at each node, being the sum of sinusoids of the same frequency, would be a sinusoidal signal. In this case, there would be no information content in the received signal other than its amplitude and phase; therefore, we would not be able to extract more than two pieces of information from it.
Iii Numerical experiments
iii.1 Version I
As a first example, we consider the problem in Version I [Eqs. (1) and (2)], and we seek to estimate the time evolution of a small network of nodes and directed links (two directed links incoming to each node). For our uncoupled dynamical system, , we consider a Rössler oscillator, and
We choose the oscillators to be coupled in the variable, i.e., , , .
The delays associated with links , are given in matrix form, , as follows:
where is defined as the time at which the autocorrelation function of , , decays to , and is a typical orbit on the chaotic attractor of the uncoupled system, . Note that for each , the for every are all unequal. For our choice in Eq. (30), we obtain .
For we assume the following network evolution,
where the are random numbers drawn from a uniform distribution between and , the are random numbers drawn from a uniform distribution between and , and the are random numbers drawn from a uniform distribution between and . We define the network evolution time scale for our example to be .
The network state vectors are initialized as follows,
Here is a randomly chosen point on the Rössler attractor, obtained by evolving the uncoupled system, , until it is on the attractor, and then taking its state at some arbitrarily chosen subsequent random time; and are zero-mean independent random numbers of unit variance drawn from a normal distribution; are the standard deviations of the time evolutions of obtained from numerical solution of an uncoupled oscillator, . For our numerical experiment we take , , and . Using these parameters we have that the inequalities in Eq. (19) are satisfied: . The computations were done over a time interval .
To measure the extent to which our strategy works we introduce the following two error measures:
Figures 1,2, and 3 show results of our computations. Figure 1 shows (Fig. 1(a)) and (Fig. 1(b)) versus at the end of the the run, . For all nodes, , and are identical to within the width of the plotted curves. Figure 2(a) shows that is typically less than over the entire run, while Fig. 2(b) shows that the synchronization error is typically less than . These results confirm that our adaptive strategy is effective in synchronizing the network and that the closely follow the evolutions of the time evolution of the true couplings (see Fig. 3).
with all the other conditions the same as those used in obtaining Figs. 1-3. As can be seen from Fig. 5, for this case our identification strategy fails, with the deviating from the evolution of the . At the same time, we observe from Fig. 4 that the network systems evolve approximately synchronously in time. Thus, in this case we can approximately synchronize the network, but we fail in correctly estimating the evolution of the network couplings.
Unfortunately, by using the strategy in Version I, we were not able to synchronize complete graphs of Rössler oscillators, as the network dimension was increased. As we will explain in the Appendix of this paper, we have an understanding of the reason for this failure. Briefly, a requirement for correctly estimating at any given node, say , is that the ’s associated with are different enough (we address this issue more rigorously in the Appendix). However, if is too large, we cannot make the too separated, because stability of the synchronization manifold (8a) (under the condition ), is lost if the ’s are too large. Thus, for our setup in Version I, achieving the required diversity in the ’s, contrasts with our requirement of having the network synchronize.
We have also tested our Version I scenario for a case as in Figs. 1-3, but with the network evolving on a faster time scale (instead of ), and for this faster evolution we have found that our strategy fails (results not shown).
iii.2 Version II
The failure of our method with large when applied to Version I is the main reason for considering in what follows the alternative set up of our problem in Version II [Eqs. (4) and (5)], which incorporates a maestro node, having the function of maintaining the network synchronizability. Thus we assume that each node/dynamical system is subject to two actions: a coupling exerted by the other nodes of the network, which in what follows we will refer to as the orchestra and a coupling exerted by a master system, which in what follows we will refer to as the maestro. By tuning the parameter in (4) and (5) between and we are able to interpolate between a situation where the nodes feel only the influence of the maestro and a situation where they feel only the influence of the orchestra (corresponding to Version I).
We assume that the maestro coupling is known and time independent, but that the other network couplings, , are time-varying and unknown at the network nodes. In order to guarantee that the synchronization manifold is stable, we take the time delays with the maestro to be small enough. Our hope is that, once synchronization is ensured through the maestro, we will potentially be able to choose the other ’s so that they are sufficiently spread that we have a much greater likelihood of being able to correctly estimate the . As we will show, by making use of the maestro node, we will be able to identify larger networks, whose time-evolution is faster than would be possible in Version I ().
To test our strategy in Version II, Eq. (4), again, as for Version I, we assume that we want to identify the time evolution of globally connected networks of Rössler oscillators, with , , . For our first experiment we take a network of nodes and directed links to be identified. We assume , at each nodes in (4); we have verified beforehand that in the case where for , for this choice of the parameters, the synchronization manifold (8b) is stable. We are interested in identifying the time evolutions of all the ’s associated with the network links.
We take the time delays along the network directed links to be
where is the remainder of the integer division of by , and is a variable parameter to be specified. Note that, for our choice of the time delays in (29), we have that , . Moreover, from (29) we also have that for any node and for any pair of delays associated with any two links incoming to , say and , we have that , and the values and are separated by an integer multiple of . Thus, we do not have equal delays along incoming links to any node.
A first experiment involving our strategy in Version II, is shown in Figures 7 and 8. For this experiment, we took (i.e., synchronization depends mostly on the coupling with the maestro), and set , indicating that the separation time between any pair of network delays is a multiple of . Moreover we assumed that the network evolves much faster than in the cases reported in Figs. 1-6 of Sec. III.A for Version I, that is, , and . The network was evolved for a long time until time , starting from synchronous initial conditions and assuming .
As shown in Fig. 7, and , throughout the length of the run. The black dots in Fig. 8 show the time evolutions of the ’s and the red lines those of the ’s. As can be seen, our adaptive strategy is effective in synchronizing the network with the closely following the evolutions of the . Remarkably, we have succeeded in identifying the time evolutions of all the links present in the network.
The same experiment in Figs. 7 and 8 has also been repeated for different values of . In Fig. 9, the time averaged errors and are plotted as is increased between and (continuous line). As can be seen, our strategy is observed to be effective in identifying the network evolution for values of up to about .
We tested our strategy for increasing network size up to . As shown in Fig. 10, we found that our strategy was effective in identifying the network evolution for . Note that for our choice of fully connected networks, the total number of unknowns to be identified is equal to , e.g., it is equal to when and is equal to when . We expect that, the task of identifying the evolution of the couplings becomes more difficult as the network dimension is increased. Consider for example a globally connected network of nodes. In such a case, each one of the nodes has to identify the time evolutions of different couplings to the remaining nodes in the network, from only one available received signal. Note that, by making the assumption of dealing with fully connected networks, and for which all the links are time-varying, we are considering a worst case scenario; in the case of our sensor application (see Sec. I), we expect that our strategy could benefit from a situation in which only a fraction of the network links might be changing as an object crosses the surveyed region.
We also tried to estimate the time evolution of globally connected networks whose dimension is larger than . For an node network, our strategy failed at identifying the couplings evolution; however, we repeated the same experiment with , but with the network evolving more slowly (i.e., ), and in this case we were able to correctly identify the time evolving couplings.
iii.3 Sensitivity to non-identicality of the network systems
In practical situations (including our reference application of sensor networks) it is impossible to guarantee that all the network systems are precisely identical. Therefore, we have also tested the robustness of our scheme to deviations of the individual systems from identicality. To this end, we replace in Eq. (30) by
where for each node the parameter is chosen randomly with uniform density in the interval . We then repeat our original experiment. The parameter characterizes the degree of non-identicality of the node dynamical systems. In Fig. 11, we show how the time averaged errors and depend on , for a network of nodes, . As can be seen, the more non-identical the systems are, the worse our strategy performs. However, for small values of , for example, for , the synchronization and identification errors are respectively less than and , i.e., and , thus indicating that good results may still be obtained when the coupled systems deviate from being precisely identical.
iii.4 Sensitivity to errors in the values of
The results presented previously rely on the assumption that the coupling with the maestro is exactly known at the network systems. Here we investigate a situation where the assumed known coupling with the maestro has a small error, and we study the effects of this error on our strategy. To this aim we replace Eq. (4) by
where , with is randomly chosen in , and the parameter characterizes the degree to which the real coupling with the maestro deviates from the coupling, , assumed at node . As can be seen from Fig. 12, the larger the error, the worse our strategy performs. However, for small values of , for example, for , the synchronization and identification errors are respectively less than and , i.e., and , thus indicating that good results may still be obtained provided that the coupling with the maestro is affected by a not too large error.
We have also investigated the dependence of our strategy on the parameter for the case in which small mismatches in the coupling with the maestro are present, i.e., . This is shown by the dotted line in Fig. 9. As can be seen, for such a case, our strategy works in a bounded interval of values of the parameter , i.e., , indicating that under these conditions, the best setup is realized when synchronization depends on both the influence of the maestro, as well as that of the orchestra.
In this paper we have dealt with the problem of identifying the time-evolutions of the couplings (assumed unknown) associated with the links of a network by making use of the phenomenon of chaos synchronization. We considered networks with time varying strengths of the couplings (due the effects of external unpredictable factors), as well as different coupling delays affecting the communication between the nodes.
We proposed an adaptive strategy based on the concept of a ‘potential’ that the network systems seek to minimize in order to identify the time evolutions of the network couplings. We considered fully connected networks of nodes and showed how our strategy could be successfully applied to identify the unknown couplings at each of the network nodes. We emphasize that, we achieved this goal by using very little information exchanged at the network nodes: namely the coupling was realized in only one state variable, and the only available information at each node at time was an aggregate incoming signal from the other nodes, (compare for instance our approach with others previously reported in the literature Yu and Cao (2006); Wu (2008)). Here we took advantage both of the phenomenon of synchronization of chaos and of the temporal signal diversity of the signal due to the presence of different time delays over the network links.
We proposed two alternative versions of our strategy. In the first one, synchronization had to be achieved at each node (as well as the identification of the couplings over the network links), based solely on the information received from the other nodes in the network. For this case, we showed that we were able to correctly identify the evolution of the couplings in the case of small networks, i.e., a network of nodes, whose time evolution was sufficiently slow.
In order to extend and improve our results, we introduced a second version of our strategy, where we made use of an additional node, termed the maestro, having the function of maintaining network synchronization. We made the assumption that the coupling with the maestro was known at the other network nodes and that the maestro was not affected by the dynamics of the other nodes in the network. We showed that such an arrangement could lead to much better results than in the case where the maestro was absent. In particular, we were able to identify the time evolution of networks that were larger and more rapidly evolving than was possible in the other case. Other possible setups of our strategy are also possible, e.g., a selective choice of the nodes to be connected to the maestro (pinning control).
In our experiment in Figs. 4, 5, and 6 (Version 1, ) we observed that we were able to synchronize the network on an approximate synchronous evolution obeying (3), but we were not able to correctly identify the time evolutions of the ’s. In what follows, we give an explanation for this finding. Our explanation is that, when the are not large enough, it is possible to obtain approximate (not exact) synchronization for a range of values of the . Although this synchronization is not exact, it is still good enough to imply that the minimum of the potential is so broad as to make our procedure ineffective. To illustrate this, assume the network systems evolve synchronously, i.e., , and that the function and are infinitely differentiable. Expanding in a Mac-Laurin series in the parameter , we have
If this expansion of the function has a finite radius of convergence, then (32) is valid under the assumption that is small enough.
Our strategy seeks to minimize the potential at each node , where may be thought of as a sliding average of the mean squared synchronization error at node , and the synchronization error is
Assuming for all and and making use of (32), we obtain
By construction, the change slowly with time (on the time scale ), as do the (on the time scale ). On the other hand, , and hence change rapidly with time. Thus to make the sliding average of as small as possible, we need
giving zero for the sliding average. However, if is large and the are not too big, we can achieve values of that are very small, even if (35) is far from being satisfied. To see this we note that if the delays are too small, for some sufficiently large , say , the summation over may be truncated with only a small resulting error in . Assuming this to have been done, and setting each term in the sum over in (34) equal to zero for each , we obtain equations for the unknowns . If the number of unknowns is greater than , then we can satisfy all the equations for a continuum of values not satisfying (35), and the potential will then be near zero.
This work was supported by an MURI contract administered through ONR (contract N000014-07-0734).
- Winfree (1980) A. Winfree, The geometry of biological time (Springer, New York, 1980).
- Kuramoto (1984) Y. Kuramoto, Chemical Oscillations, Waves and Turbolence (Springer, Berlin, 1984).
- Ermentrout and Kopell (1984) G. Ermentrout and N. Kopell, SIAM J. Math. Analysis pp. 215–237 (1984).
- Strogatz (2001) S. Strogatz, Nature 410, 268 (2001).
- Strogatz and Stewart (1993) S. Strogatz and I. Stewart, Scientific American 269, 102 (1993).
- Pecora and Carroll (1998) L. Pecora and T. Carroll, Phys. Rev. Lett. 80, 2109 (1998).
- Barahona and Pecora (2002) M. Barahona and L. Pecora, Phys. Rev. Lett. 89, 054101 (2002).
- Boccaletti et al. (2006) S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. U. Hwang, Physics Reports 424, 175 (2006).
- Nishikawa et al. (2003) T. Nishikawa, A. Motter, Y. Lai, and F. Hoppensteadt, Phys. Rev. Lett. 91, 014101 (2003).
- Motter et al. (2005) A. Motter, C. Zhou, and J.Kurths, Phys. Rev. E 71, 016116 (2005).
- Moreno and Pacheco (2004) Y. Moreno and A. F. Pacheco, Europhys. Lett. 68 (2004).
- Oh et al. (2005) E. Oh, K. Rho, H. Hong, and B. Kahng, Phys. Rev. E 72, 047101 (2005).
- Lee (2005) D.-S. Lee, Phys. Rev. E 72 (2005).
- Yook and Meyer-Ortmanns (2005) S.-H. Yook and H. Meyer-Ortmanns, Physica A p. 781 (2005).
- Zhou and Kurths (2006) C. Zhou and J. Kurths, Phys. Rev. Lett. 96, 164102 (2006).
- Wang (2002) X. Wang, International Journal of Bifurcation and Chaos 12, 885 (2002).
- Wang and Chen (2002) X. F. Wang and G. Chen, IEEE Transactions on Circuits and Systems-I 49, 54 (2002).
- Parlitz (1996) U. Parlitz, Phys. Rev. Lett. 76, 1232 (1996).
- Maybhate and Amritkar (2000) A. Maybhate and R. E. Amritkar, Phys. Rev. E 61, 6461 (2000).
- Hegazi et al. (2002) A. S. Hegazi, H. N. Agiza, and M. M. El-Dessoky, Int. J. Bif. Chaos 12, 1579 (2002).
- Chen et al. (2004) S. H. Chen, J. Hu, C. Wang, and J. Lu, Phys. Lett. A 321, 50 (2004).
- Lu and Cao (2005) J. Lu and J. Cao, Chaos 15, 043901 (2005).
- Atay (2003) F. M. Atay, Phys. Rev. Lett. 91, 094101 (2003).
- Atay et al. (2004) F. M. Atay, J. Jost, and A. Wende, Phys. Rev. Lett 92, 144101 (2004).
- Sales-Prado et al. (2007) M. Sales-Prado, R. Guimerá, A. A. Moreira, and L. A. N. Amaral, PNAS 104, 15224 (2007).
- Jin et al. (1999) C. W. Jin, I. Marsden, X. Q. Chen, and X. B. Liao, J. Mol. Biol. 289, 683 (1999).
- Goto et al. (1998) S. Goto, T. Nishioka, and M. Kanehisa, Bioinformat. 14, 591 (1998).
- Zhou and an Lu (2007) J. Zhou and J. an Lu, Physica A 386, 481 (2007).
- Sorrentino and Ott (2008) F. Sorrentino and E. Ott, Phys. Rev. Lett. 100, 114101 (2008).
- Yu and Cao (2006) W. Yu and J. Cao, Chaos 16 (2006).
- Wu (2008) X. Q. Wu, Physica A 387 (2008).