Predicting the separation of time scales in a heteroclinic network

Abstract We consider a heteroclinic network in the framework of winnerless competition, realized by generalized Lotka-Volterra equations. By an appropriate choice of predation rates we impose a structural hierarchy so that the network consists of a heteroclinic cycle of three heteroclinic cycles which connect saddles on the basic level. As we have demonstrated in previous work, the structural hierarchy can induce a hierarchy in time scales such that slow oscillations modulate fast oscillations of species concentrations. Here we derive a Poincaré map to determine analytically the number of revolutions of the trajectory within one heteroclinic cycle on the basic level, before it switches to the heteroclinic connection on the second level. This provides an understanding of which parameters control the separation of time scales and determine the decisions of the trajectory at branching points of this network.


Introduction
Nonlinear dynamics of heteroclinic networks is frequently found in ordinary differential equations under certain constraints like symmetries [1] or delay [2]. It is predicted in models of coupled phase oscillators [2,3], vector models [2], pulse-coupled oscillators [4] and models of winnerless competition [5] that we consider in the following. Applications are manifold and range from social [6,7] and ecological [5,8] systems, to computation [4] and neuronal [9][10][11][12][13][14][15][16] networks. As it was emphasized in the work of V. Afraimovich and his coworkers [5,9,10,12,17], heteroclinic sequences in models of winnerless competition predict transient dynamics that shares features with cognitive dynamics: being simultaneously sensitive to the input and robust against perturbations. According to [17], heteroclinic dynamics may describe the binding between different information modalities in the brain (without an intrinsic hierarchy) [9], or chunking dynamics, which the brain uses to perform information processing of long sequences by splitting them into shorter information items [10]. There the hierarchy in time scales of chunking dynamics results from additional equations to the generalized Lotka-Volterra equations.
In our previous work [18] we considered a hierarchical heteroclinic network composed of a (large, superordinated) heteroclinic cycle (LHC) of three (small) heteroclinic cycles (SHCs). Each of these SHCs incorporates three saddles, corresponding to single species that temporarily survive out of the nine species in total. In [18] we have demonstrated how an appropriate choice of predation rates can steer the trajectory through the heteroclinic network in a desired way. The choice of predation rates is determined by conditions on the eigenvalues of the Jacobian at the saddles. In particular we observed a separation of time scales between the (fast) oscillations within the SHCs and the (slow) oscillations in the LHC that modulate the fast oscillations such that the scales differ by an order of magnitude. It is this aspect that is of interest in view of applications to neuronal dynamics, where a modulation of fast oscillations via slow ones is commonly observed among the multiple time scales in the brain.
In this paper we derive a Poincaré map that allows to predict how many revolutions the trajectory stays within an SHC before it switches to a heteroclinic connection of the LHC, provided that we start already in the vicinity of a saddle. The subsequent evolution of the trajectory is not covered by this map: after the switch towards a heteroclinic connection of the LHC it may either turn to a new, second SHC, or stay within the LHC. In section 2 we present the model with the choice of predation rates, in section 3 we derive the Poincaré map and present its predictions in section 4 before concluding in section 5.

The model
We consider a hierarchical heteroclinic network as described in [18], which consists of three small heteroclinic cycles (SHCs) that are the saddles of a large heteroclinic cycle (LHC), c.f. fig. 1. The defining equation where s i ∈ R + 0 in the original Lotka-Volterra context denotes the concentration of species i, γ the death rate, while the set of A i, j constitutes the predation matrix A, which determines the topology of the heteroclinic network. For simplicity, we define the index j for given i of A i j directly by the five index functions j = es(i) = ((i + 2) mod 9) + 1 (that returns the index of the expanding small direction es(i)), in analogy el(i) = (i mod 3) + 1 + 3 i−1 3 (expanding large), cs(i) = es 2 (i) (contracting small), cl(i) = el 2 (i) (contracting large) and tv(i) = {1, . . . , 9} \ {i, es(i), el(i), cs(i), cl(i)} for transverse directions at saddle σ i . The resulting topology of the hierarchical heteroclinic network is displayed in fig. 1. We set A i,i = 0, A i,cs(i) = c, A i,cl(i) = d, A i,es(i) = e, A i,el(i) = f , and A i,tv(i) = r. By this degeneracy of parameter values we guarantee Z 3 × Z 3 permutation symmetry (both between species within one SHC and between the SHCs). As we have shown in [18], we can steer the trajectory along selected paths of the heteroclinic network and guarantee its stability via conditions on the eigenvalues of the Jacobian at the corresponding saddles. These conditions imply for the remaining rate parameters: 0 An important property of eq. (1) is that the coordinate planes are invariant sets. This guarantees the existence of heteroclinic orbits. Usually such invariant planes are the result of reflection symmetries in the equations of motion. Here, such a symmetry is not manifest. However, since we restricted s i ≥ 0, we can extend our definition to include negative values by making the coordinate change s i → r 2 i . This yields symmetric dynamics in the other 2 9 − 1 hyperoctants, related by the reflections ρ i : r i → −r i . In the new coordinates (and after rescaling time by a factor of 2 for convenience) the equations (1) read They are now equivariant under the nine reflections ρ i , i.e., ∂ t (ρ k r i ) = ρ k (∂ t r i ). Due to this, the heteroclinic network is structurally stable to perturbations as long as they do not break these symmetries. Another effect of this coordinate change is a simplification of the linearized dynamics at the saddles. The Jacobian of eq. (3), evaluated at a saddle (r i = 0 ∧ r j = 0 ∀ j = i), has only entries on its diagonal, while the Jacobian of eq. (1) has additional entries in one row. In the following section we shall see why this simplification is useful.
While in all of the following we work with eq. (3), our results are still valid for the special case s i ≥ 0 ∀i, corresponding to eq. (1) after reversing the coordinate change.

Poincaré sections and return map
Our goal is to construct a Poincaré map to determine the number of revolutions in an SHC before the trajectory turns to the heteroclinic connection of the LHC. We derive the return map by proceeding in the standard way, e.g. similarly to [19]. First, we determine local maps φ abc that characterize the dynamics in the vicinity of each saddle. Then, we argue about the form of the global maps ψ bc , which transport the trajectory between the saddles' vicinities. Next, we compose local map and global map to find the transition map that describes one third of a revolution. Finally, three applications of the transition map yield the full return map.
For the local map, we linearize the dynamics at the example of saddle 1 (r * . Therefore we use local coordinates, i.e., we introduce u 1 = r 1 − γ −1/2 and consider {u 1 , r 2 . . . r 9 } small, so that terms of order 2 will be neglected. Thus, the linearized equations are At saddle 1, the trajectory takes one of two paths: either to saddle 2 (staying inside the SHC) or onwards to saddle 4, c.f. fig. 2(a). Each possibility is described by its own local map, φ 312 and φ 314 , respectively. The domain and codomains of these maps are the cross sections H in,3 The value of h > 0 thereby is the distance of the Poincaré sections from the saddles. It must be chosen small enough that linearizing the dynamics is justified and the local maps φ i jk approximate the true dynamics sufficiently well. However, if h is chosen too small, the dynamics at the other side of the section is not well approximated by the linearized dynamics of the global map (see below). With respect to the Z 3 × Z 3 symmetry, we set h the same for all cross sections without loss of generality.
As a result, taking only lowest order terms into account we find for some constants a k , b 3 , . . . b 9 . Thus, there is no "mixing" of variables due to the reflection symmetries (apart from the radial coordinate). Apparently Z n 2 symmetries in general support the assumption that the global map acts as the identity [20]. We can thus derive the transition map g i = ψ i,es(i) • φ cs(i),i,es(i) by describing its action on each coordinate, to find the variables r j at time T 2 , i.e., r j,2 ∈ H in,i es(i) .
with f (r i,0 , r es(i),0 ) = a 0 + a 1 (r i,0 − γ −1/2 )( r es(i),0 h ) 2γ/E + O(r 2 j ). The second indices 0 and 2 represent in general the time instants before and after application of the map g i . Upon the evolution in time, the role of the variables r i , i ∈ {1, . . . , 9} changes with respect to their directions (es, el, cs, cl). What was a cs(i)-direction before the map, becomes the direction i at the subsequent saddle. Thus we apply cs(i) to the indices at time T 2 (terms in the middle of eqs. (12) to (14)), using cs(es(i)) = i, cs(cs(i)) = es(i), etc. This way we trace back the origin of the new directions. The arguments of the double primed terms show how the arguments transform under a generic g i . For example, r es(i),0 transforms to b cs(i) r c−γ γ−e es(i),0 h, which is the new variable r cs(i),2 assigned to the contracting small direction at time T 2 .
From eqs. (5) and (8) it seems convenient to rescale the variables r i by h −1 and time by h. We thereby remove the h-dependence of the condition entering eq. (5). A short calculation with eq. (3) shows that this transformation requires to rescale all parameters (γ, c, d, e, f and r) by a factor of h 2 . To simplify the notation we introduce the new parameters Furthermore, it is convenient to rename the variables in a way that they refer at each saddle to the same kind of direction (es, el, cs, cl, etc.). We thus identify x = r es(i),0 , y 0 = r el(i),0 , y n = r es(el(i)),0 , y p = r cs(el(i)),0 .
The reason why we ignore the remaining variables will become clear in eq. (20) where B 0 , B n , B p and A remain as constants. For example, the components resulting from an application of g 1 can be read off from the transition map eqs. (12) to (14) as follows: the transformation of x follows from the left equation in eq. (14), that of y 0 from the right equation in eq. (12), using the index relations for index 0 = el(i) = cs( j) = j − 1, so that j = el(i) + 1 = index n, of y n again from the right equation in eq. (12), now with n = es(el(i)) = cs( j) = j − 1, so that j = el(i) + 2 = index p, and for y p using the right equation from eq. (13) and p = cs(el(i)). In the same way the actions of g 2 and g 3 can be derived.
The condition x F E > y 0 must be fulfilled at every saddle. In total this yields as conditions for the domain of g(x, y 0 , y n , y p ). If any condition of eq. (20) is violated, the trajectory will leave the SHC during this revolution at the respective saddle. At this point it is clear why the restriction to four variables is justified, as no other variables enter the conditions (20). All information about whether (and when) the switch to the next SHC happens is contained in x, y 0 , y n and y p .

Results on the ratio between time scales
In the following, we estimate how long (in terms of revolutions) the trajectory spends near one SHC before switching to the next. Roughly three times this duration defines the time scale of an LHC. This relation defines a kind of ratio between time scales.
The idea behind the subsequent calculation is the following. Let us assume that we start at point 0 of fig. 3(a). By applying the return map g, we find the coordinates of the next point 1. The system is still in the domain of the SHC return map (below the line x F E ), so that a further application of g is in order. After the nth application, the line x F E is crossed. Here, the system switches to the next SHC. We start by recalling the condition for following the heteroclinic connections of the SHC, given by the domain of the local map eq. (5) together with the definitions eq. (15): Note that the relevant variable z is the ratio between orders of magnitude of the two "deviations" x and y 0 (more precisely, the distances in expanding directions from the saddles in the SHC and LHC). Every revolution corresponds to applying the return map once. From eq. (19) we find where α = C 3 E 3 and β = − F E − RC E 2 − RC 2 E 3 , B = B 0 and A as before. Note that the undetermined constants A and B can be traced back to eq. (19) to the parameters b k in eq. (11), which remain undetermined. Commonly, these constants from the global maps are either irrelevant for the result (when only stability or convergence statements should be made, as e.g. in [19,20] (as e.g. in [21]). If we numerically determine A and B, for our choice of parameters they turn out to be of order 1. Therefore neglecting their logarithms in eq. (22), we obtain with α and β as given above. We now repeatedly apply eq. (23) to eq. (21): which generalizes to where we used the geometric series. Assuming equality and solving for n we find a closed expression for the approximate number of revolutions that the system stays near the SHC for given parameters and initial condition z > F E : Note that z ∈ H in,cs(i) i is an initial value on a cross section already near the heteroclinic network. For arbitrary initial conditions, first transient dynamics lead towards the heteroclinic network and then determine the initial values of z.
Up to now we have only taken into account y 0 and neglected the remaining conditions in eq. (20). Thus, any departure from one of the other two saddles will be perceived as the subsequent departure at saddle 1, introducing an error in n of up to 2 3 . To account for the other two conditions we may either repeat the above derivation accordingly and separately, starting from eq. (20); or, as we proceed in the following, transforming the initial condition so that we start g not at σ 1 , but at σ 2 or σ 3 . We thus apply g 1 and g 2 • g 1 to (x, y 0 , y n , y p ) and deduce As a result, for any initial condition (x, y 0 , y n , y p ) ∈ H in,3 1 we find the number of revolutions n during which the system stays near the SHC as In summary, the following approximations entered our calculation: at this distance h is best located in view of applying the local map on one side of the section and the global map on the other side of the section. Otherwise the domains are not optimally located with respect to the approximations to apply. Yet, the error that is introduced by A = B = 1 is small (n pred. ≈ 4.97 vs. n num. ≈ 5 in fig. 3(a)). Since setting A = B = 1 enables the derivation of a closed expression at all, it is certainly useful to estimate the number of revolutions in this way. (It should noticed that this number anyway fluctuates of the order of one as a function of the very initial conditions, chosen in the beginning of the time evolution, which then -after transient dynamics -lead to different initial values for z [18].) So far the description is complete up to the point of leaving the current SHC. For the subsequent time evolution we are interested in the questions of (i) whether the system will enter the next SHC or rather continue along an LHC connection and (ii) how long (in terms of revolutions) the system will stay near the second SHC if it goes there. Again let us start at saddle σ 1 w.l.o.g. Leaving the first SHC corresponds to applying the local map φ 314 and the global map ψ 14 thereafter. However, the form of this escape map G 1 = ψ 14 • φ 314 is not covered by the generic form of g i if it should be used for answering the former questions, though locally the condition at saddle σ 4 is the same as the first one in eq. (20) (x F E > y 0 ): When the trajectory approaches H out,5 4 at distance r 5 = h, r 7 is required to be smaller than h to continue along the LHC. However, if this prediction should be made at earlier times, from coordinates before the trajectory escapes at saddle 1 to saddle 4, the map has to keep track on r 7 ≡ w 0 that enters condition eq. (31) This condition can be derived from the (locally valid) condition (x F E > y 0 ) at saddle 4 by taking the inverse of G 1 = ψ 14 • φ 314 to obtain the coordinates prior to the escape from saddle 1. The condition (31) makes manifest the need for r 7 .
To answer the second question as to the number of revolutions within the second SHC, further coordinates r 8 and r 9 must be pursued from the beginning in order to satisfy analogous conditions to the latter two in eq. (20). Accordingly the escape map G i = G 1 depends on seven coordinates G i (x, y 0 , y n , y p , w 0 , w n , w p ) = Ā y of which we worked out the four locally relevant ones: apart from y 0 , playing now the role of the former x at saddle 1, it depends on three more variables, w 0 , w n and w p corresponding to r 7 , r 8 and r 9 respectively. These are taking over the role of the former y 0 , y n and y p , respectively, once the trajectory arrives at H out,5 4 . The escape from the first SHC, represented by G i , is visualized in fig. 3(b), together with the following revolutions within the second SHC consisting of σ 4 , σ 5 and σ 6 . It should be further noticed that revolutions within the second SHC start with φ 145 , unlike subsequent visits of σ 4 where φ 645 applies. Moreover, from fig. 3(b) it can be seen that the "distances" y i of the trajectories to the LHC increase during revolutions within a single SHC, but decrease from one SHC to the next SHC, while the distances x to the heteroclinics of both SHCs decrease upon further revolutions, both results in accordance with the expectation.

Conclusions and Outlook
In summary, we have demonstrated how to construct a Poincaré map for a hierarchical heteroclinic network. We worked out the map in detail for predicting the number of revolutions within a small heteroclinic cycle, before the trajectory turns to the heteroclinic connection in the large cycle. The map reproduced the sensitivity to the choice of certain rate parameters that we had numerically identified in previous work as the essential ones for tuning the time scale separation between slow and fast oscillations. It also showed explicitly that the variation of the number of revolutions within an SHC depends on the initial conditions just before entering the SHC. We also indicated how to include more coordinates in the Poincaré map if its prediction should comprise further decisions of the future time evolution of the trajectory (e.g., to turn to a second and third SHC).
A detailed understanding of the time scale separation may not only be relevant for applications in neuronal dynamics, but also useful in connection with heteroclinic computing [4]. When heteroclinic computing is realized with winnerless competition, a tuned number of revolutions in a heteroclinic cycle may cause a certain amount of delay. Independently of whether the delay is desired or not, it is calculable via the Poincaré map.
If we furthermore include noise in our model of winnerless competition, it is an open and interesting question what determines the dwell time for a stay near a saddle if the saddle is a heteroclinic cycle itself. Our preliminary simulations indicate that a naive generalization of results of [2] is not applicable. According to the results of [2], the time would be determined by the logarithm of the noise intensity and the eigenvalue of the most unstable direction. Our results obtained via the Poincaré map -so far without noise -may provide a starting point for calculating the dwell time under inclusion of noise.