Predicting the separation of time scales in a heteroclinic network

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.

We consider a hierarchical heteroclinic network as described in [18], which consists of three small hetero-clinic cycles (SHCs) that are the saddles of a large heteroclinic cycle (LHC), c.f. fig. 1. The defining equation is

Fig. 1

where Si∈ℝ0+ ${{S}_{i}}\in \mathbb{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 hetero-clinic network. For simplicity, we define the index j for given i of A_ij 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)=(imod3)+1+3⌊ i−13 ⌋ $el\left( i \right)=\left( i\,\bmod 3 \right)+1+3\left\lfloor \frac{i-1}{3} \right\rfloor$(expanding large), cs(i) = es²(i) (contracting small), cl(i) = el²(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 Ai,i=0,Ai.cs(i)=c,Ai.cl(i)=d, ${{A}_{i,i}}=0,{{A}_{i.cs\left( i \right)}}=c,{{A}_{i.cl\left( i \right)}}=d,$ Ai,es(i)=e,Ai,el(i)=f,and Ai,tv(i)=r. ${{A}_{i,es\left( i \right)}}=e,{{A}_{i,el\left( i \right)}}=f,\,\,\,\text{and }{{A}_{i,tv\left( i \right)}}=r.$By this degeneracy of parameter values we guarantee ℤ3×ℤ3 ${{\mathbb{Z}}_{3}}\times {{\mathbb{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:

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 Si↦ri2. ${{S}_{i}}\mapsto r_{i}^{2}.$This yields symmetric dynamics in the other 29−1 hyperoctants, related by the reflections ρi:ri↦−ri. ${{\rho }_{i}}:{{r}_{i}}\mapsto -{{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(ρ_kr_i) = ρ_k(∂_tr_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 (ri≠0∧rj=0∀j≠i), $\left( {{r}_{i}}\ne 0\wedge {{r}_{j}}=0\forall j\ne i \right),$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.

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 (r1*=γ−1/2,rj*=0∀j≠1). $\left( r_{1}^{*}={{\gamma }^{-{1}/{2}\;}},r_{j}^{*}=0\forall j\ne 1 \right).$Therefore we use local coordinates, i.e., we introduce u₁ = r₁ −γ^−1/2 and consider {u₁, r₂ . . . r₉} 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, ϕ₃₁₂ and ϕ₃₁₄, respectively. The domain

Fig. 2

and codomains of these maps are the cross sections H1in,3={ (u1,r2…r9):r3=h },H1out,2={ (u1,r2…r9):r2= $H_{1}^{\operatorname{in},3}=\left\{ \left( {{u}_{1}},{{r}_{2}}\ldots {{r}_{9}} \right):{{r}_{3}}=h \right\},H_{1}^{\text{out,2}}=\left\{ \left( {{u}_{1}},{{r}_{2}}\ldots {{r}_{9}} \right):{{r}_{2}}= \right.$h} and H1out,4={ (u1,r2…r9):r4=h }. $H_{1}^{\text{out,4}}=\left\{ \left( {{u}_{1}},{{r}_{2}}\ldots {{r}_{9}} \right):{{r}_{4}}=h \right\}.$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 ℤ3×ℤ3 ${{\mathbb{Z}}_{3}}\times {{\mathbb{Z}}_{3}}$symmetry, we set h the same for all cross sections without loss of generality.

We now integrate the linearized equations (4) from time T₀ of crossing H1in,3 $H_{1}^{\operatorname{in},3}$to time T₁ of hitting either H1out,2 $H_{1}^{\text{out,2}}$or H1out,4 $H_{1}^{\text{out,4}}$(whichever is first determines the subsequent path). We use the notation r_i,₀ = r_i(T₀) and r_i,₁ = r_i(T₁) (and equally u_i,₀ = u_i(T₀) and u_i,₁ = u_i(T₁)). First we focus on the local map ϕ₃₁₂ (at σ₁, coming from σ₃ and going to σ₂), i.e., staying in the SHC. Figure 2(b) shows a sketch. The integration yields

The map ϕ₃₁₄ looks similar, but with its fourth coordinate replaced by h (it maps to H1out,4 ) $\left. H_{1}^{\text{out,4}} \right)$and its second coordinate by (r4,0h)−γ−eγ−fr2,0<h. ${{\left( \frac{{{r}_{4,0}}}{h} \right)}^{-\frac{\gamma -e}{\gamma -f}}}{{r}_{2,0}}<h.$Note that the condition r2,0γ−fγ−e>r4,0hγ−fγ−e−1, $r_{2,0}^{\frac{\gamma -f}{\gamma -e}}>{{r}_{4,0}}h\frac{\gamma -f}{\gamma -e}-1,$distinguishing the domains of both maps, reflects the fact that while the second coordinate of ϕ₃₁₂ has already reached a distance h to the saddle, the fourth coordinate of ϕ312,(r2,0h)−γ−fγ−er4,0, ${{\phi }_{312}},\left( \frac{{{r}_{2,0}}}{h} \right)-\frac{\gamma -f}{\gamma -e}{{r}_{4,0}},$lags behind, so it is still smaller than h.

We abstract from this concrete case to the general local mappings ϕ_{cs(i),i,es(i)} utilizing the index functions from above. From eq. (5) we observe that

Next, we turn to the global map ψ1,2:H1out,2→H2in,1. ${{\psi }_{1,2}}:H_{1}^{\text{out,2}}\to H_{2}^{\operatorname{in},1}.$We will assume that H2in,1 $H_{2}^{\operatorname{in},1}$is reached at time T₂, and define r _j(T₂) = r _j,₂. Since the original equations (3) are equivariant under the reflections ρ_i, the global maps must as well fulfill this property, i.e., ψ₁,2(ρ_ir) = ρ_iψ1,2(r). More specifically: if

then by equivariance we also have

As a result, taking only lowest order terms into account we find

for some constants a_k,b₃, . . .b₉. Thus, there is no “mixing” of variables due to the reflection symmetries (apart from the radial coordinate). Apparently ℤ2n $\mathbb{Z}_{2}^{n}$symmetries in general support the assumption that the global map acts as the identity [20]. We can thus derive the transition map gi=ψi,es(i)∘ϕcs(i),i,es(i) ${{g}_{i}}={{\psi }_{i,es\left( i \right)}}\circ {{\phi }_{cs\left( i \right),i,es\left( i \right)}}$by describing its action on each coordinate, to find the variables r _j at time T₂, i.e., r _j,₂ ∈Hes(i)in,i. $\in H_{es\left( i \right)}^{\operatorname{in},i}.$

with f(ri,0,res(i),0)=a0+a1(ri,0−γ−1/2)(res(i),0h)2γ/E+O(rj2). $f\left( {{r}_{i,0}},{{r}_{es\left( i \right),0}} \right)={{a}_{0}}+{{a}_{1}}\left( {{r}_{i,0}}-{{\gamma }^{-{1}/{2}\;}} \right){{\left( \frac{{{r}_{es\left( i \right),0}}}{h} \right)}^{{2\gamma }/{E}\;}}+O\left( r_{j}^{2} \right).$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₂ (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 bcs(i)res(i),0c−γγ−eh, ${{b}_{cs\left( i \right)}}r_{es\left( i \right),0}^{\frac{c-\gamma }{\gamma -e}}h,$which is the new variable r_cs(i),2 assigned to the contracting small direction at time T₂.

From eqs. (5) and (8) it seems convenient to rescale the variables r_i by h⁻¹ 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². 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

The reason why we ignore the remaining variables will become clear in eq. (20) below. The transition map, restricted to the four variables, gi(x,y0,yn,yp):Hiin,cs(i)→Hes(i)in,i, ${{g}_{i}}\left( x,{{y}_{0}},{{y}_{n}},{{y}_{p}} \right):H_{i}^{\operatorname{in},cs\left( i \right)}\to H_{es\left( i \right)}^{\operatorname{in},i},$can now be applied repeatedly. In this way we construct the return map g(x,y0,yn,yp):Hiin,cs(i)→Hiin,cs(i) $g\left( x,{{y}_{0}},{{y}_{n}},{{y}_{p}} \right):H_{i}^{\operatorname{in},cs\left( i \right)}\to H_{i}^{\operatorname{in},cs\left( i \right)}$by stepwise composition:

where B₀,B_n,B_p and A remain as constants. For example, the components resulting from an application of g₁ 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₀ 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₂ and g₃ can be derived.

The condition xFE>y0 ${{x}^{\frac{F}{E}}}>{{y}_{0}}$must be fulfilled at every saddle. In total this yields

as conditions for the domain of g(x, y₀, 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₀, y_n and y_p.

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 xFE ${{x}^{\frac{F}{E}}}$), so that a further application of g is in order. After the nth application, the line xFE ${{x}^{\frac{F}{E}}}$is crossed. Here, the system switches to the next SHC.

(b)

Fig. 3

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₀ (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 α=C3E3 $\alpha =\frac{{{C}^{3}}}{{{E}^{3}}}$and β=−FE−RCE2−RC2E3, $\beta =-\frac{F}{E}-\frac{RC}{{{E}^{2}}}-\frac{R{{C}^{2}}}{{{E}^{3}}},$B = B₀ 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]), or they are determined numerically if their exact values matter for the result (as e.g. in [21]). If we numerically determine A and B, for our choice of parameters they turn out to be of order

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>FE: $z>\frac{F}{E}:$

Note that z∈Hiin,cs(i) $z\in H_{i}^{\operatorname{in},cs\left( i \right)}$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₀ 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 23. $\frac{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 σ₁, but at σ₂ or σ₃. We thus apply g₁ and g₂ ◦ g₁ to (x, y₀, y_n, y_p) and deduce

As a result, for any initial condition (x,y0,yn,yp)∈H1in,3 $\left( x,{{y}_{0}},{{y}_{n}},{{y}_{p}} \right)\in H_{1}^{\operatorname{in},3}$we find the number of revolutions n during which the system stays near the SHC as

In summary, the following approximations entered our calculation: Linearization of the dynamics in the vicinity of the saddles (resulting in local maps), linearization of the global maps, and fixing the remaining constants A and B to the value 1. The justification of the linearizations as well as setting A = B = 1 depend on the choice of h. Note that while the distance h of Poincaré sections from the saddles cancels out of the final results (eqs. (25) and (30)), the constants A and B do depend on h. How h should be chosen so that A = B = 1 is a good approximation is subtle. As our tests have shown, if h is chosen too small, the approximation works less well. There is an intermediate value of h (such as h = 0.1) so that the Poincaré section Hiin,cs(i) $H_{i}^{\operatorname{in},cs\left( i \right)}$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 σ₁ w.l.o.g. Leaving the first SHC corresponds to applying the local map ϕ₃₁₄ and the global map ψ₁₄ thereafter. However, the form of this escape map G₁ = ψ₁₄ ◦ ϕ₃₁₄ 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 σ₄ is the same as the first one in eq. (20) (xFE>y0): $\left( {{x}^{\frac{F}{E}}}>{{y}_{0}} \right):$When the trajectory approaches H4out,5 $H_{4}^{\text{out,5}}$at distance r₅ = h, r₇ 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₇ ≡ w₀ that enters condition eq. (31)

This condition can be derived from the (locally valid) condition (xFE>y0) $\left( {{x}^{\frac{F}{E}}}>{{y}_{0}} \right)$at saddle 4 by taking the inverse of G₁ = ψ₁₄ ◦ ϕ₃₁₄ to obtain the coordinates prior to the escape from saddle 1. The condition (31) makes manifest the need for r₇.

To answer the second question as to the number of revolutions within the second SHC, further coordinates r₈ and r₉ 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₁ depends on seven coordinates

of which we worked out the four locally relevant ones: apart from y₀, playing now the role of the former x at saddle 1, it depends on three more variables, w₀,w_n and w_p corresponding to r₇, r₈ and r₉ respectively. These are taking over the role of the former y₀, y_n and y_p, respectively, once the trajectory arrives at H4out,5. $H_{4}^{\text{out,5}}.$

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 σ₄,σ₅ and σ₆. It should be further noticed that revolutions within the second SHC start with ϕ₁₄₅, unlike subsequent visits of σ₄ where ϕ₆₄₅ 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.

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.