Complex transitions to synchronization in delay-coupled networks of logistic maps
aa r X i v : . [ n li n . AO ] M a y Complex transitions to synchronization in delay-couplednetworks of logistic maps
Cristina Masoller ∗ Fatihcan M. Atay † Preprint.
Final version in
European Physical Journal
D 62:119–126, 2011.DOI: 10.1140/epjd/e2011-10370-7
Abstract
A network of delay-coupled logistic maps exhibits two different synchronization regimes,depending on the distribution of the coupling delay times. When the delays are homoge-neous throughout the network, the network synchronizes to a time-dependent state [Atay etal., Phys. Rev. Lett. , 144101 (2004)], which may be periodic or chaotic depending onthe delay; when the delays are sufficiently heterogeneous, the synchronization proceeds to asteady-state, which is unstable for the uncoupled map [Masoller and Marti, Phys. Rev. Lett. , 134102 (2005)]. Here we characterize the transition from time-dependent to steady-statesynchronization as the width of the delay distribution increases. We also compare the two tran-sitions to synchronization as the coupling strength increases. We use transition probabilitiescalculated via symbolic analysis and ordinal patterns. We find that, as the coupling strengthincreases, before the onset of steady-state synchronization the network splits into two clusterswhich are in anti-phase relation with each other. On the other hand, with increasing delayheterogeneity, no cluster formation is seen at the onset of steady-state synchronization; how-ever, a rather complex unsynchronized state is detected, revealed by a diversity of transitionprobabilities in the network nodes. PACS:
A fascinating and intriguing feature of spatially extended systems composed of many interact-ing units, like chanting crowds, tropical Malaysian flashing fireflies, pacemaker heart cells, cellsgoverning the circadian rhythms, pedestrians crossing the Millennium Bridge, etc., is that theycan synchronize even when the units are spread over wide spatial areas [1, 2, 3]. In order tounderstand their synchronization phenomena, these systems have been modeled by networksof coupled phase oscillators, like the Kuramoto model [4], and by networks of coupled maps[5, 6, 7], such as circle maps [8, 9], Bernoulli maps [10], logistic maps [11, 12, 13], Rulkov maps[14, 15, 16] etc.In systems of coupled units, communication delays naturally arise from a realistic consid-eration of the finite speed of information transmission between pairs of units, and can have agreat impact on their collective behavior. In particular, in networks of coupled maps, synchro-nization phenomena in the presence of time-delays has received considerable attention and isstill an active research area [10, 11, 12, 13, 14, 15, 17, 18]. Networks of delayed-coupled mapsare popular for studying the effects of delayed interactions because one can simulate large en-sembles of coupled units, even in the presence of heterogeneous and long delays, with a greatreduction of computational time and memory requirements, as compared to delay-differential ∗ Departament de F´ısica i Enginyeria Nuclear, Universitat Polit`ecnica de Catalunya, Colom 11, Terrassa 08222,Barcelona, Spain. [email protected] † Max Planck Institute for Mathematics in the Sciences, 04103 Leipzig, Germany. [email protected] ate-equations. The logistic map has been a popular choice because is a prototype exampleof how chaotic dynamics and universal scaling laws [19, 20, 21, 22] arise in simple non-linearsystems.In networks of delayed-coupled logistic maps, when the delays are heterogeneous the net-work exhibits a synchronized collective behavior that is qualitatively different from that ofinstantaneously interacting units, or by units interacting with homogeneous delays [23, 24].Heterogeneous delays can enhance the synchronizability of the network, but they can also affectits synchronized dynamics. A network of delayed coupled logistic maps displays two qualita-tively different synchronization regimes, depending on the delay distribution. When the delaysare homogeneous throughout the network, the network synchronizes to a time-varying state[23], and the synchronizability depends mainly on the network architecture; when the delaysare sufficiently heterogeneous, the network synchronizes to a steady-state, which is unstablefor the uncoupled maps [24], and the synchronizability depends mainly on the average numberof neighbors per node.The stability of the steady-state of delay-coupled maps is well-understood when the delay ishomogeneous (delta-distributed): Ref. [25] gave exact conditions for stability and showed thatthe largest eigenvalue of the Laplacian matrix determines the effect of the network structure onstability. Such precise results are unavailable for arbitrary delay distributions. Nevertheless, itis known that distributed delays can induce or improve stability of the steady-state in coupledlimit-cycle oscillators [17], or in more general delay-differential equations in the vicinity of aHopf instability [26]. A recent example is reported in [27], for an integro-differential equationdescribing the collective dynamics of a neural network with distributed signal delays: WithGamma distributed delays, which are less dispersed than the exponential distribution, thesystem exhibits reentrant phenomena (i.e., the stability is lost but then recovered as the meandelay is increased), while with delays that are more highly dispersed than exponential, thesystem does not destabilize.The aim of this paper is to characterize the transition to the two synchronized regimesof delayed coupled logistic maps (time-dependent for homogeneous delays and steady-statefor heterogeneous ones) as the coupling strength or as the width of the delay distributionincreases. The degree of synchronization is measured in terms of the transition probabilitiesin the network nodes, which are calculated via symbolic analysis and ordinal patterns. Thesymbolic method is based in dividing the state space of a given node into two regions andconsidering the relative frequencies of the transitions between those regions [28]; the ordinalpatterns method is based in defining patterns in the time-series of a node that result fromordering relations in consecutive values in the series [29], and computing the relative frequenciesof the transitions between those patterns. The paper is organized as follows: Section II presentsthe network model and the magnitudes employed to quantify the degree of synchronization.Section III presents the results, and Sec. IV contains a summary and the conclusions. We consider N logistic maps coupled as x i ( t + 1) = f [ x i ( t )] + ǫk i N X j =1 w ij ( f [ x j ( t − τ ij )] − f [ x i ( t )]) , (1)where t is a discrete time index, i is a discrete spatial index, f ( x ) = ax (1 − x ) is the logisticmap with parameter a , ǫ is the coupling strength, τ ij denotes the delay in the link from node j to i , w ij are the elements of the adjacency matrix w whose values equals 1 if there is a linkfrom node j to node i and 0 otherwise, and k i is the in-degree of the node i , k i = P j w ij .Here, τ and w are not restricted to be symmetric matrices.When the delays are sufficiently heterogeneous, the solution in the spatially homogeneoussteady-state, x i ( t ) = x ∀ i, (2)is stable in a certain range of coupling strengths [24], where x is the fixed point of theuncoupled logistic map, x = f ( x ) = 1 − /a. (3) e will refer to this solution as “steady-state synchronization”. In contrast, when the delaysare homogeneous throughout the network ( τ ij = τ ∀ i, j ) the network synchronizes to a time-dependent state [23], x i ( t ) = x ( t ) ∀ i, (4)where x ( t ) is a solution of x ( t + 1) = f [ x ( t )] + ǫ ( f [ x ( t − τ )] − f [ x ( t )]) , (5)and thus, the dynamics can be periodic or chaotic depending on τ . We will refer to thissituation as “time-dependent” synchronization.Clearly, other “out of phase” synchronization regimes, where the different nodes maintaincertain lag-times among them, are also possible. For example, a 1D linear globally-couplednetwork with distance-dependent delays, τ ij = | i − j | /v , where v is the speed of informationtransmission, synchronizes to a state in which the nodes evolve along a periodic orbit of theuncoupled logistic map (i.e., x i ( t ) is a solution of x i ( t + 1) = f [ x i ( t )]), while the spatialcorrelation of the nodes along the network is such that x i ( t ) = x j ( t − τ ij ) ∀ i, j (i.e., each map“sees” all other maps in his present, current, state) [30, 31]. In the following we only focus on“steady-state” and “time-dependent” synchronization.To capture the degree of synchronization and to distinguish between steady-state syn-chronization, Eq. (2), and time-dependent synchronization, Eq. (4), we use the followingmeasures:1) The variance of the nodes’ states, σ = 1 N h N X i =1 ( x i ( t ) − h x i s ) i t (6)where h . i s denotes an average over the nodes of the network, and h . i t denotes an average overtime.2) The variance of the distance to the steady state, σ ′ = 1 N h N X i =1 ( x i ( t ) − x ) i t , (7)where x is the fixed point of the uncoupled logistic map, Eq. (3).One can notice that σ = 0 if and only if x i = x j = h x i s ∀ i, j , while σ ′ = 0 if andonly if x i = x ∀ i . Thus, σ ′ allows to distinguish synchronization in the steady state fromsynchronization in a time dependent state. In the former case, both σ and σ ′ are zero, inthe latter case, only σ is zero.We note that both σ and σ ′ are “global” indicators that give no information about themicroscopic local dynamics in the nodes of the network. To gain inside into this local dynamics,the transition probabilities in individual nodes can be computed via symbolic dynamics [28]or ordinal patterns [29], as follows.3) Transition probabilities computed via symbolic dynamics: At each node i , a two-symboldynamics is generated by the partition of the phase space as s i ( t ) = α if x i ( t ) ≤ x ∗ s i ( t ) = β otherwise, (8)where x ∗ is a threshold value, which in the following is chosen equal to the fixed point of theuncoupled logistic map, x . The transition probability in node i , P i,sd ( α, α ), is calculated as P i,sd ( α, α ) = P Lt =1 n ( s i ( t ) = α, s i ( t + 1) = α ) P Lt =1 n ( s i ( t ) = α ) , (9)where n is a count of the number of times of occurrence in a time-series of length L . The globalproperties of the network can be quantified by the variance of P i,sd ( α, α ) over the network[28], ζ sd = 1 N N X i =1 ( P i,sd ( α, α ) − h P sd ( α, α ) i s ) , (10) here h P sd ( α, α ) i s = (1 /N ) P Ni =1 P i,sd ( α, α ) is the average transition probability.4) In addition, in each node i , a sequence of symbols can be generated via a comparison ofconsecutive values (“ordinal patterns” of dimension two, as proposed by Brandt and Pompe[29]) s i ( t ) = α if x i ( t ) ≤ x i ( t + 1) s i ( t ) = β otherwise. (11)A nice advantage of this procedure is that it does not require the definition of a threshold. Asbefore, the transition probability in node i , P i,BP ( α, α ), can be calculated as P i,BP ( α, α ) = P Lt =1 n ( s i ( t ) = α, s i ( t + 1) = α ) P Lt =1 n ( s i ( t ) = α ) , (12)and its variance, ζ BP = 1 N N X i =1 ( P i,BP ( α, α ) − h P BP ( α, α ) i s ) , (13)where h P BP ( α, α ) i s = (1 /N ) P Ni =1 P i,BP ( α, α ), can be used to capture global properties of thenetwork. In the following we present the results for an Erd¨os-Renyi random network [32] of N nodeswith an average degree h k i s such that the network has a single component. Unless otherwiseexplicitly stated, N = 200, h k i s = 20 and the coupling delays are Gaussian distributed witha mean delay h τ i s = 5. The parameter that controls the delay heterogeneity is the standarddeviation of the delay distribution, normalized by the mean delay, c ∗ = σ τ / h τ i s . The param-eter of the logistic map is taken to be a = 4 and the simulations start with random initialconditions. Unless otherwise explicitly stated, the quantifiers σ , σ ′ , ζ sd and ζ BP are com-puted over time series of length L = 500, after the first 3000 iterations are disregarded, andthey are averaged over 20 stochastic trajectories, where the random initial conditions ( x i (0)),delay distribution ( τ ij ), and adjacency matrix ( w ij ) are varied.First we consider the transition to “steady-state” synchronization as the coupling strength ǫ increases, while the delay heterogeneity c ∗ is kept constant. The delays are sufficientlyheterogeneous such that, for large enough ǫ , the network synchronizes as x i = x ∀ i .Figure 1 displays σ , σ ′ , ζ sd and ζ BP vs. the coupling strength ǫ . It also displays the fourtransition probabilities, for one typical stochastic trajectory, in 20 randomly selected nodes,as computed via symbolic dynamics (circles) and ordinal patterns (squares). It can be seenthat before the onset of synchronization there is a formation of two clusters, as the transitionprobabilities P i ( α, α ), P i ( α, β ) and P i ( β, α ) are 0 in some nodes and 1 in others. One canalso notice that P i ( β, β ) is very small in all the nodes, and that the transition probabilitiescalculated with symbolic dynamics are very similar to those calculated with ordinal patterns.Further insight into the networks’ dynamics near the synchronization transition can ob-tained by examining the time evolution of the quantifiers, of the transition probabilities (nowcomputed over a moving time-window of length 500), and the dynamics of a few, randomlyselected nodes. These are shown in Fig. 2, where the coupling strength is slightly smallerthan that needed for ”steady-state” synchronization. In Fig. 2(d) the network configurationat a fixed time (i.e., a ’snapshot’ of the states of the nodes) is also shown. One can noticethat the nodes form two clusters that oscillate in anti-phase: when one cluster is above thefixed-point solution, the other one is below, and at the next time step, the two clusters switchtheir positions.Next, we consider the situation where the delay heterogeneity c ∗ increases, starting witha delay distribution that is a delta function ( c ∗ = 0), while the coupling strength ǫ is keptconstant. The coupling strength is strong enough that, for homogeneous delays, the networksynchronizes as x i = x j ∀ i, j (time-dependent synchronization), while for sufficiently heteroge-neous delays, the network synchronizes as x i = x ∀ i (steady-state synchronization). Figure 3 σ σ ´ α→α ζ α→β ζ β→α ζ −4 (i) β→β ζ ε α→α P i α→β P i β→α P i β→β P i ε Figure 1: Transition to “steady-state” synchronization for fixed delay heterogeneity ( c ∗ = 0 . ǫ . The quantifiers σ , σ ′ , and ζ are plotted vs. the couplingstrength. The transition probabilities P i in 20 randomly selected nodes are also shown. In pan-els (c)-(j) the transition probabilities are computed via symbolic dynamics (circles) and ordinalpatterns (squares; red online). 5 isplays the quantifiers vs. the delay heterogeneity and also displays the four transition prob-abilities for one typical stochastic trajectory, in 20 randomly selected nodes. In this scenario,for small delay heterogeneity the time-dependent synchronization is gradually lost, and as thedelay heterogeneity increases, there is a smooth transition to the steady-state synchronization.No cluster formation can be observed at the onset of “steady-state synchronization”, since thetransition probabilities are within a certain range of values.The dynamics of the network near “steady-state” synchronization is examined in Fig. 4,with parameters such that the heterogeneity of the delays is slightly smaller than that neededfor “steady-state” synchronization. Here one can notice that the nodes evolve together, ina single cluster, displaying slow oscillations around the steady state [compare the oscillationfrequencies in Figs. 2(f),(h),(j) with 4(f),(h),(j)]. The period and shape of these oscillationsvary with c ∗ . One should keep in mind that the scenario we are considering is with strongcoupling, such that, for c ∗ = 0 the array synchronizes in a time-dependent state; the networkdynamics near this state (with the presence of a small delay heterogeneity), is shown in Fig.5. The approach towards “steady-state synchronization” reveals ‘critical slowing down’ inthe sense that the amplitude of the oscillations in Figs. 2(f),(h),(j) and 4(f),(h),(j) graduallydecreases with increasing ǫ or c ∗ , and there is a slow approach towards the fixed-point solution.The main differences being that for sufficiently heterogeneous delays and small coupling, thenetwork splits in two clusters which display fast anti-phase oscillations, while for large enoughcoupling but not sufficiently heterogeneous delays, the network approaches the fixed pointsolution as a single cluster and slow oscillations.One can then interpret the diversity of transition probabilities seen at the boundary ofsteady-state synchronization as “noise amplification”. When the network is almost or nearlysynchronized, for all the nodes we have x i ( t ) ∼ x and therefore very small variations near x result in a diversity of transition probabilities. This occurs when both ǫ or c ∗ is varied.However, because of the different way the network approaches the homogeneous solution,increasing ǫ yields two clusters and the transition probabilities are either close to 0 or to 1,while, increasing c ∗ yields a single cluster and the transition probabilities are within an intervalof values.Two-dimensional plots in the parameter space (coupling strength, delay heterogeneity),shown in Fig. 6, provide a more complete picture of the various dynamical regimes. We canrecognize two synchronization regions occurring for large coupling: steady-state synchroniza-tion for large delay heterogeneity [top-right corner in Figs. 6(a),(b), where both σ and σ ′ are zero], and time-dependent synchronization, for homogeneous delays [bottom-right cornerin Figs. 6(a),(b), where only σ is zero]. In addition, there is a narrow window of synchroniza-tion for weak coupling strength and almost homogeneous delays [ η ∼ . − .
2, bottom-leftcorner in Figs. 6(a),(b), where σ is zero and σ ′ is small]. This region was reported in [23]for homogeneous and odd delay values, and it can be seen from Fig. 6 that it is also robust tosmall delay heterogeneities.In Fig. 7 we consider finite-size and time-dependent effects during the onset of steady-statesynchronization, Figs. 7(a),(b) and of time-dependent synchronization, Figs. 7(c),(d). Weplot the time-evolution of the instantaneous values of σ and σ ′ [i.e., σ and σ ′ are computedas in Eqs. (6)–(7) but without time-averaging] for various network sizes N , while the averagenumber of neighbors per node is kept constant. Approaching the steady-state synchronization,there is a gradual decrease of the quantifiers, and initially their time-evolution is independent ofthe network size. In contrast, the approach to time-dependent synchronization, Figs. 7(c),(d)occurs abruptly, at a time that is nearly independent of the network size.For parameters close to “steady-state” synchronization critical slowing down occurs duringthe approach to the homogeneous steady state, as can be seen in Fig. 8, where we display thetime-variation of σ for various values of ǫ and c ∗ .We have checked the robustness of the above observations by considering delays that areexponentially distributed, and very similar results were found: the formation of two clustersbefore the onset of steady-state synchronization for increasing ǫ , while there is a single clusterfor increasing c ∗ . The small synchronization region that occurs for weak coupling strength isalso robust to exponentially distributed delays, as long as the width of the distribution is nottoo wide. The difference with Gaussian delays is that, with exponentially distributed delays,for strong coupling ( ǫ ≈ σ ′ is small and positive) but x 10 −3 (a) σ t σ ´ t0 2000 4000 6000 α→α P i t0 2000 4000 60000.80.91 (e) α→β P i t0 2000 4000 60000.951 (g) β→α P i t0 2000 4000 600000.020.04 (i) β→β P i t 0 100 200 x ( t ) , x ( t + ) (d) i x i ( t ) (f)t0 50 1000.60.81 x j ( t ) (h)t0 50 1000.60.81 x k ( t ) (j)t Figure 2: Time-variation of the quantifiers σ (a), σ ′ (b), and of the transition probabilities (c),(e), (g), (i), calculated in a moving time-windows of length L=600. The circles indicate transitionprobabilities computed via symbolic dynamics; the squares (red online), via ordinal patterns. Theparameters are such that the coupling strength is slightly below that required for synchronizationin the steady-state ( ǫ = 0 . c ∗ = 0 . x of the map f .7 −3 (a) σ σ ´ −3 (c) α→α ζ −3 (e) α→β ζ −3 (g) β→α ζ −3 (i) β→β ζ c ∗ α→α P i α→β P i β→α P i β→β P i c ∗ Figure 3: As Fig. 1, but keeping the coupling strength fixed ( ǫ = 0 .
9) and increasing the delayheterogeneity c ∗ . 8 x 10 −3 (a) σ t 0 2000 4000 6000 σ ´ t0 2000 4000 6000 α→α P i t0 2000 4000 60000.160.180.2 (e) α→β P i t0 2000 4000 60000.160.180.2 (g) β→α P i t0 2000 4000 60000.80.85 (i) β→β P i t x ( t ) , x ( t + ) (d) i0 50 1000.70.80.9 x i ( t ) (f)t0 50 1000.70.80.9 x j ( t ) (h)t0 50 1000.70.80.9 x k ( t ) (j)t Figure 4: As Fig. 2 but with the delay heterogeneity slightly below that required for synchronizationto the steady-state ( ǫ = 0 . c ∗ = 0 . −3 (a) σ t (b) σ ´ t0 2000 4000 6000 α→α P i t0 2000 4000 600000.51 (e) α→β P i t0 2000 4000 600000.51 (g) β→α P i t0 2000 4000 600000.51 (i) β→β P i t 0 100 200 x ( t ) , x ( t + ) (d) i0 50 10000.51 x i ( t ) (f)t0 50 10000.51 x j ( t ) (h)t0 50 10000.51 x k ( t ) (j) t Figure 5: As Fig. 4 but for small delay heterogeneity ( ǫ = 0 . c ∗ = 0 . c * (a) σ ε (b) σ ´ ε c * (c) ζ α → α ε (d) ζ α → α ε c * (e) ζ α → β ε (f) ζ α → β ε c * (g) ζ β → α ε (h) ζ β → α ε c * (i) ζ β → β ε (j) ζ β → β Figure 6: Synchronization quantifiers in the parameter space (coupling strength on the horizontalaxis, delay heterogeneity on the vertical). 11 −20 t σ (a) −20 t σ ´ (b) −20 t σ (c) −20 t σ (d)50 Figure 7: (a), (b) Time-evolution σ and σ ′ , during the onset of steady-state synchronization( ǫ = 0 . c ∗ = 0 . σ during the onset of time-dependent synchronizationfor homogeneous delays ( ǫ = 0 . c ∗ = 0). σ ′ remains finite and is not shown. (d) Time-evolutionof σ during the onset of time-dependent synchronization, in the window for weak coupling existingonly for homogeneous and odd delays ( ǫ = 0 . c ∗ = 0). σ ′ remains finite and is not shown. σ and σ ′ were computed for the various network sizes indicated in panel (a). −30 −20 −10 t σ (a) −30 −20 −10 t σ (b) c * =0.57c * =0.47c * =0.40c * =0.35 ε =0.9 ε =0.7 ε =0.67 ε =0.65 Figure 8: Time-evolution of σ during the transition to steady-state for (a) various values of ǫ , c ∗ = 0 .
57 (b) for ǫ = 0 . c ∗ . σ ′ exhibits similar behavior (not shown).12 he network remains synchronized, as σ = 0 and the transition probabilities in the nodes areall equal. To summarize, we have studied the transition to synchronization in a network of delay-coupledlogistic maps. When the coupling delays are homogeneous throughout the network, the net-work synchronizes to a time-dependent state; when the delays are sufficiently heterogeneous,the synchronization occurs in a steady-state. We employed global and local measures to char-acterize the synchronization transitions. The global measures are the standard deviation ofthe distance to the synchronized state, as well as the standard deviation of the transition prob-abilities in the nodes. The transition probabilities were computed using symbolic analysis andordinal patterns. We have found that, as the coupling strength increases or as the width of thedelay distribution grows, there is a gradual approach to the synchronized state, as seen withthe global indicators. An inspection of the local dynamics in the individual nodes, measuredby the transition probabilities, reveals that for increasing coupling there is the formation oftwo clusters before the steady-state synchronization, detected by the fact that the nodes ex-hibit two qualitatively different transition probabilities. For increasing delay heterogeneity, nocluster formation is seen at the onset of steady-state synchronization, but there is a diversityof values of transition probabilities.
This research was supported in part by the Spanish Ministerio de Educacion y Ciencia throughproject FIS2009-13360-C03-02, the Agencia de Gestio d’Ajuts Universitaris i de Recerca(AGAUR), Generalitat de Catalunya, through project 2009 SGR 1168, and the ICREA foun-dation. The authors acknowledge the hospitality of The Max Planck Institute for Physics ofComplex Systems, Dresden, where the initial steps of the research were taken.
References [1] L. Glass and M. C. Mackey,
From Clocks to Chaos: The Rhythms of Life , PrincetonUniversity Press, Princeton, NJ, 1988.[2] S. H. Strogatz, D. M. Abrams, A. McRobie, B. Eckhardt, and E. Ott, Nature , 43(2005).[3] A. Arenas, A. Diaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Physics Reports ,93-153 (2008).[4] J. A. Acebron, L. L. Bonilla, C. J. Perez Vicente, F. Ritort, and R. Spingleri, Rev. Mod.Phys.
137 (2005).[5]
Theory and Applications of Coupled Map Lattices , edited by K. Kaneko, Wiley, New York,1993.[6] R. E. Amritkar and S. Jalan, Physica A , 118-128(2003).[9] A. R. Sonawane, Phys. Rev. E
15] Q. Y. Wang, M. Perc, Z. S. Duan, and G. R. Chen, Phys. Rev. E , 094101 (2003).[18] F. M. Atay, S. Jalan and J. Jost Phys. Lett. A , 130–135 (2010).[19] M. J. Feigenbaum, J. Statist. Phys. , 25-52 (1978).[20] P. Coullet and C. Tresser, J. Phys. Colloque C , C5-25 (1978).[21] C. Tresser and P. Coullet, C. R. Acad. Sci. Paris , 577-580 (1978).[22] M. J. Feigenbaum, J. Statist. Phys. , 669-706 (1979).[23] F. M. Atay, J. Jost and A. Wende, Phys. Rev. Lett. , 144101 (2004).[24] C. Masoller and A. C. Marti, Phys. Rev. Lett. , 134102 (2005).[25] F. M. Atay and ¨O. Karabacak, SIAM J. Appl Dyn. Syst. , 508–527, (2006).[26] F. M. Atay, Disc. Cont. Dyn. Syst. S , 197–205 (2008).[27] T. Omi and S. Shinomoto, Phys. Rev. E , 033124 (2006).[29] C. Bandt and B. Pompe, Phys. Rev. Lett. , 056219 (2003).[31] C. Masoller, A. C. Marti and D. H. Zanette, Physica A , 186 (2003).[32] P. Erd¨os and A. R´enyi, Publ. Math. Inst. Hung. Acad. Sci. , 17-60 (1960)., 17-60 (1960).