Heteroclinic units acting as pacemakers: entrained dynamics for cognitive processes

When the unstable manifold of a saddle equilibrium intersects the stable manifold of another saddle, the intersection is called a heteroclinic orbit. A heteroclinic network is a set of vertices, corresponding to saddles, connected by edges, which are heteroclinic orbits. Thus heteroclinic networks are networks in phase space. Saddles in these networks are not restricted to saddle equilibria, but refer to any invariant set with non-trivial stable and unstable manifolds. Early seminal papers on heteroclinic networks are [23, 24], see, for example, also [25]. A simple special case is a heteroclinic cycle between saddle equilibria, in which the unstable direction of one saddle becomes the stable direction of the subsequent saddle. The trajectory associated with the heteroclinic cycle approaches the vicinity of the saddles in a cyclic way. Upon its revolutions, it comes closer and closer to the saddles, but slows down and reaches the heteroclinic cycle in the infinite-time limit. The slowing down may be weakened or entirely prevented by exposing the system to a small amount of noise.

In view of our applications we couple heteroclinic networks on a spatial grid, as considered in [12]. To disentangle heteroclinic networks (which are networks in phase space) from network aspects of the spatial grid, we term heteroclinic networks assigned to a single site of the grid a heteroclinic unit. Thus the coupling of heteroclinic units refers to their coupling in coordinate space. This construction allows to share typical features of brain dynamics: according to [26], the overall organization of the brain is in terms of parallel processing streams with complementary properties, hierarchical interactions within each stream, and parallel interactions between the streams. Our assignment of heteroclinic units to different sites of the grid realizes such a performance that is ongoing in parallel at different locations, while the imposed hierarchy of individual heteroclinic networks in phase space implements processes with hierarchies in time scales such as the modulation of fast oscillations via slow oscillations.

The system of heteroclinic units is defined on regular grid topologies with nodes labelled by k. To each node, we assign a heteroclinic unit that obeys GLV-equations according to

$$\begin{equation}{\partial }_{t}{s}_{k,i}=\rho {s}_{k,i}-{\gamma }_{k}{s}_{k,i}^{2}-\sum\limits _{j\ne i}\,{A}_{i,j}{s}_{k,i}{s}_{k,j}+\sum\limits _{l}\,{K}_{k,l}({s}_{l,i}-{s}_{k,i})+{\sigma }_{k}\vert {\xi }_{i}(t)\vert .\end{equation} \tag{ 1 }$$

The meaning of the variables s_k,i depends on the intended application. In an ecological context, s_k,i represent the concentration of the ith species of the kth GLV unit, ρ is the net reproduction rate, chosen as 1 to set the time scale, γ_k alters the net reproduction rate of the kth unit in a density-dependent way. Throughout the paper we focus on possible applications to brain dynamics, in particular transient cognitive processes. Here we interpret s_k,i as densities of neural sub-populations. During the winnerless competition, a sequence of saddles is approached. In the vicinity of individual saddles various specific sub-populations get temporarily dominantly excited. Information (for example from sensory input) is believed to be encoded in the specific spatiotemporal excitation pattern. In the vicinity of a saddle, the variable s_k,i (t) tells us which subpopulation (i) gets dominantly excited at what location (k) and at what time. Thus, information is contained in the very selection of the specific saddle and its time coordinate within the ordered temporal sequence of approached saddles as well as in the sequence itself. In short, we term the variables s_k,i the information items or just the 'items' in the following, without specifying the neural populations any further.

The parameter ρ is again the net reproduction rate and γ_k the item-density dependent decay rate, σ_k the strength of the additive noise to the kth unit, and ξ is Gaussian white noise with zero mean. A is called the rate matrix with A_i,j being the competition rate with which the item j of a unit acts on item i of the same unit. The interaction between the units is determined by the coupling matrix K: K_k,l is finite if unit l influences unit k, otherwise it is zero.

2.2.1. GLV-equations for three items

When each unit comprises of three items and the items play a rock-paper-scissors game, the rates characterize intransitive (non-hierarchical) competition. For the brain dynamics they characterize the 'competition' strength in a winnerless game of items, the winnerless competition among its items then results in a single hierarchy level (one time scale), and A is a 3 × 3 matrix given by

$$\begin{equation}A=\left(\begin{matrix}\hfill 0\hfill & \hfill c\hfill & \hfill e\hfill \\ \hfill e\hfill & \hfill 0\hfill & \hfill c\hfill \\ \hfill c\hfill & \hfill e\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ 2 }$$

As studied previously [10, 11], there is a critical value of the decay rate γ : γ_crit = (c + e)/2, at which a single unit undergoes a Hopf bifurcation. For γ_k < γ_crit, the dynamics of the kth unit is characterized by a stable heteroclinic cycle connecting the saddle equilibria ${\xi }_{i}=\left\{\vec{{s}_{k}}\in {\mathbb{R}}^{3}:{s}_{k,i}=\rho /{\gamma }_{k},{s}_{k,j}=0\ \forall \,j\ne i\right\}$. For γ_k > γ_crit, the coexistence equilibrium ${\xi }^{\ast }=\left\{\vec{{s}_{k}}\in {\mathbb{R}}^{3}:{s}_{k,i}=\frac{\rho }{c+e+{\gamma }_{k}}\ \forall \,i\right\}$ is stable. Throughout the manuscript, we fix c = 2.0 and e = 0.2, and therefore, γ_crit = 1.1. At γ = γ_crit, due to an emerging conservation of the sum and product of three items for γ_crit, the system has an infinity of periodic orbits, depending on the choice of initial conditions, with typical signatures of critical behavior [27].

2.2.2. GLV-equations for nine items

When each unit contains nine items, the rate matrix A is chosen in such a way that the winnerless competition among its items has an attractor that is a (large) heteroclinic cycle (LHC) between three saddles that are themselves (small) heteroclinic cycles (SHCs) between three saddle equilibria. As we have shown in [10, 11], this structure of the attractor induces a hierarchy in time scales of slow oscillations (due to LHCs) modulating fast oscillations (due to SHCs). The long (short) time scale amounts to one revolution in the long (short) heteroclinic cycle, respectively, which is well defined as long as the slowing down is suppressed by the application of a small amount of noise, otherwise it is itself a function of time. Such a modulation of oscillations is a typical feature that is also experimentally observed in brain dynamics (see, for example, [28, 29]). The rate matrix A for the units with two hierarchy levels is then given as the following block matrix

$$\begin{equation}A=\left(\begin{matrix}\hfill {m}_{0}\hfill & \hfill {m}_{d}\hfill & \hfill {m}_{f}\hfill \\ \hfill {m}_{f}\hfill & \hfill {m}_{0}\hfill & \hfill {m}_{d}\hfill \\ \hfill {m}_{d}\hfill & \hfill {m}_{f}\hfill & \hfill {m}_{0}\hfill \end{matrix}\right),\quad \text{where}\;{m}_{0}=\left(\begin{matrix}\hfill 0\hfill & \hfill c\hfill & \hfill e\hfill \\ \hfill e\hfill & \hfill 0\hfill & \hfill c\hfill \\ \hfill c\hfill & \hfill e\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ 3 }$$

The matrix m_d is also 3 × 3 and has rates d on the diagonal elements and r on the remaining elements. Similarly, the matrix m_f has a rate f on the diagonal with all off-diagonal elements being r. Again, on tuning the decay rate γ, the system undergoes a sequence of Hopf bifurcations, whose order and values depend on the choice of the other parameters [10, 11]. Two types of Hopf bifurcations are relevant to us here. For low values of γ, the system has an intricate heteroclinic dynamics with two hierarchy levels. When γ is increased towards a bifurcation point γ_c, the system undergoes a Hopf bifurcation, at which the highest hierarchy level is gone. What remains is an LHC between three three-items saddle equilibria as a remnant from the SHCs. Further increasing γ, the system undergoes a second Hopf bifurcation which drives the unit into stable global coexistence of nine items with the same concentrations. In extension of such a single hierarchical heteroclinic network, the dynamics on a spatial network (a two-dimensional square lattice) of identical such heteroclinic units was analyzed in [12], where the units were coupled via diffusion.

The very construction of the adjacency matrix A in phase space may appear rather peculiar. It serves to enforce hierarchical heteroclinic motion just by the choice of rates. The rates are chosen out of certain intervals, so that the eigenvalues of the corresponding Jacobians at the saddle equilibria lead to a preference of the trajectory for one (desired) direction out of two unstable directions when escaping from the saddles. This construction was based on results of [30]. A hierarchy in time scales can be implemented differently, see e.g. [5]. In the brain, external input may be responsible for selecting certain paths in the (counterparts of) heteroclinic networks. We do not consider external input here. Our choice is just for convenience to generate an intricate motion that should be entrained, independently of how it was generated.

In chemical oscillatory systems, the role of pacemakers seems to be played by defects and impurities [31]. Accordingly, in the simplest case we treat the heteroclinic unit which should play the role of a pacemaker on the same footing as all other units (the driven units), but assign a 'defect' in the sense that we distinguish its individual parameters from those of all other units: for given parameters c and e we choose γ in the range of the highest (here at most second) hierarchy level of motion, so that this unit performs either a heteroclinic cycle of heteroclinic cycles (nine items) or a heteroclinic cycle between saddle equilibria (three items), while the other units are chosen in the range of γ for which the dynamics converges to a stable equilibrium point with coexisting items. In the following we call the corresponding parameter regimes the heteroclinic regime (HC) or the coexistence regime (CE). In addition, the pacemaker may be distinguished by its location and its coupling to the driven units. Along a chain, for local coupling, the pacemaker will be distinguished by its position at one of the edges of the chain and directed local coupling along all units with a smaller back-coupling than forward coupling. When we close the chain to a ring, we assign a smaller back-coupling to the pacemaker at least along the last link that leads to the closure. In [19] we discuss an implementation also with distance-dependent couplings of driven units to the pacemaker https://stacks.iop.org/JPCOMPLEX/3/035003/mmedia. The chains (and rings in [19]) are studied for varying size, in general with stronger forward than backward coupling. On the two-dimensional grid, the coupling is bidirectional (diffusive), in order to compare with earlier results of [12], defects are either randomly assigned or confined to a region outside a small disc of 'non-defects' at the center of the grid. The different options are displayed in figure 1.

Zoom In Zoom Out Reset image size

The nomenclature which we use in this paper is summarized in table 1. Some remarks to the table are in order. We call oscillations to be large (LHCs) or small (SHCs) heteroclinic cycles, even if the characteristic slowing down of the HCs is suppressed by a small amount of noise and it would be more precise to call the motion a quasi-limit cycle in the vicinity of the heteroclinic orbit. This is to distinguish such oscillations from limit cycles remote from the heteroclinic orbit. Such cycles are also encountered in this work, they are called small limit cycles (SLCs). The attributes small or large do not reflect the distance in phase space, but are merely due to the visualization and correspond to hierarchy levels with slow (large cycle) or fast (small cycle) time scales.

For systems with noise and on a two-dimensional grid, we use XMDS2 software [32] for the numerical integration; everywhere else Matlab is used. For analytical calculations, we make use of the Mathematica software.

We start with a system of two 1L-units and three items, with one unit acting as a pacemaker to the other (driven) unit. The pacemaker is at location k = 1 with decay rate γ_P, the driven unit at k = 2 with decay rate γ_D. The governing equation (equation (1)) takes the following form:

$$\begin{align}\hfill {\partial }_{t}{s}_{1,i}& =\rho {s}_{1,i}-{\gamma }_{\mathrm{P}}{s}_{1,i}^{2}-{s}_{1,i}(c{s}_{1,i+1}+e{s}_{1,i+2})+{\delta }_{b}({s}_{2,i}-{s}_{1,i}),\hfill \\ \hfill {\partial }_{t}{s}_{2,i}& =\rho {s}_{2,i}-{\gamma }_{\mathrm{D}}{s}_{2,i}^{2}-{s}_{2,i}(c{s}_{2,i+1}+e{s}_{2,i+2})+\delta ({s}_{1,i}-{s}_{2,i}),\hfill \end{align} \tag{ 4 }$$

where we denote the interaction strength from the pacemaker to the driven unit, that is K_2,1, by δ and the back-coupling from the driven unit to the pacemaker (K_1,2) by δ_b.

We first consider δ_b = 0. The γ_P of the pacemaker is chosen such that it is characterized by a stable heteroclinic cycle that connects three saddles of the pacemaker (ρ/γ_P, 0, 0), (0, ρ/γ_P, 0), (0, 0, ρ/γ_P) via heteroclinic orbits. Due to the characteristic slowing-down of the heteroclinic cycle near the saddles, the trajectory spends most of the time in the vicinity of one of the saddles. During this dwell time, only one of its items is dominating (this is winnerless competition with one temporary winner), while the concentration of the other two items approaches zero. Thus the driven unit feels a constant forcing to one of its items during that time.

This leads us to an effective modeling of the driven unit exposed to constant forcing that we present in detail in appendix A. The result is that depending on the forcing, parameterized by the forward coupling δ, the driven unit approaches a one-item equilibrium (1S), a two-items equilibrium (2S), or a three-items equilibrium (3S), where the numbers (1, 2, 3) refer to the number of coexisting dominating items. In the original model, as the pacemaker approaches one of its saddles, the driven unit approaches one of these equilibria. If the equilibrium approached is a stable equilibrium of the effective model, it stays there until the pacemaker moves towards another of its saddles, thus changing forcing. Otherwise, if it approaches an unstable equilibrium, the driven unit moves towards the stable equilibrium of the effective model while the pacemaker continues to reside in the vicinity of its saddle.

We next consider larger systems of coupled units in a chain configuration as shown in figure 1(a). If we take one 1L-unit as the pacemaker and all other N − 1 units as CE-units, under what conditions will the pacemaker be able to entrain the other units to heteroclinic motion? We consider first unidirectional couplings, bidirectional ones with a strong forward and a weak back-coupling are discussed in [19]. The unidirectional nearest-neighbour couplings along the chain are chosen according to the coupling matrix K, given as

$$\begin{equation*}{K}_{k,l}=\begin{cases}\delta ,\quad \hfill & \text{if}\;k\in \left\{2,\dots ,N\right\}\quad \text{and}\quad l=k-1\hfill \\ 0,\quad \hfill & \text{otherwise}\hfill \end{cases}.\end{equation*}$$

We place a 1L-unit as the pacemaker at the edge of the chain at k = 1 with γ₁ = γ_P < 1.1. For all other CE-units we choose (γ_k = γ_D > 1.1 ∀ k ∈ {2, ..., N}). Note that in this configuration, a unit at location k − 1 acts effectively as a pacemaker to the kth unit which itself drives the unit at k + 1. Therefore, each unit directly or indirectly feels the effect of all the units to its left but is not affected by the dynamics of the units to its right side in the schematic figure 1(a).

In appendix B.1 we zoom into the proliferation of equilibria in such a simple chain as it can be pursued to some extent what causes the proliferation. The dynamics of the nth unit along the chain depends on the dynamics of all the n − 1 units to its left. Even though there are 3n(n + 1)/2 possible equilibria in the phase space of the nth unit that it can visit in principle as shown in table 2, it effectively gets to visit only a few of them. Which path in the branching tree of possible visits it actually visits depends on a number of factors such as

• (a)  
  Which equilibria its (predecessor) pacemaker visits and for how long;
• (b)  
  Whether the approached equilibrium is a stable or unstable equilibrium of the corresponding system with constant forcing. If it gets kicked to an unstable equilibrium, it can move to a more stable one if its experienced forcing does not change during that time interval.
• (c)  
  How closely it approaches an unstable equilibrium. The closer it comes, the longer it stays in the vicinity of the saddle equilibrium.
• (d)  
  The dwell time of the pacemaker at k = 1 in the vicinity of its equilibria.

All these factors in turn depend on the bifurcation parameter of the pacemaker γ_P, the bifurcation parameter of the driven units γ_D, and the coupling along the chain δ. This explains why a precise prediction of the path that is finally taken is practically not feasible for these chains.

Moreover, in appendix B.2 we present the impact of δ, γ_P and γ_D on the entrainment length in some detail. Here is the summary: the stronger the forcing δ, the more units can be entrained, as expected. The more remote a given driven unit is from the original pacemaker at the end of the chain, the more the trajectory wiggles on its approach of the heteroclinic cycle due to the proliferated number of saddles, felt by the trajectory. Furthermore, the entrainment length increases with γ_P, which seems counterintuitive at a first view, and also decreases with an increase in γ_D as expected. For further details we refer to appendix B.

The effect of noise on heteroclinic dynamics can be subtle [30, 33–35]. However, for a certain range of noise strengths (for our other parameter choice this range is between 10⁻¹²–10⁻⁸), its main effect on the heteroclinic dynamics considered here is to prevent the slowing down of heteroclinic dynamics in the neighborhood of the saddles ⁴ . Without noise, we have analyzed the proliferation of possible equilibria in appendix B.1 and table 2. If we add noise to the dynamics of the pacemaker, the slowing down is prevented, thus the forcing that is experienced by the driven units is continually changing. These results are shown for two units in figures 2(a)–(f) for σ = 10⁻¹². For low values of the coupling strength, the driven unit continually switches between its three three-items focus nodes (figure 2(a)), but the stronger the coupling (that is, the forcing), the faster the entrainment to the pacemaker's trajectory and the larger the amplitude, without any further visits to the two- or three-items equilibria after a transient. Finally in (d) and (e), after a very short transient, the driven unit only switches between its three 1S with the same frequency as the pacemaker, being entrained to the heteroclinic dynamics of the pacemaker.

Zoom In Zoom Out Reset image size

The zoom in figure 2(f) clearly shows how an increasing noise strength avoids to reach the vicinity of the saddle equilibrium. While the blue curve (σ = 0) seems to spiral into the saddle before it escapes to the next one, the red curve (σ = 10⁻¹⁰) makes half a revolution of the saddle before it turns off, while the green one (σ = 10⁻⁸) surrounds the saddle at the largest distance, resolving the vicinity of the saddle the least.

A similar effect is seen for larger systems with N > 2 units. In figure 2(g), we show the results for a chain of unidirectionally coupled 1L-units in the presence of noise, applied to the pacemaker. Without noise, the pacemaker would be able to entrain around 100 units before its dynamics slows down, with noise it entrains all 256 units (figure 2(g)). Similar is the result for a small amount of back-coupling along the chain [19], no de-entrainment of the units after a peak of entrainment is reached.

Thus, in summary, both stronger coupling and some amount of noise facilitate the entrainment to the pacemaker dynamics. Also a small back-coupling to the pacemaker can have a similar effect as noise in supporting collective heteroclinic oscillations.

In view of brain dynamics, our result illustrates how a constructive role of noise comes about. The large number of bifurcation points in the phase space of coupled heteroclinic networks may naively suggest that much fine-tuning of the parameters is needed to achieve a state of reproducible partial synchronization. When noise impedes a detailed exploration of the phase space in all its possible bifurcation points, the driven units easily follow a pacemaker, and as 'noise' is ubiquitous, this is a convenient way of facilitating synchronization by exploiting the noise.

In our modeling we do not take into account the thermodynamic cost for achieving precision of the entrainment. In [22] the relation between a stochastic Turing pattern in a Brusselator model and the thermodynamic cost have been addressed. According to the results of [22], the precision of the pattern is maximized for an intermediate thermodynamic cost that is paid for suppressing fluctuations; this means that larger fluctuations seem to have a positive effect on the precision of the pattern. A possible explanation for this observation may be the constructive role of noise on synchronization that we see here for heteroclinic dynamics, while the underlying mechanism is not restricted to heteroclinic dynamics: a certain substructure of the attractor space is not accessible in the presence of noise so that the trajectories of different units in phase space become more similar and synchronize more easily. They cannot explore the close vicinities of the multitude of saddles.

N = 2 units: a pacemaker and a driven unit. For the pacemaker we choose the 2L-unit with γ_P < γ_c, for the driven unit a CE-unit with γ_D > γ_g (γ_g the bifurcation point of the second Hopf bifurcation), unidirectional coupling and sufficiently strong forcing δ, so that the pacemaker (figure 3(a)) can entrain the driven unit towards slow oscillations of fast oscillations. An example is shown in figure 3(b) where the driven unit entrains to the dynamics of the pacemaker till around ten thousand time units. The three items of the first (second, third) SHC of both the pacemaker and the driven unit are shown in shades of red (green, blue), respectively. Beyond time = 10 000, the trajectory of the pacemaker spends so much time near the first item saddle of the second SHC that the trajectory of the driven unit has enough time to go through a sequence of equilibria upon its pursuit to reach the stable equilibrium. This stable equilibrium would be approached in the equivalent system of a unit driven by a constant force. The dynamics is the same as we discussed for the systems with a single hierarchy level. However, if the pacemaker applies a stronger forcing to the driven unit, its trajectory is pushed towards the corresponding saddle and hence spends a long time together with the pacemaker as shown in figure 3(c).

Zoom In Zoom Out Reset image size

As an alternative to applying noise to prevent the slowing down, a small amount of back coupling (δ_b) to the pacemaker has the same effect, as we saw already for the systems with a single hierarchy [19]. Both units then also entrain to heteroclinic dynamics with two levels of hierarchy. Also, the number of revolutions within an SHC increases with an increase in the back coupling. We show an example in figure 3(d) for the driven unit. The pacemaker has similar dynamics, with its amplitude being around 1.

If we increase δ_b beyond a certain value, the first level of hierarchy collapses, both units entrain to an HC between their corresponding three-items saddle equilibria as shown in figure 3(e). If the driven unit is deeper in the CE-region, the second hierarchy level also breaks down as shown in figure 3(f).

N > 2 units without back-coupling. Also a larger system of units with N > 2 can be entrained towards slow oscillations of fast oscillations. We show the results for a unidirectionally coupled chain of 2L-units of size N = 128 in figures 3(g)–(i). The pacemaker is again placed at the edge position k = 1. The units closer to the pacemaker get entrained first, for example, the 30th at around 1000 time steps (figure 3(h)), the 100th unit at around 2000 time steps, again there is a pronounced transient time for the entrainment process to propagate along the chain as revolutions along the cycles of successive units take time. The different shades in panel (g) beyond the three main colors indicate the presence of two hierarchy levels.

So far our main interest was in the option of whether a single unit with oscillatory motion in the vicinity of heteroclinic orbits can entrain other units in a resting state, therefore acting as a pacemaker. The directed coupling, acting more strongly in forward than back-direction and the position of the pacemaker at the edge of a chain were obviously important when a single unit should entrain a whole set of driven units. (For bidirectional (diffusive) coupling, as we consider in section 5, either more pacemakers are needed, or again an asymmetry in the allocation within a whole disc.) Thus it is natural to ask what happens in the inverse case, a single CE-unit at the edge of a chain, directed coupling as before, but 1L or 2L-units as the driven units along the chain. Is a single CE-unit able to stop the heteroclinic oscillations of the entire set of driven units? The answer is positive. Thus it is the combination of the location with the directed coupling and a certain parameter set that enables a single heteroclinic unit to control a larger set of other such units and entrain this set to its own dynamics.

As known from previous results [10, 11] and the last section, the predation rates of a 2L-unit can lead to a trajectory in the attractor space that is schematically indicated in figure 4(a). Corresponding to (a), a possible sequence of saddles, whose neighborhood the trajectory may visit subsequently, would be (1, 2, 3) → (6, 4, 5) → (8, 9, 7) with cyclic repetition; an alternative sequence according to figure 4(b) would be (1, 4, 7) → (8, 2, 5) → (6, 9, 3), cyclically repeated ⁵ .

Zoom In Zoom Out Reset image size

N = 2 units. We consider again two units, the pacemaker with γ_P = 1.05, the driven unit with γ_D = 1.47, δ = 0.1 and δ_b = 0.001.

For most initial conditions, the trajectories of both units follow an LHC of SLCs along the first path as shown in figures 4(c) and (d). ⁶ However, for some initial conditions, the trajectories follow an LHC of SLCs along the second path as shown in figures 4(e) and (f). The LHC of SLCs along the second path is seen for a range of back-couplings which also depends on the choice of the initial conditions. For example, for the initial conditions chosen in (e) and (f), the LHC of SLCs along the second path is seen in the range of δ_b ∈ [0.000 05, 0.023 79]. Right outside this interval, the trajectories instead follow the first path.

N > 2 units. The choice of the second path in the hierarchical heteroclinic network is also seen for larger systems with N > 2 units. Results are obtained for a system of unidirectionally coupled units along a ring of size N = 16 with δ = 1.5 for all links apart from the last one which closes the chain to a ring. Since this coupling amounts to a back-coupling to the pacemaker, it is chosen less strong as δ_P = 0.01. The time evolutions of the trajectories for the pacemaker and the driven units look similar to those in figures 4(c)–(f) for two sets of initial conditions, one choosing the first path and the other choosing the second path. The amplitudes of entrained units are smaller than those of the pacemaker, both when we have to deal with a single driven unit (in case of N = 2) and also as a function of the distance from the pacemaker (visible for larger N, not displayed). The attractor corresponding to the second path seems to have a smaller basin of attraction. We ran the system of rings of heteroclinic units for 100 initial conditions, of which three led to the path according to (b). The corresponding spatiotemporal plots are shown in figures 4(g) and (h), respectively, where we plot the dominant item at each site as a function of time.

In view of possibly shared features with brain dynamics, figures 4(g) and (h) display synchronized units along the entire chain for both paths in the attractor space, and it is the temporal order of switches between subsequently visited saddles that differs between both paths. This means that the corresponding sequence of dominant information items (coded in respective colors) is different, as visible in the spatiotemporal representation. Since the very temporal order of collective patterns is said to encode information [3, 36], our results indicate possibly many more options for encoding information by choosing different paths in a given heteroclinic network. In our implementation, different paths are selected as a result of the rates, the (back)-coupling and the initial conditions (in case of multistability). In the brain, the counterpart of different paths may be selected by external input.

In summary, depending on the selected path in the phase space of hierarchical heteroclinic motion, the encoded information string varies. When it is transferred from the pacemaker along the chain, as a result of entrainment different hierarchical processes are synchronized and going on in parallel at different sites of the grid. This looks like a mathematical realization of the typical features of brain dynamics as conjectured in [26].

As we know from a single heteroclinic unit with two hierarchy levels [10], noise can act similarly to the bifurcation parameter γ. It acts as a control parameter to drive the dynamics from two levels of hierarchy to one level with an effective limit cycle between three three-items saddle equilibria, and further to global coexistence of all items. Similarly here, for example for a chain of N = 64 units, with the pacemaker located at the edge position with γ_P = 1.05 (from the 2L-regime) and all other units with γ_D = 1.47 (from the CE-regime), we pursue the entrainment as a function of an increasing noise strength. Full entrainment to two levels of oscillatory motion is possible for low noise of the order of σ = 10⁻¹²; noise then prevents the slowing down of the pacemaker and supports entrainment to two levels for driven units in the neighborhood of the pacemaker. However, at a more remote distance to the pacemaker, entrainment is achieved only to one hierarchy level, as the influence of the pacemaker dampens with the distance. Stronger noise reduces the levels of the pacemaker from two to one and along with that, the levels of the driven units to a single one at all sites. A further increase of noise then drives the pacemaker to the resting state together with the driven units, whose parameters favor this state anyway. The corresponding figures and more details are presented in [19].

For the first case we assign to each site a value from a uniform distribution in the interval [0, 1] resulting in defects at sites with values less than or equal to p_d . The second case will be discussed subsequently. The decay rate of all possible pacemaker units is taken as γ_x = 0.8, for the defective units we choose a uniform distribution from the interval [1.5, 1.6] of the CE-regime. For different values of the fraction of defects, we observe different patterns on the grid as shown in figure 5. The subfigures (a)–(c) in the upper panel display the color of the dominant item at each location of the grid with species 1–3, 4–6, and 7–9 in different shades of red, green, and blue, respectively.

Zoom In Zoom Out Reset image size

For p_d = 0.1, most of the units are in the 2L-region. The pattern formation in figure 5(a) amounts to nested spirals on the scale of the LHCs, visible in spiral arms of different color, while different shades of a certain color represent the chasing inside an SHC. This means, at each location, two hierarchy levels are present. We show an example in figure 5(d) with the time series of all items of the unit at location x = {32, 15}. The decay rate of this unit is γ_{32,15} = 1.531, so it is a 'defect', but gets entrained to join the collective hierarchical motion. We focus on such a particular unit as its decay rate does not change when we change the fraction p_d in figures 5(b) and (c), because we have used the same random number seeds for the uniform distributions in each case.

Until around p_d = 0.4, the dynamics on the grid is able to maintain two levels of hierarchy, visible in the corresponding time series of winnerless competition between nine items in figure 5(d), while for a larger fraction of defects the lower level of hierarchy is eliminated. We observe LHC spirals in figure 5(b), that is, the spirals of clusters. The LHC connects the three three-items saddle equilibria. As an example, figure 5(e) gives the corresponding time series of all items of the unit at location x = {32, 15}.

For an even larger number of defects, the second hierarchy level also collapses and the global coexistence state stabilizes. In figures 5(c) and (f), we present the results for p_d = 0.95. Though the items of the unit shown in figure 5(f) are moving towards global coexistence, there is still some switching with a rather small amplitude between the three SHCs as shown in the inset.

Zoom In Zoom Out Reset image size

If the 2L-units as pacemakers are confined to a circular region of radius R on the grid, we observe target waves. The results are displayed in figure 6. The size of the grid is again 128 × 128, 2L-units have γ_x = 0.8 inside a circular region of radius R = 8 at the center. All units outside this region are 'defective' (in a resting state) with γ_x chosen randomly from a uniform distribution in the interval [1.5, 1.6]. This distribution of the decay rates and therefore of the defects is color coded in figure 6(a). It leads to the target wave pattern of figure 6(b). The wave fronts on the coarse scale of LHCs propagate inwards, whereas the fronts corresponding to the chasing inside an SHC propagate outwards, see the movie in [19]. In figure 6(c), we plot the concentration of the dominant-item (temporary winner) at each site. Therefore, even though the items at sites in the purple region have apparently the same concentration close to 0.09 (that is, close to their global coexistence equilibria), there are small oscillations as a remnant of the entrainment. For the units closer to the pacemaker disc, the two hierarchy levels are manifest in the switching between different shades of green, blue, and red in the inner rings. A zoom into the time evolution of the nine items at the inner site (60, 60) (figure 6(d)), or slightly more remote from the pacemakers at site (55, 55) (figure 6(e)), reveals the sequences of nine different temporary winners, partially differing in their amplitudes. Outside the red rings (outer red ring in (b)) with a time evolution according to (f), the dynamics seems to have collapsed to a 1L-unit with a heteroclinic cycle between coexistence equilibria, but a closer look into the motion shows still nine oscillating items with very small amplitudes and even smaller differences between their amplitudes. At the outmost site at (5, 5) in the blue area of (b), the units seem to have approached a coexistence equilibrium (g), but also here remnants of a hierarchical heteroclinic motion are visible in a zoom. However, for the cases (f) and (g) the strong suppression of amplitudes and amplitude differences mean nevertheless a loss of information that was encoded in the winnerless competition sequence of information items closer to the disc of pacemakers.

Finally we would like to briefly mention another observation. So far it is preliminary and only based on numerical simulations. When we first let the target waves fully evolve and then quench the parameters of the pacemakers (2L-units) to parameters of CE-units, the whole grid of CE-units keeps on oscillating with hierarchical heteroclinic cycles for a rather long time. For example, if the quench is applied at time t = 1000 time steps as in the movie of [19], the waves, manifest in heteroclinic switching of items (colors), go on until at least t = 2500. This means that the system effectively has a rather high 'inertia'. For future work it seems interesting to pursue the very path in phase space that the system takes to the coexistence equilibrium after a quench in the bifurcation parameter. Moreover, our system should be complemented to account for the energetic cost of maintaining synchronized heteroclinic motion in comparison to other modes of synchronization, and to consider the cost as a function of the system parameters and the type and precision of the generated patterns.

We consider the pacemaker to be very close to its saddle (ρ/γ_P, 0, 0). Then we substitute s_1,1 with ρ/γ_P and s_1,2 = s_1,3 = 0 in equation (4) and get the reduced system of equations

$$\begin{align}\hfill {d}_{t}{s}_{2,1}& =\rho {s}_{2,1}-{\gamma }_{\mathrm{D}}{s}_{2,1}^{2}-{s}_{2,1}(c{s}_{2,2}+e{s}_{2,3})+\delta (\rho /{\gamma }_{\mathrm{P}}-{s}_{2,1}),\hfill \\ \hfill {d}_{t}{s}_{2,2}& =\rho {s}_{2,2}-{\gamma }_{\mathrm{D}}{s}_{2,2}^{2}-{s}_{2,2}(c{s}_{2,3}+e{s}_{2,1})-\delta {s}_{2,2},\hfill \\ \hfill {d}_{t}{s}_{2,3}& =\rho {s}_{2,3}-{\gamma }_{\mathrm{D}}{s}_{2,3}^{2}-{s}_{2,3}(c{s}_{2,1}+e{s}_{2,2})-\delta {s}_{2,3}.\hfill \end{align} \tag{ A.1 }$$

This driven system has three equilibria of interest:

• (a)  
  A one-item equilibrium: (S_a , 0, 0) which is a stable node for forcing $\delta > {\delta }_{{c}_{1}}$. In the original system, where the pacemaker is a heteroclinic cycle, the driven system has two additional equilibria, (0, S_a , 0), and (0, 0, S_a ), corresponding to the two other saddles of the pacemaker at (0, ρ/γ_P, 0), and (0, 0, ρ/γ_P), respectively.
• (b)  
  A two-items equilibrium: (S_b, S_c, 0). This equilibrium is a saddle equilibrium for $\delta > {\delta }_{{c}_{1}}$ and becomes a stable node for ${\delta }_{{c}_{2}}< \delta < {\delta }_{{c}_{1}}$. Again, when the pacemaker performs a heteroclinic cycle, there are two additional two-items equilibria: (0, S_b , S_c ), and (S_c , 0, S_b ).
• (c)  
  A three-items equilibrium: (S_d , S_e , S_f ). For $\delta > {\delta }_{{c}_{2}}$, this equilibrium is a saddle equilibrium, for ${\delta }_{{c}_{3}}< \delta < {\delta }_{{c}_{2}}$, the equilibrium becomes a stable node. For ${\delta }_{{c}_{4}}< \delta < {\delta }_{{c}_{3}}$, the equilibrium is a stable focus-node and for $\delta < {\delta }_{{c}_{4}}$, it is an unstable focus node. The unit driven by a pacemaker has two additional equilibria: (S_f , S_d , S_e ), and (S_e , S_f , S_d ).

The expressions of S_a , ..., S_f , ${\delta }_{{c}_{1}}$, and ${\delta }_{{c}_{2}}$ can be obtained analytically and are presented in [19]. The expressions for ${\delta }_{{c}_{3}}$ and ${\delta }_{{c}_{4}}$ are only numerically accessible, thus we have plotted their numerical values along with ${\delta }_{{c}_{1}}$ and ${\delta }_{{c}_{2}}$ in figure A1 for a set of parameters. In the following we restrict our discussion to the right of the vertical line in figure A1 with a unique correspondence of the equilibria between the critical lines, while multistability is observed to the left of this line that we do not further discuss. The vertical line corresponds to the choice of γ_D = γ_P.

Zoom In Zoom Out Reset image size

If we increase δ along a vertical line to the right of the dashed line in figure A1, we observe the following features. Below ${\delta }_{{c}_{4}}$, the driven unit has a stable limit cycle only if it is also in the regime of heteroclinic dynamics not considered further. If the driven unit is in the CE-regime, ${\delta }_{{c}_{4}}=0$. For ${\delta }_{{c}_{4}}< \delta < {\delta }_{{c}_{3}}$, the three-items coexistence equilibrium (3S) is stable, it is a focus-node and therefore the trajectory spirals towards the equilibria (as shown in figure A2(a) for δ = 0.05 and γ_D = 1.11). For ${\delta }_{{c}_{3}}< \delta < {\delta }_{{c}_{2}}$, the three-items equilibrium is still stable, but now it is a stable node (as shown in figure A2(b) for δ = 0.17). For ${\delta }_{{c}_{2}}< \delta < {\delta }_{{c}_{1}}$, the two-items equilibrium (2S) becomes a stable node (as shown in figure A2(c) for δ = 0.2). Finally, above ${\delta }_{{c}_{1}}$, the one-item equilibrium (1S) becomes a stable node and the two-items equilibrium loses stability (as shown in figure A2(d) for δ = 0.8). To the left of the vertical line, in the region bounded by ${\delta }_{{c}_{1}}$ from below and ${\delta }_{{c}_{5}}$ from above, two-items equilibria coexist with one-items equilibria.

Zoom In Zoom Out Reset image size

Now let us come back to our original system for which the pacemaker performs a heteroclinic cycle placed at the spatial index k = 1. As the heteroclinic trajectory goes from saddle to saddle, the driven unit (at location k = 2) goes from one of the above equilibria to the next equilibrium under the influence of the pacemaker. Let us slowly increase the strength of the heteroclinic forcing, δ, from zero on, and study its effect on the dynamics of the driven unit if γ_D is chosen to the right of the vertical line in figure A1.

${\delta }_{{c}_{4}}< \delta < {\delta }_{{c}_{3}}$ and γ _D = 1.11 . For figures A3(a)–(j), we consider the driven unit to be in the CE-regime by taking γ_D = 1.11. When the forcing δ is between ${\delta }_{{c}_{3}}$ and ${\delta }_{{c}_{4}}$, the three 3S-equilibria are focus nodes. When the trajectory of the pacemaker goes to the vicinity of one of its saddles, the trajectory of the driven unit spirals to the corresponding focus-node. As the trajectory of the pacemaker goes to another saddle, the trajectory of the driven unit also moves to the corresponding focus node and stays there until the pacemaker moves to its next saddle as shown in figure A3(a).

Zoom In Zoom Out Reset image size

${\delta }_{{c}_{3}}< \delta < {\delta }_{{c}_{2}}$. For this range, the system with constant forcing goes to the 3S-equilibrium (figure A2(b)), but the trajectory of the unit driven by a heteroclinic cycle tries to first visit the 1S- and 2S-equilibria, pushed by the pacemaker. However, since the 1S and 2S-equilibria are unstable (under constant forcing), if the driven unit has enough time while the pacemaker dwells on the corresponding saddle, it relaxes finally to one of the 3S-equilibria (which are stable as long as constant forcing is applied) (see figure A3(b)).

${\delta }_{{c}_{2}}< \delta < {\delta }_{{c}_{1}}$. Similarly, the trajectory of the driven unit visits the 1S-equilibria in this range before relaxing to the stable 2S-state in figure A3(c) ('stable' again referring to constant forcing), or it may also keep on residing in the unstable 1S-equilibrium as shown in figure A3(d) for δ = 0.4. However, if we add even a very small amount of noise to the dynamics of the driven unit, the trajectory relaxes to the 2S-equilibrium (not shown here).

$\delta > {\delta }_{{c}_{1}}$. Finally, for $\delta > {\delta }_{{c}_{1}}$, the trajectory of the driven unit goes to one of its 1S-stable nodes when the pacemaker is in the vicinity of the corresponding saddle (see figure A3(e)). So for large values of δ (above ${\delta }_{{c}_{1}}$), the pacemaker is able to keep the driven unit entrained to its dynamics.

In summary, based on the effective modeling of the driven unit under stepwise constant forcing, the driven unit directly coupled to the pacemaker approaches a 1S, 2S, or 3S-saddle equilibrium.

If we introduce a back-coupling δ_b to the pacemaker, a first effect is a slowdown of the slowing down of the heteroclinic oscillations. Thus, increasing δ_b is similar to increasing γ_P, thus effectively shifting the pacemaker closer to the coexistence regime with less pronounced slowing down. Beyond a certain strength of δ_b, the heteroclinic dynamics loses stability to limit cycle oscillations, and after a further increase, the limit cycles lose stability in a Hopf bifurcation, and the dynamics of both units stabilize to their respective coexistence equilibria. Further details are presented in [19].

Consider a driven unit. Depending on whether the trajectory of its respective nearest-neighbor pacemaker (its driving unit) is in the vicinity of a 1S, 2S, or 3S-equilibrium, the driven unit experiences forcing to either one, two, or all three of its items. As discussed before for a single driven unit, if the driven unit experiences a constant force at one of its items it has three types of equilibria: a 1S, 2S, or 3S-equilibrium. Here, in addition, when the driven unit at k ⩾ 3 is forced at two items at the same time, it has two equilibria: a 2S and a 3S. If it is forced at all three of its items, it has only one equilibrium of interest, a 3S.

As long as there is a pronounced slowing down of the heteroclinic motion, a driven unit is exposed for a long time interval to a constant forcing corresponding to the trajectory of its pacemaker approaching one of its equilibria and dwelling there, while for a comparatively very small time interval it is exposed to time-varying forcing when the pacemaker quickly transits from one of its equilibria to the next. Therefore the forcing experienced by the driven unit resembles a sequence of step functions with different constant forcing while different equilibria are visited by the driving unit, interrupted by jumps, when the driving unit transits between different equilibria. The transitions are visible in spikes of the trajectories, the more pronounced, the farther the distance from the original pacemaker at k = 1, as seen in figure B1.

Zoom In Zoom Out Reset image size

This figure shows the time series and phase-space trajectories of a system of five units unidirectionally coupled along a chain at sites k = 1 to k = 5. The time series of the three items of the kth unit are color-coded in red, green, blue. The items of the k − 1th unit, acting as a pacemaker to the kth unit, are shown in shades of gray. The amplitudes decrease with the distance to the pacemaker, as well as the differences in the dominance of different items. In the lower panel the phase space trajectories change from clear heteroclinic motion at k = 1 and k = 2 to increasingly more involved trajectories, exposed to more saddles and concentrated on a smaller part of phase space.

The possible stable or unstable equilibria are listed in table 2 of the main text. The approaches to these equilibria is determined by the following rules:

$$\begin{equation}1S\ \to \ 1S\ \text{or}\ 2S\ \text{or}\ 3S\ \text{or}\ (1S\to 2S)\ \text{or}\ (1S\to 2S\to 3S)\ \text{or}\ (2S\to 3S),\end{equation} \tag{ B.1a }$$

$$\begin{equation}2S\ \to \ 2S\ \text{or}\ 3S\ \text{or}\ (2S\to 3S),\end{equation} \tag{ B.1b }$$

$$\begin{equation}3S\to 3S.\end{equation} \tag{ B.1c }$$

For example, for unit-2 the rules should be read as follows: while unit-1 is at one of its three 1S-equilibria, unit-2 either dwells directly at one of its 1S, 2S, or 3S equilibria, or it moves to a more stable equilibrium according to the rules in the brackets of equation B.1(a). The options for unit-3 can be read off from equation (B.1) in the same way, if unit-2 is in the vicinity of 1S, or 2S, or 3S. Similarly for more remote distances from the pacemaker. Thus, the number of equilibria that are available to a unit to visit in its phase space increases, as the distance from the pacemaker at k = 1 increases.

The entrainment length corresponds to the number of units that get entrained to heteroclinic dynamics in a certain time window. For our numerical results, we consider the length of the chain to be N = 256 units. The results of varying δ are shown in figure B2. In panels figures B2(a)–(c), we plot the spatiotemporal plots of the first items s_k,1 of each unit for three values of the coupling strength δ = 0.5, δ = 0.6, and δ = 0.75, respectively. In the red region of the plots, the concentration of items s_k,1 is dominant at the location k, in the dark blue region it approaches zero while either s_k,2 or s_k,3 are dominant (winnerless competition). The entrainment process takes some transient time that is not negligible. The stronger the forcing δ, the more units synchronize to the heteroclinic motion. In figure B2(d), we plot the phase-space trajectories of a few units along the chain for δ = 0.5. While the pacemaker follows a heteroclinic cycle (blue line in (d)), a driven unit at position 50 wiggles on its approach of the heteroclinic cycle (red line) and even more so a unit at position 70 (yellow line). These features can be traced back to the proliferated saddle equilibria, felt by the respective trajectories. Here we also see that the effective modeling in terms of constant forcing to one, two, or three items loses its applicability, the more, the less pronounced the time scale separation between the dwell times in the vicinity of the saddles and the transits to the subsequent saddles is.

Zoom In Zoom Out Reset image size

Another parameter that affects the entrainment length is the bifurcation parameter γ_P of the pacemaker. Here, the entrainment length increases with increasing γ_P, this means, the closer the pacemaker is to the bifurcation towards the CE-regime. At a first view this may seem counterintuitive, as one may expect that the deeper the pacemaker is in the regime of heteroclinic dynamics, the stronger is its impact on the driven units. However, what prevents an easy entrainment with a pacemaker far off the bifurcation point is the fact that the heteroclinic motion slows down much faster for parameters deeply in the heteroclinic region.

The bifurcation parameter of the driven units γ_D also affects the entrainment length. Here, as expected, the entrainment length decreases with an increase in γ_D for which the driven units move deeper into the CE-regime and are harder for the pacemaker to entrain.

Moreover it should be mentioned that the entrainment length also depends on the initial conditions. The reason is that depending on the initial conditions, the transient time it takes the driving pacemaker to approach the vicinity of the heteroclinic trajectory differs and accordingly the time which is spent in the vicinity of the saddles. These differences propagate along the chain and get amplified, so that in a given time window the numbers of entrained units may slightly differ, even more so for a 2L-unit as a pacemaker, as there the initial conditions determine the number of revolutions in the SHCs.

If we close the unidirectionally coupled chain to a ring (figure 1(b)) such that only the Nth unit applies a weak forcing δ_P to the pacemaker (weaker in comparison to the unidirectional forcing δ between all other units), the entrainment depends on the size of the ring. Larger rings may not be synchronized at all to heteroclinic dynamics. As δ_P is increased and other parameter are kept fixed, the system shows quasi-periodic dynamics coexisting with heteroclinic dynamics, for further increased δ_P, limit cycles and finally a stable coexistence equilibrium, details are described in [19].

In [19] we discuss the role of bidirectional, but strongly asymmetric nearest neighbor coupling between a pacemaker and N − 1 driven units. From a physical point of view some back-coupling seems more realistic than unidirectional coupling. The forward coupling K_k,l = δ if k ∈ {2, ..., N} and l = k − 1, the back-coupling K_k,l = δ_b if k ∈ {1, ..., N − 1} and l = k + 1, K_k,l = 0 otherwise. The main difference is that each unit is now sensitive to the dynamics on both sides along the chain, the entrainment length becomes sensitive to the size of the chain. Initially more and more units get entrained, but after reaching a peak value over time in the number of entrained units, some recently entrained ones start de-entraining when they receive the feedback from the back-coupling, and then the entrainment length settles to a value less than the peak. The stronger the back-coupling, the smaller the entrainment length.