A hierarchical heteroclinic network: Controlling the time evolution along its paths
EEPJ manuscript No. (will be inserted by the editor)
A hierarchical heteroclinic network
Controlling the time evolution along its paths
Maximilian Voit a and Hildegard Meyer-Ortmanns b Physics and Earth Sciences, Jacobs University Bremen, P.O. Box 750561, 28725 Bremen,Germany
Abstract.
We consider a heteroclinic network in the framework of win-nerless competition of species. It consists of two levels of heterocliniccycles. On the lower level, the heteroclinic cycle connects three saddles,each representing the survival of a single species; on the higher level, thecycle connects three such heteroclinic cycles, in which nine species areinvolved. We show how to tune the predation rates in order to generatethe long time scales on the higher level from the shorter time scales onthe lower level. Moreover, when we tune a single bifurcation parameter,first the motion along the lower and next along the higher-level het-eroclinic cycles are replaced by a heteroclinic cycle between 3-speciescoexistence-fixed points and by a 9-species coexistence-fixed point, re-spectively. We also observe a similar impact of additive noise. Beyondits usual role of preventing the slowing-down of heteroclinic trajectoriesat small noise level, its increasing strength can replace the lower-levelheteroclinic cycle by 3-species coexistence fixed-points, connected byan effective limit cycle, and for even stronger noise the trajectoriesconverge to the 9-species coexistence-fixed point. The model has appli-cations to systems in which slow oscillations modulate fast oscillationswith sudden transitions between the temporary winners.
A heteroclinic network is a sequence of trajectories connecting saddles σ , ..., σ n in atopological network. In particular a heteroclinic cycle is an invariant set consistingof the union of a set of saddles and the corresponding trajectories, which are back-ward asymptotic to each saddle σ i and forward asymptotic to σ i +1 [1]. Although thisparticular nonlinear dynamics looks exceptional, it is frequently found in ordinarydifferential equations under certain constraints like symmetries [2] or delay [3]. Ac-cordingly it is predicted in models of coupled phase oscillators [3,4], vector models[3], pulse-coupled oscillators [5] and models of winnerless competition [6]. The appli-cations range from social systems [7,8] to ecological systems like foodwebs [6,9], fluidmechanics [10], chemostats [11], computation [5] to neuronal networks [12,13,15,16,17,18,19,20]. In relation to neuronal systems, heteroclinic cycles (and, more generally,heteroclinic sequences) in models of winnerless competition provide mechanisms for a e-mail: [email protected] b e-mail: [email protected] a r X i v : . [ n li n . AO ] J un Will be inserted by the editor generating transient dynamics, which can be sensitive to the very input and robustagainst perturbations at the same time [21]. The transient feature is appreciated andfavored as compared to infinite-time limits, as realized in attractors like simple fixedpoints or limit cycles. The reason is that cognitive processes in the brain and episodesin ecological or social systems are inherently transient themselves.More specifically, for neuronal networks heteroclinic dynamics was proposed asa mechanism of binding between different information modalities in the brain [12].Here the winnerless competition via inhibitory connections is between active brainmodes, representing the temporal processing of different kind of input. The hetero-clinic network stands for multi-dimensional binding options, without implementingany hierarchy. An explicit hierarchy in time scales is implemented in the so-calledchunking dynamics of the brain [13]. Chunking refers to the phenomenon that thebrain uses to perform information processing of long sequences by splitting them intoshorter information items, called chunks. In [13] the generalized Lotka-Volterra equa-tions only serve as an elementary building block for different levels of the chunkinghierarchy, describing the competition between informational items to produce sta-ble sequences of metastable states. The very generation of chunks and the controlof the performance of tasks such as a temporally ordered sequence are described byadditional equations in [13].In contrast to this description, the hierarchy in time scales, which we considerhere, is merely implemented via the choice of rates in the generalized Lotka-Volterraequations (GLV) on the price that the choice of rates has to be tuned towards someintervals and in a specific order of sizes. From the physics point of view our maininterest is in a possible modulation of fast oscillations by slow oscillations as it is acommonly observed feature in brain dynamics. In particular we are interested in howa structural hierarchy in the attractor space (heteroclinic cycles of heteroclinic cycles)can induce such a dynamical generation of time scales. On the two levels of hierarchywhich we consider in this paper, the elementary items correspond to the dominanceof single species in 1-species saddles. The short time scales refer to fast oscillationsbetween different saddles within one small heteroclinic cycle (SHC) (analogous to one“chunk”), while the slow time scale is generated by a large heteroclinic cycle (LHC)between three SHCs. Our GLV-dynamics is less rich than the chunking dynamicsas considered in [13], nevertheless it already allows for a tuning of time scales overorders of magnitude. So the goal is then to control the path of the species trajectoriesthrough a desired heteroclinic network.From the structural point of view a remarkable property of this system is thatthe invariant sets that are the vertices of the LHC are themselves heteroclinic cycles.Such constructions have been considered before in [14] in the context of depth-twoheteroclinic networks. Our work extends the study of [14] by explicitly providing thevery construction of such a network, in view of illustrating the option of dynamicallygenerating different time scales and analyzing the role of noise.The paper is organized as follows. In section 2 we present the model with themain focus on the construction of the predation matrix that ensures the desired LHCof SHCs. Section 3 gives the results on the dynamical generation of a hierarchy intime scales (section 3.1) and the reduction of hierarchy levels via tuning the deathrate (section 3.2) or increasing the noise strength (section 3.3). Section 4 provides thesummary and conclusions. ill be inserted by the editor 3
We study a system of generalized Lotka-Volterra equations, given by ∂ t s i = ρs i − γs i − (cid:88) j (cid:54) = i A i,j s i s j i ∈ { , ..., } , (1)where s i denotes the concentration of species i , ρ is the reproduction rate, γ thedeath rate, and A i,j the rate by which species j preys on species i . The set of A i,j constitutes the predation matrix A . Without loss of generality we fix the time scaleby setting the reproduction rate ρ = 1.Inspired by the predation matrix for rock-paper-scissors(RPS)-games, we choosethe predation matrix A as a block matrix. Its diagonal consists of 3 × m , similar to the RPS-predation matrix for which c = 1 and e = 0. Theoff-diagonal blocks m d have d on their diagonal and r for the remaining elements, m f is chosen accordingly. A = m m d m f m f m m d m d m f m , where (2) m = c ee cc e , m d = d r rr d rr r d , m f = f r rr f rr r f . Note the similarity of the block form of A to the form of m . In section 2.4 weshall see why we choose this block form and why we keep the entries in the blockmatrices as independent parameters c, e, d, f and r . Although this choice may lookrather peculiar, as it stands it still allows a large variety of transients and stationarybehavior. Our interest is in how much it is possible to control the heteroclinic dynamics by a suit-able choice of predation rates c, e, d, f and r so that the vector of nine species evolvesbetween 1-species saddles along a prescribed trajectory in the nine-dimensional phasespace, as indicated in Figure 1. The vertices of the heteroclinic network are 1-speciessaddles, at which a single species has a non-vanishing concentration, being the tem-porary “winner of the game”. The links of the network are heteroclinic connections.The rates shall be chosen in such a way that each saddle has two incoming contractingdirections and two outgoing repulsive directions of different strength towards whichthe trajectory can escape. These directions are coming from and going to two othersaddles, respectively, while the remaining five directions in phase space should onlystabilize the desired trajectories, while moving along one of the indicated heteroclinicconnections.As indicated in Figure 1, we search for small heteroclinic cycles (SHCs) betweenspecies 1 , ,
3, as well as 4 , , , ,
9. The species within the SHCs play thenRPS, for which species 1 preys on 3, 3 on 2, and 2 on 1 etc.. (Note that the arrows,indicated along the heteroclinic cycles and following the time evolution of a desiredtrajectory, are reverse to the predation relations of who preys on whom.)
Will be inserted by the editor
In addition, we want to have heteroclinic connections also between saddles be-longing to different clusters, for example, from saddle σ to σ , and σ to σ andback to σ , making up a large heteroclinic cycle (LHC). As indicated in Figure 1,the directed connections between different SHCs are allowed from each saddle via asingle heteroclinic connection. This means that an LHC amounts to a cycle betweenSHCs, where the connections between the different SHCs may be realized via any ofthe saddles from an SHC. The distinction between small and large heteroclinic cyclesshould not be understood as referring to the distances along heteroclinic connectionsin phase space, which are all of the same order of magnitude due to the choice ofsaddle coordinates. In Figure 1 they appear different only due to the projections ontwo-dimensional space. Instead, it refers to short and long time scales, as we shallchoose the rates in a way that typically the trajectory spends a long time withinone of the SHCs, before it escapes to the next SHC and altogether repeatedly orbitsalong LHCs. If we define the long time scale as one complete revolution of an LHC(“complete” in the sense that each of the three SHCs has been visited once), thistime is then roughly three times the number of revolutions within one SHC, as thetransition between two SHCs happens almost instantaneously.Due to the characteristic slowing-down effect of heteroclinic cycles (as long asno noise is applied), it is not meaningful to talk about a characteristic period of anSHC or an LHC. Yet it makes sense to compare the speed with which the winnerswithin a single cluster alternate, with the speed with which whole clusters alternatein dominating the species of the other clusters if both are compared in the samefixed time interval. In this time interval one will then observe fast oscillating speciesconcentrations, modulated by slowly oscillating cluster dominance. This combinationof slow and fast oscillations is shared with chunking dynamics, here, however, to begenerated merely by the GLV-equations via a heteroclinic cycle of heteroclinic cycles.Note that while the system considered in [14] is similar in the overall structureto the present one, there each saddle connects to all other saddles of the successiveSHC rather than to one as in our case. The examples discussed in [14, sections 5 and6] are related to our system as they also possess a Z × Z symmetry and a similartopology. The Z symmetry of the cubic equations in [14] is absent in our system, asare certain degeneracy breaking terms.In the following we shall see how we can control an actual trajectory along thehypothetical network of Figure 1. For the calculation of fixed points, note that either s i itself must be zero in equa-tion (1), or the bracket term that remains when s i is factored out. Thus, there aremaximally 2 = 512 fixed points in total, which can be partitioned into classes ac-cording to their permutation symmetry. In Table 1 we list the classes which arerelevant for the further considerations, together with their general form, coordinatesand eigenvalues of the Jacobian, evaluated at these fixed points. F P stands for thefixed point, corresponding to global extinction of species. It is characterized by s i = 0for i ∈ { , ..., } , the nine eigenvalues are equal to ρ >
0, so that it is unstable.
F P stands for single-species fixed points. Their eigenvalues as indicated in thetable are in the main focus of our interest. They will be appropriately chosen toguarantee the desired stability properties.There are multiple classes of fixed points involving three species. Most relevantis the class containing only those species whose single-species fixed points belongto the same three-cycle, i.e. s i = s i +1 = s i +2 , s j = 0 ∀ j / ∈ { i, i + 1 , i + 2 } for i ∈{ , , } . These species form what we will call a “cluster”. Accordingly, we denote ill be inserted by the editor 5 σ σ σ σ σ σ σ σ σ Fig. 1.
Topology of a hierarchical heteroclinic network. The vertices are marked as blackdots. For the GLV-system of equation (1) they correspond to single-species fixed points σ i ,where species i has a population of s i = ργ . Links are heteroclinic connections, arrows marktheir direction. the class of the (local) three cluster-coexistence fixed points F P c . Each of them hasone pair of complex conjugated eigenvalues, corresponding to two unstable directionsfor parameters which satisfy equations (7) below. The remaining eigenvalues are − ρ , ρ ( c + γ − d + e − r ) c + γ + e , and ρ ( c + γ + e − f − r ) c + γ + e . The latter ones, both three times degenerate, arestable for c + γ + e > d + 2 r and c + γ + e > f + 2 r , respectively.Finally, there is the class F P g , containing only the global coexistence fixed point, s i = s ∀ i . Here, all species have the same concentration. It has four pairs of complexconjugated eigenvalues. For parameters according to equation (7) below, at least oneof them has a positive real part, making this fixed point unstable. The remainingeigenvalue is − ρ .While both of the latter fixed point classes are not directly involved in the hete-roclinic network, their stability properties affect the dynamics. All other fixed points(of which one up to seven coordinates are zero) are either not physical or irrelevantin the sense that the desired trajectory does not “care” about them due to the choiceof predation rates that we derive in the next section. Table 1.
Selected relevant fixed points of the nine-species generalized Lotka-Volterra system. class & name form / location eigenvaluesFP0: extinc-tion s i = 0 ∀ i ρ (9 × )FP1: singlespecies s i = ργ , s j = 0 ∀ j (cid:54) = i − ρ, ρ − cργ , ρ − dργ , ρ − eργ , ρ − fργ , ρ − ρrγ (4 × )FPc: clustercoexistence s i = s i +1 = s i +2 = ρc + γ + e , s j =0 ∀ j / ∈ { i, i + 1 , i + 2 } for i ∈ { , , } ρ (cid:16) ±√ √ − ( c − e )2+ c − γ + e (cid:17) c + γ + e ) , − ρ , ρ ( c + γ − d + e − r ) c + γ + e (3 × ), ρ ( c + γ − f + e − r ) c + γ + e (3 × )FPg: globalcoexistence s i = ρc + γ + d + e + f +4 r ∀ i ρ (cid:16) ±√ √ − ( c + d − e − f )2+ c + d + e + f − γ + r ) (cid:17) c + γ + d + e + f +4 r ) (4 × ), − ρ Will be inserted by the editor
The choice of the predation rates will be determined by the eigenvalues of the Ja-cobian, evaluated at the single-species fixed points. Let us start with the Jacobian,evaluated at the saddle σ as prototype for 1-species saddles. It is given by the matrix − ρ − eργ − cργ − fργ − rργ − rργ − dργ − rργ − rργ ρ − cργ ρ − eργ ρ − dργ ρ − rργ ρ − rργ ρ − fργ ρ − rργ
00 0 0 0 0 0 0 0 ρ − rργ . Its eigenvalues are listed in Table 1. It is, however, instructive to first write themin terms of general matrix elements A i,j in order to establish their desired properties,which can be read off from the heteroclinic network of Figure 1. From the perspectiveof vertex σ , the direction of the eigenvector to vertex σ should correspond to anexpanding direction within the small cycle SHC (es), towards σ a contracting direc-tion within the small cycle (cs), towards σ an expanding direction within the largecycle (el) and a contracting direction towards σ within the LHC (cl). We call thevector along the 1-axis the “radial”-direction, while the remaining four eigenvectordirections “transverse”.Before we formulate the conditions on the eigenvalues to achieve the desired mo-tion, we introduce a notation that is independent of the actual index values. Cor-responding to the former distinction between expanding and contracting directionswithin the small and large heteroclinic cycles and the remaining ones, we define func-tions es , el , cs , cl and r between indices, such that es ( i ) = j ( cs ( i ) = j ) with j theindex of the node, towards which the trajectory from σ i expands (from which themotion contracts towards σ i ), respectively, both within the small cycle. For exam-ple, es (1) = 2 and cs (1) = 3. Note that es ( i ) = i and es ( i ) = cs ( i ) for a 3-cycle.Similarly, el ( i ) and cl ( i ) yield the indices of the outgoing and incoming connectionsin the large cycle, while r ( i ) denote the indices assigned to i in the four transversedirections ( r ( i ) := { , . . . , } \ { i, es ( i ) , el ( i ) , cs ( i ) , cl ( i ) } ).Let us start with the expanding directions, that is, the outgoing heteroclinics from σ i to σ es ( i ) and σ el ( i ) . Both directions must be unstable, their eigenvalues positive. Inaddition we require that the small cycle is preferred over the large one in the sensethat the trajectory spends some time within the SHC, before it escapes to the LHC.The inequality between the corresponding eigenvalues translates into0 < A es ( i ) ,i < A ec ( i ) ,i < γ (3)as the first condition. The coefficient of e i in the corresponding eigenvector − A i,j γ − A j,i e i + e j for j ∈ { es ( i ) , el ( i ) } is then negative, so that the species concentration s i decreaseswhen following these directions.Next we turn to the contracting directions, associated with the incoming hete-roclinics that reach σ i from σ cs ( i ) and σ cl ( i ) . As they are stable, the correspondingeigenvalues must be negative. In addition, we impose a condition that is based on aconjecture made in [22]. Accordingly, the absolute values of eigenvalues of contracting ill be inserted by the editor 7 directions should be larger than that of the strongest expanding direction to achieveasymptotic stability in the heteroclinic network. This means (cid:12)(cid:12)(cid:12)(cid:12) ρ γ − A cs ( i ) ,i γ (cid:12)(cid:12)(cid:12)(cid:12) > (cid:12)(cid:12)(cid:12)(cid:12) ρ γ − A es ( i ) ,i γ (cid:12)(cid:12)(cid:12)(cid:12) ∧ (cid:12)(cid:12)(cid:12)(cid:12) ρ γ − A cl ( i ) ,i γ (cid:12)(cid:12)(cid:12)(cid:12) > (cid:12)(cid:12)(cid:12)(cid:12) ρ γ − A es ( i ) ,i γ (cid:12)(cid:12)(cid:12)(cid:12) , (4)which results in the conditions A cs ( i ) ,i > γ − A es ( i ) ,i ∧ A cl ( i ) ,i > γ − A es ( i ) ,i . (5)Concerning the transverse directions, in [22] it is conjectured that a sufficient(but not necessary) condition for asymptotic stability of the embedding heteroclinicnetwork is that in particular all transverse eigenvalues are negative. In our systemthis leads to A j,i > γ ∀ j ∈ r ( i ) . (6)At this point it is then suggested to introduce as many different parameters in thepredation matrix as are needed for satisfying conditions (3), (5), and (6). We definethe parameters e = A es ( i ) ,i , f = A el ( i ) ,i , c = A cs ( i ) ,i , d = A cl ( i ) ,i , and r = A j,i forall i , with j ∈ r ( i ) . Thereby we guarantee permutation symmetry both within thesmall cycles and the large cycles, i.e., between the three-cycles among each other.This finally yields the predation matrix equation (2). In terms of these parameters,conditions (3), (5), and (6) take the form0 < e < f < γ ∧ c > γ − e ∧ d > γ − e ∧ r > γ . (7)Note that this uniform choice of parameters (both c, d, e, f , and r within the predationmatrix and ρ and γ in (1)) introduces a Z × Z symmetry. This may look ratherartificial in view of any applications, where such a finetuning of parameters is highlyunlikely. We verified that the original dynamics near a heteroclinic cycle of heterocliniccycles persists when breaking this symmetry, see the next section.Table 2 summarizes the eigenvalues of the Jacobian at the single-species fixed point σ . The function ν ( x ) is defined as ν ( x ) = ρ − ρxγ . Numerical values are rounded to twodecimals and evaluated for our standard choice of parameters: ρ = 1 , γ = 1 . , r =1 . , e = 0 . , f = 0 . , c = d = 2. As abbreviations for the directions we use “rad”=radial, “exp.” = expanding, “ctr.” = contracting, ”l” = large, ”s” = small, ”trns.” =transverse. Table 2.
Table of eigenvalues at the single species fixed point σ . For further explanationssee the text. direction rad exp.s ctr.s exp.l trns. trns. ctr.l trns. trns.general − ρ ν ( A , ) ν ( A , ) ν ( A , ) ν ( A , ) ν ( A , ) ν ( A , ) ν ( A , ) ν ( A , )parameters − ρ ν ( e ) ν ( c ) ν ( f ) ν ( r ) ν ( r ) ν ( d ) ν ( r ) ν ( r )numeric − . − .
87 0 . − . − . − . − . − . In order to modulate fast oscillations by slow oscillations, we are primarily interestedin heteroclinic connections as indicated by the red (long-dashed) trajectory in Figure
Will be inserted by the editor , ,
6, it escapes to the LHC at σ and later returns to σ , closing the LHC afterspending some time within the SHCs. While Figure 2 a) is schematic, Figure 2 b)shows a projection of a simulated LHC of SHCs on a three-dimensional hyperplane ofphase space, where the color codes time. Initially the system starts from a randomlychosen location and performs a first LHC (black). At later times, blue, red, yellowcolors overwrite the black one in repeated orbits along both cycles. The projectionis chosen so as to represent the hierarchy similar to the sketch. Note that the size ofthe cycles in phase space does not differ as suggested by the Figure. All saddles areequidistant from the origin, so that inter- and intra-cluster distances between saddlesare equal in phase space. σ σ σ σ σ σ σ σ σ a) s + s + s + . ( s + s + s ) b) s + s + s + . ( s + s + s ) s + s + s + . ( s + s + s ) s + s + s + . ( s + s + s )
0 1000 2000time
Fig. 2.
Sketch of trajectories in the heteroclinic network with two sample trajectories. a)schematic: Solid grey and black lines mark heteroclinic connections between the saddles σ i .The red curve is a sketch of a trajectory performing LHCs of SHCs on both levels of thehierarchy. The blue dashed curve connects cluster-coexistence fixed points (dotted), when thelowest hierarchy level is broken as discussed in section 3.3. b) Plot of a simulated trajectory(cf. fig. 3) in a projection of phase space. The color codes time, from black (early) to late(yellow). For further explanations see the text. Figure 3 shows the actual time evolution of the nine species according to equation(1) with oscillating concentrations. Panel a) displays the typical time characteristicsof heteroclinic cycles, both for the LHC and the SHCs. Clusters performing the LHCare colored in blue red, green, and species, performing SHCs within one cluster, indifferent shades of the corresponding cluster color. We see an increase of time intervalsboth between different shades of one color and between colors. The first one reflectsthe slowing-down within an SHC, the latter one within an LHC. So, as time goes on,the system spends more and more time in the vicinity of an SHC, for which the speciesof a certain cluster are dominant, and within that SHC, more and more time close tosingle-species saddles. So far the simulations are performed without noise. Moreover,it is evident that the switching between saddles and between clusters occurs almostinstantaneously. The similar system in [14, Fig. 9] shows comparable dynamics. Panelb) shows the same data with a logarithmic s i axis to reveal the activities also ofspecies, which do not belong to the dominant cluster (the red one for t > s , ..., s ) are strongly suppressed during the dominanceof s , ..., s , they go on chasing each other rather than being completely frozen. Thisis in contrast to chunking dynamics as considered in [13], where the species get frozen. ill be inserted by the editor 9 The same qualitative dynamics as in Figure 3 is also observed if we break the Z × Z symmetry in the following ways: We replace the predation matrix A → ˜ A = A ◦ (1+ bW ) where ◦ is the entrywise product, b the “strength” of variation and W ∈ R × a random matrix (elements chosen from a Gaussian distribution with zero mean andunit variance). In a similar way we modify the reproduction and death rates for eachspecies i individually, defining: ρ i = ρ · (1 + bu i ) and γ i = γ · (1 + bv i ), where u, v ∈ R are random vectors (chosen from a Gaussian distribution as above). For both ways ofsymmetry breaking (either both individually or combined) the original dynamics neara heteroclinic cycle of heteroclinic cycles persists as long as the variation b is smallenough. Explorative investigations yield for b = 0 .
01 almost unaffected dynamics, andnoticeable quantitative changes for b = 0 .
1. Strong variations (e.g. b = 0 .
2) lead tosevere changes of the dynamics when the conditions (7) are violated and bifurcationsoccur. s i time a) s s s s s s s s s SHCLHC s i time b) Fig. 3.
Time evolution according to equation (1), as a) linear plot and b) log-linear plot. a)The times during which species 1 , , ρ = 1 , γ = 1 . , r = 1 . , e = 0 . , f = 0 . , and c = d = 2. Initial conditions are sampled uniformly randomly from [0 , . The time scales of slow oscillations between clusters and short oscillations betweensingle-species saddles can differ by orders of magnitude. In spite of the slowing-downeffect without noise, we can compare the time scale say of the first performance of aLHC in Figure 3 a) that takes of the order of 3000 time units, while the latest SHCof the red cluster takes some hundred time units.To further quantify the preference of SHCs over LHCs, we plot in Figure 4 ahistogram of transitions between all saddles, for a given fixed choice of parameters.The histogram should be read from columns to rows. So most transitions occur withinthe small cycles (e.g. 1 → →
4, 9 → σ i as being reached when the corresponding species concentration s i exceedsa threshold chosen as half of the concentration the fixed point, where it is located.Transitions between three-cycles happen with equal frequency to all saddles of thethree-cycles (e.g. 1 → →
5. Here, possibly a transition from 1 to 2 was already in progress,but not detected, so that the actual sequence would have been 1 → → Fig. 4.
Histogram of transitions between the saddles of the nine species GLV system (tobe read from row to column). Most transitions occur within the small cycles (e.g. 1 → → ρ = 1 , γ = 1 . , r = 1 . , e = 0 . , f = 0 . , and c = d = 2. The statistics is collected over 1000 runs of 1000 time steps each. Next we vary the parameters e and f , as it is particularly the choice of these ratesthat determines the dwell times within an SHC as compared to those at individualsaddles and therefore allows a tuning of the preference of SHCs over LHCs. Figure 5shows the ratio R of the number of transitions within small cycles to transitions toother clusters within large cycles, where f was changed, while e was kept fixed. R (f-e) Fig. 5.
Ratio R of the number of small over large heteroclinic cycles. R is plotted as afunction of the difference of parameters f and e . Simulations are sampled over 10000 runs of1000 time steps per data point. The other parameters that were used are ρ = 1 , γ = 1 . , r =1 . , e = 0 . , and c = d = 2. The reason why we have to collect a statistics in spite of the deterministic dynam-ics is due to the random choice of initial conditions, uniformly chosen from [0 , . If we tune the death rate γ , the system undergoes a sequence of Hopf bifurcations,whose order depends on the concrete choice of parameter values. We distinguish twotypes of Hopf bifurcations. One of them occurs simultaneously at the different cluster ill be inserted by the editor 11 coexistence fixed points at γ ( c ) , where the real part in the pair of complex conjugateeigenvalues of the cluster-coexistence fixed points changes sign, so that these fixedpoints become attracting. As a consequence, the lowest level of the heteroclinic cyclegets extinct, and the dynamics follows a heteroclinic cycle between three 3-speciescoexistence fixed points F P c , as displayed in Figure 6 b) and the blue trajectory inthe schematic plot of Figure 2 a).The other type of Hopf bifurcations happen all at the global coexistence fixed point F P g . Each of them stabilizes eigen-directions, corresponding to one pair of complexconjugated eigenvalues of F P g . In the order of increasing γ the last one occurs at γ ( g )4 = max { ( c + d + e + f − r ) , ( c − d + e − f + 4 r ) , ( − c + d − e + f + 4 r ) } . Asa result of this last bifurcation, the global coexistence fixed point becomes attractingand the system dynamics converges to this fixed point, as shown in Figure 6 c). Thiseliminates the second hierarchy level (slow time scale) from the system. s i time a) s i time b) s s s s s s s s s s i time c) Fig. 6.
Reduction of the hierarchy levels as a function of the death rate γ . Upon increasing γ , the system goes through a sequence of Hopf bifurcations. a) γ = 1 .
05 with both hierarchylevels still present; b) γ = 1 .
3, where the lowest level of the hierarchy is eliminated and thesystem is in a heteroclinic cycle between three-species coexistence fixed points. c) γ = 1 . ρ = 1 , r = 1 . , e = 0 . , f = 0 . , and c = d = 2. Initial conditions are sampleduniformly randomly from [0 , . The order of Hopf bifurcations depends on the concrete parameter choice. For ρ = 1, r = 1 . e = 0 . c = d = 2, we find0 < γ ( g )1 , = 1 . < γ ( c ) = 1 . < γ ( g )3 = 1 . < γ ( g )4 = 1 . . (8)This means, by tuning one parameter γ , we are able to reduce the hierarchy levels ofthe heteroclinic network from two levels in 6 a) (from γ ≥ . γ ≥ .
45 on) to none in c). Panels a) and b) correspond to Figures 9 and 10 in[14], respectively.
First we consider additive noise. When the dynamics of a systems follows a heterocliniccycle, a typical effect of noise is to prevent the slowing-down of the switching eventsbetween saddles. Instead, the switching between saddles continues with a period thatscales with the logarithm of the noise strength [3,23]. As a first result on the impactof noise we confirm this scaling behavior on both levels, the switching between specieswithin a cluster and between species belonging to different clusters. -12 -10 -8 -6 m ean d w e ll t i m e noise strength σ speciesclusterscl. broken Fig. 7.
Scaling of the mean dwell-times. Mean dwell-times scale logarithmically with thenoise strength σ . Plotted are the dwell-times of single species (blue) and of the three-cycles(clusters) (green). Dwell times are also shown for the case where the first hierarchy levelof SHC is broken and the species coalesce into 3-species coexistence fixed points (red). Thedwell-times within an SHC (green) and at a coexistence-fixed point are remarkably similarand drawn at separate noise values only for better readability. Parameters are ρ = 1 , γ =1 . , r = 1 . , e = 0 . , f = 0 . , and c = d = 2. For the broken hierarchy system γ = 1 . We study the influence of noise on equation (1) by including additive noise ξ i withstrength σ . Since species concentrations must be positive, we take its absolute value,so the dynamics is given by ∂ t s i = ρs i − γs i − (cid:88) j (cid:54) = i A i,j s i s j + σ | ξ i ( t ) | i ∈ { , ..., } , (9)where ξ i itself is Gaussian white noise with zero mean.The dynamics of equation (9) leads to periodic switching both on the lowesthierarchy level and the second level. To further quantify the switching events, wedefine as dwell-time the time the system stays in the neighborhood of a saddle. Here,a viable definition of the neighborhood is the range, where the species concentrationexceeds a threshold θ . For θ = 0 . S i = s i − + s i − + s i for i ∈ { , , } . The system is near cluster i when S i > θ (here we use θ = 0 . γ in the previous section, when we tune itsstrength to higher values, see Figure 6. For σ ≥ . · − noise leads to the extinctionof the lowest level of heteroclinic cycles; instead the noisy counterparts of the 3-speciescoexistence-fixed points get connected in an effective limit cycle, see Figure 8 b), whilefor σ ≥ . · − the effective limit cycle gets eliminated as well, and the dynamics isconstrained to the vicinity of the global coexistence-fixed point. The threshold valuesof noise, at which the crossover in its effect happens, (that is, from Figure 8 a) to b)and b) to c), respectively,) have been determined by measuring the variance in sumsof cluster concentrations in comparison to the sum of variances within single cluster ill be inserted by the editor 13 s i time a) s i time b) s s s s s s s s s s i time c) Fig. 8.
Reduction of hierarchy levels by means of an increasing noise strength σ . Uponincreasing noise from σ = 10 − in a) to σ = 10 − in b), the SHC is replaced by three3-species coexistence fixed points, connected by an effective limit cycle. Further increasingthe noise strength from σ = 10 − in b) to σ = 10 − in c) drives the system close to theglobal coexistence-fixed point. Parameters are ρ = 1 , γ = 1 . , r = 1 . , e = 0 . , f = 0 . , and c = d = 2. Initial conditions are sampled uniformly randomly from [0 , . concentrations. First the sum of variances should drastically drop, when the threespecies of the three clusters coalesce to 3-species coexistence-fixed points; next thevariance of the sum of different cluster concentrations should drop, when all speciescoalesce to the global coexistence-fixed point. This is precisely what we observed. Itshould be remarked that this effect of noise results from an implementation accordingto equation (9), that is, with effectively nonzero mean. If we only discard negativenoise contributions from unconstrained Gaussian white noise, once it leads to negativeconcentrations, an intermediate phase similar to Figure 8 b) is absent.Next we discuss the role of multiplicative noise as it would be naturally imple-mented in view of biological applications. If we replace the noise term in equation (9)by σs i ξ i ( t ), up to a certain strength its main effect is to smear out the single speciesfixed points, which can be understood as follows. The usual slowing down of thedynamics near a heteroclinic cycle is due to the trajectory coming closer and closerto the saddles, where the dynamics would get stuck. Noise increments parallel tothe unstable direction of the saddle prevent this slowing down. However, with mul-tiplicative noise such noise increments are extremely small, as the concentrations ofnon-dominant species are small. As a result, multiplicative noise neither prevents theslowing down nor does it lead to the collapse of hierarchy levels. Within the generalized Lotka-Volterra model with nine species we have demonstratedhow a suitable choice of predation rates can lead to a heteroclinic cycle of heterocliniccycles, each of which connects three saddles of single species dominance. The differentlevels of heteroclinic cycles can go along with a dynamical generation of a hierarchy intime scales, as the time it takes a full revolution on the upper level is determined bythe dwell-times within the lower-level cycles. In turn, this dwell time depends on thedwell-times in the vicinity of the individual saddles on the lowest level and the numberof revolutions within the small heteroclinic cycle. The time scales can mutually differby an order of magnitude. This feature may be of interest in view of hierarchical timescales observed in the brain, where slow oscillations modulate fast oscillations andexternal (sensory) input may select a pattern of “predation” rates in our effectivedescription of transient dynamics such that, for example, chunking is induced. Ourchoice of rates is specific but not singular as some explicit symmetry breaking in the predation matrix and some perturbation around the choice of reproduction and deathrates still support the heteroclinic dynamics.In view of ecological systems this kind of heteroclinic dynamics would mean thatmembers within a population compete with each other, while simultaneously also clus-ters of populations compete, resulting in particular in alternating, suddenly changingdominance of a whole cluster over other clusters of populations. It is well known thatrock-paper-scissors is played on many scales, starting from the genetic to the cellularand macroscopic level. To our knowledge it is an interesting and open question whenin real systems the same game on different scales can be traced back to one and thesame underlying microscopic dynamics (as in our artificial system). In ongoing workwe generalize the predation rules towards further levels of hierarchy in heterocliniccycles and explore the impact of time scales on spatial scales, when we assign theheteroclinic networks to a spatial grid and couple them.
Acknowledgments
We would like to thank Darka Labavi´c (University of Lille) for valuable discussions in thebeginning of this work. Also we thank an anonymous referee for drawing our attention toreference [14]. Financial support from the German Research Foundation (DFG, contractME-1332/28-1) is gratefully acknowledged.
Author contributions
H.M-O. and M.V. designed the model, discussed the results and contributed to the manuscript.M.V. carried out the numerical experiments and analytical calculations. H.M.-O. proposedthe project and wrote the main manuscript text in its final form. Both authors approved thefinal version.
References
1. C.M. Postlethwaite and J. H. P. Dawes, Nonlinearity , 1477 (2005).2. M. Krupa, J. Nonlinear Sci. , 129 (1997).3. H. Kori and Y. Kuramoto, Phys. Rev. E , 06214 (2001).4. J. Wordsworth and P. Ashwin, Phys. Rev. E , 066203 (2008).5. F. Schittler Neves and M. Timme, Phys. Rev. Lett. , 018701 (2012).6. V.S. Afraimovich, I. Tristan, R. Huerta, and M. I. Rabinovich, Chaos , 043103 (2008).7. C. Hauert and G. Szab´ o , Am. J. Phys. , 405 (2005).8. G. Szab´ o , J. Vukov, and A. Szolnoki, Phys. Rev. E , 047107 (2005).9. M. A. Nowak and K. Sigmund, Nature (London) , 138 (2002).10. Y. Tu and M. C. Cross, Phys. Rev. Lett. , 2515 (1992).11. S-B. Hsu and L-J. W. Roeger, J. Math. ANal. Appl. , 599 (2009).12. M. I, Rabinovich, V.S.Afraimovich, and P. Varona, Dyn. Systems (3), 433 (2010).13. M.I. Rabinovich, P. Varona, I. Tristan, and V.S. Afraimovich, Front.Comp.Neurosci. (22), 1 (2004).14. P. Ashwin and M. Field, Arch. Ration. Mech. Anal. , 107 (1999).15. M. I. Rabinovich, R. Huerta, and P. Varona, Phys. Rev. Lett. , 014101 (2006).16. M. I. Rabinovich, R. Huerta, P. Varona, and V. S. Afraimovich, PLoS Comput. Biol. ,e1000072 (2008).17. M. A. Komarov, G. V. Osipov, J. A. K. Suykens, and M. I. Rabinovich, Chaos ,015107 (2009).18. A. Szucs, R. Huerta, M. Rabinovich, and A. Selverston, Neuron , 439 (2009).19. M. A. Komarov, G. V. Osipov, and J. A. K. Suykens, Europhys. Lett. , 60006 (2009).20. M. A. Komarov, G. V. Osipov, and J. A. K. Suykens, Europhys. Lett. , 20006 (2010).ill be inserted by the editor 1521. V. S. Afraimovich, M. I. Rabinovich, and P. Varona, Int. J. Bifurcation Chaos (4),1195 (2004).22. P. Ashwin and C. M. Postlethwaite, Physica D , 26 (2013).23. D. Hansel, G. Mato, and C. Meunier, Phys. Rev. E48