Emergence of target waves in paced populations of cyclically competing species

We investigate the emergence of target waves in a cyclic predator-prey model incorporating a periodic current of the three competing species in a small area situated at the center of the square lattice. The periodic current acts as a pacemaker, trying to impose its rhythm on the overall spatiotemporal evolution of the three species. We show that the pacemaker is able to nucleate target waves that eventually spread across the whole population, whereby three routes leading to this phenomenon can be distinguished depending on the mobility of the three species and the oscillation period of the localized current. First, target waves can emerge due to the synchronization between the periodic current and oscillations of the density of the three species on the spatial grid. The second route is similar to the first, the difference being that the synchronization sets in only intermittently. Finally, the third route towards target waves is realized when the frequency of the pacemaker is much higher than that characterizing the oscillations of the overall density of the three species. By considering mobility and the frequency of the current as variable parameters, we thus provide insights into the mechanisms of pattern formation resulting from the interplay between local and global dynamics in systems governed by cyclically competing species.


Introduction
Spatially extended excitable systems are investigated frequently due to their omnipresence in several biological, chemical and physical settings [1,2,3,4,5,6], as reviewed comprehensively in [7,8]. Especially pattern formation and the related self-organized emergence of target waves have received ample attention in the past [9]. Interestingly, although typically lacking an excitable character, systems describing competing interactions in natural and social systems may exhibit similar phenomena as well [10,11,12,13,14], including the formation of target waves [15,16,17,18,19,20]. Previous works have revealed that patterns in excitable systems emerge primarily due to the instabilities induced by the interplay between the fast excitatory and the slow recovery variable [21,22,23]. For example, the Bromous acid diffuses much faster than Ferroin in the Belousov-Zhabotinsky reaction, and the cyclic adenosine monophosphate (cAMP) diffuses much faster than membrane receptors in the Dictyostelium discoideum. However, many patterns, especially those occurring in mobile populations such as migrating animals or moving bacteria, cannot be explained by the mechanism that is applicable for the emergence of pattern formation in excitable systems. This is mainly because the excitability as a dynamical property is absent in systems where the moving of individuals makes all species equivalent in terms of their diffusion [24,25].
In the social amoeba Dictyostelium discoideum periodic currents of cAMP have been introduced in order to study the movement of the cells controlled by the emerging propagating waves of cAMP as during aggregation in the mound [26,27,28]. When Dictyostelium discoideum cells were left starving, they begun to emit pulses of cAMP which were relayed to more distant cells by radially divergent waves. Periodic injections of cAMP in slugs were found leading to the chemotactic attraction of anterior-like cells to the tip of the micropipette [26,27]. Dormann and Weijer reported that propagating chemoattractant waves are thus able to coordinate the periodic movement of the cells, and moreover, that synchronization between the frequency of cAMP injection and cell movement may occur [28]. While recently Cyrill et al. proposed a model that generates coherent target waves in recurrent single species populations [6], the mechanisms behind the formation of target waves within multi-species populations in the presence of an external periodic current have not yet been investigated.
The cyclic predator-prey model provides a concise description of competition in a multi-species environment. An experimental study by Kerr et al. [24] revealed that the mechanism of the rock-paper-scissors game can promote biodiversity of the three strains of Escherichia coli. Recently, Reichenbach et al. proposed a rock-paper-scissors game with mobile players [10,32], where the mobility was found to have a critical effect on species diversity. In particular, when mobility was below a critical value all species have been found to coexist by forming spiral waves, whereas above this threshold biodiversity was jeopardized. Motivated by the experimental studies on the three strains of Escherichia coli [24,25] and the species of Dictyostelium discoideum [26,27,28], we here investigate the pattern formation in a cyclic predator-prey model incorporating a localized periodic current of the three competing species. The periodic current essentially acts as a pacemaker on the population, trying to impose its rhythm on the spatiotemporal evolution of the species, but it can also be regarded as a localized inhomogeneity affecting the overall system dynamics [29]. We report the emergence of coherent target waves, in general emerging due to the nucleation induced by the localized current. More precisely, three routes to target wave formation can be distinguished based on the interplay between local and global dynamics. Local thereby refers to the properties associated with the oscillations within the paced region of the spatial grid whereas global refers to the overall dynamics of the three species. In particular, target waves can emerge either due to the constant or due to the intermittent synchronization between the local and global dynamics, or when the frequency of the pacemaker, reflecting directly the so-called local dynamics, is much higher than the frequency characterizing the oscillations of the overall density of the three species, i.e. the global dynamics. Interestingly, although the presently employed model is based on microscopic interactions among individuals only, the target waves emerge on a macroscopic level. Furthermore, our work indicates possibilities of pattern control and selection in systems governed by cyclical interactions.

Three-species cyclic predator-prey model
Based on previous works of Reichenbach et al. [10,31,32], we employ a three-species cyclic predator-prey model with the following specifications. Nodes of a L × L square lattice present mobile individuals belonging to one of the three species, which we denote by A, B and C. Each node can either host one individual of a given species or it can be vacant. Vacant sites, which we denote by ⊗, are also the so-called resource sites. Within the model three processes are possible; namely predation, reproduction and exchange, whereby these occur only between neighboring nodes. Predation.-Species A eliminates species B at a rate 1, whereby the node previously hosting species B becomes vacant. In the same manner species B can eliminate species C, and species C can eliminate species A, thus forming a closed loop of dominance between them. Reproduction.-Individuals can place an offspring to a neighboring vacant node ⊗ at a rate 1. Exchange.-Two individuals, including vacant sites, can exchange their position at a rate α, thus introducing mobility of participants. According to this description, predation (upper row), reproduction (middle row) and exchange (lower row) can be described by the reactions: where X, Y ∈ {A, B, C, ⊗}. According to the random walk theory [35], the mobility of individuals M can be defined in terms of the exchange process as M = 2α/L 2 , meaning it is proportional to the typical area explored by a mobile individual per unit time. Here the radius is equal to R = 3.5. The blue area is subjected to a periodic injection of the three species, in the recurrent order A → B → C → A, each every multiple of T 0 (see also main text for details). The gray area is initially populated either by A, B, C and ⊗ with equal probability [initial condition 1 (IC1)] or solely with vacant sites [initial condition 2 (IC2)]. Throughout this paper we will use R = 10.5 and L such that 0.01 ≤ R/L ≤ 0.02 unless stated otherwise. We note that the results reported below are robust to variations of this ratio as long as the pacing area (blue) is large enough to nucleate target waves by appropriate M and T 0 .
Unlike the deterministic approach, which regards the time evolution as a continuous process, here we make use of a stochastic simulation algorithm whereby the temporal evolution can be considered as a random walk process. The most commonly applied stochastic simulation algorithm was developed by Gillespie [33,34], where reactions occur in a random manner. In particular, predation and reproduction occur with probability 1/(α + 2), whereas exchange (moving) occurs with probability α/(α + 2). If one ignores the spatial structure and assumes the system to be well mixed, the model can be described by partial differential equations. The results presented below, however, are obtained via Monte Carlo simulations of the L × L square lattice, whereby we use no-flux boundary conditions and the algorithm of Gillespie [33,34] for determining the probabilities for the three possible processes. An elementary Monte Carlo step consists of randomly choosing an individual who interacts with one of its four nearest neighbors, which is also selected randomly, and then executing the process as determined by the Gillespie's algorithm. One full Monte Carlo step consists of N = L 2 elementary steps, during which, in accordance with the random sequential update, each player is selected once on average. The resolution of time is thus measured in full Monte Carlo steps.
Furthermore, motivated by the application of periodic currents of cAMP in the experiments with Dictyostelium discoideum [26,27,28], we introduce a localized periodic current of the three competing species, as described and schematically presented in Fig. 1. In particular, the periodic current is applied over a small (compared to the overall system size) area located at the center of the spatial grid. The periodic current thus acts as a pacemaker on the population, trying to impose its rhythm on the spatiotemporal evolution of the three species. The current is defined as follows: at time t = 0 the nodes inside R (see Fig. 1) are populated by species A, at time t = T 0 these nodes are populated by species B, at time t = 2T 0 these nodes are populated by species C, and at time t = 3T 0 these nodes are again populated by species A, and continuing further in this manner. Importantly, during nT 0 < t < mT 0 , where n ∈ Z + and m = n + 1, the evolution of species inside R (blue area in Fig. 1) is governed by the same Monte Carlo updating as outside (gray area in Fig. 1). Moreover, from the definition of the current it follows that its period equals T in = 3T 0 . In the following, we will consider the mobility M and the time T 0 between successive replacements of a species inside R as the two crucial parameters effecting the emergence of target waves in the examined model.
Before we start presenting the results, we note that in the absence of the periodic current the model returns identical results in terms of the coexistence of species in dependence on M as reported earlier by Reichenbach et al. [10]. In particular, while low M ensure coexistence of all three species, there exists a critical M ≈ 10 −3 beyond which diversity is no longer sustained.

Emergence of coherent target waves
We start by presenting snapshots of the spatial grid obtained at different times t in Fig. 2(a). It is evident that the initial random configuration (t = 0) is gradually replaced by supremely ordered target waves (t = 600T 0 ), which emerge after a transient period that is characterized by disordered turbulence at intermediate times (t = 200T 0 and t = 400T 0 ). The emergence of target waves depicted in the rightmost panel of Fig. 2(a) is a direct consequence of the localized periodic current, i.e. the pacemaker, which periodically introduces an inhomogeneous nucleus at the center of the spatial grid. These inhomogeneities, in form of regions dominated by a single species, spread and eventually self-organizes into a beautiful manifestation of pattern formation in a population governed by cyclical interactions. A more precise temporal evolution of target wave formation in the system can be visualized by means of the corresponding space-time plot, as is shown in Fig. 2(b). From the latter it is indeed evident that the target waves begin to grow around the impact site of the periodic current (centered at X = Y = 250; note that L = 500), which periodically introduces inhomogeneous regions governed by a single species. Following a period of disordered patters and very turbulent waves (t < t 1 ), as well as an intermediate period during which target waves gradually begin to dominate (t 1 ≤ t ≤ t 2 ), the complete dominance of ordered target waves finally prevails (t > t 2 ) in a resonance-like manner due to the interplay between the localized periodic current and the overall spatiotemporal dynamics of the three species. It is worth noting that the two times t 1 and t 2 in Fig. 2(b), denoting the start of target waves formation and their complete dominance, respectively, depend somewhat To ensure clarity some frames within t 1 ≤ t ≤ t 2 were discarded; t 1 and t 2 corresponding to 33T 0 and 466T 0 , respectively. In both panels yellow, red and blue squares denote species A, B and C, respectively, whereas the grey squares depict vacant sites. Employed parameter values are: M = 10 −4 , T 0 = 150 and L = 500. Initial conditions: IC1 (see the caption of Fig. 1). The full spatiotemporal evolution of the system is recorded in the file video2.avi, which is part of the supplementary material pertaining to this paper.
on the random realization of initial conditions and the location of the periodic current. However, such deviations are negligible and do not affect our results vitally.
Further elaborating on the emergence of target waves depicted in Fig. 2, we present in Fig. 3(a) oscillations of the overall density of the three species after coherent target waves start dominating the spatial grid. In accordance with the spatiotemporal dynamics of target waves, the temporal outlays of the density of all three species are predominantly periodic, which manifests as sharp peaks in the corresponding Fourier transformation of the series, occurring at multiples of the main oscillation frequency. It is more instructive, however, to compare the oscillations of the density of a given species measured only within the paced (injection) area and its overall density on the spatial grid. Figure 3(b) features these two temporal plots for species A, and it can be observed that the two are in synchrony. This confirms that the target waves emerge as a consequence of the interplay between local (encompassing only the paced area) and global (encompassing the whole spatial grid) system dynamics, which during the dominance of ordered target waves [see rightmost panel of Fig. 2(a)] can be synchronized.
Nevertheless, above results pertain to a single combination of M and T 0 , and thus it is of interest to investigate the impact of other values of these two parameters on pattern formation as well. Figure 4(a) shows typical snapshots of the spatial grid obtained for different mobility M of individuals after a long simulation time. Evidently, low mobility fails to evoke target patterns. In the leftmost panel of Fig. 4(a), where M = 10 −5 , turbulent spirals dominate over the entire lattice. For M = 10 −4 , on the other hand, the locally introduced inhomogeneities due to the periodic current are appropriately enhanced to eventually result in the emergence of target waves, as depicted in the middle panel of Fig. 4(a). Increasing the mobility further to M = 10 −3 evokes a transition from target to spiral waves, as depicted in the rightmost panel of Fig. 4(a), which is due to the strong mixing of individuals prohibiting the stability of the relatively stationary target wave pattern and instead favoring the more dynamically evolving spiral waves.  Fig. 4(b), which results from a systematic study of pattern formation over the whole 2D parameter plane. Depending on the combination of M and T 0 , three different evolutionary outcomes can be distinguished. Namely, the formation of coherent target waves (region T), the formation of turbulent or ordered spiral waves (region S), or the formation of a uniform phase (region U) following extinction of two species. From the phase diagram presented in Fig. 4(b) it follows that the region of coherent target waves is relatively small, and that thus the phenomenon results from a rather subtle interplay between the localized periodic current and the overall dynamics of the three species. It is worth noting that these results are largely independent of the system size as long as the ratio R/L remains the same, which ensures that the local dynamics induced by the periodic current is sufficiently represented in the system. Obviously, increasing the value of R by a given system size L facilitates the formation of target waves, whereas smaller R are less likely to have a noticeable impact on the overall evolution of the three species. In general, however, quantitative differences due to such variations are small, affecting the borders of the hatched region (T) [ Fig. 4(b)] in the M − T 0 plane slightly. We also note that the border separating turbulent and ordered spiral waves in the white region (S) is located at M = 6 · 10 −5 . We have observed that smaller values of M lead to turbulent spiral waves [see the leftmost panel of Fig. 4(a)], while larger values of M within the S region tend to yield ordered spiral waves [see the rightmost panel of Fig. 4(a)].

Three routes towards target waves
In order to develop a better understanding of the emergence of coherent target waves evoked by the localized periodic current for different values of M and T 0 , we change the initial setup in that we use IC2 type initial conditions, i.e. apart from the paced area, the lattice is initially populated solely with vacant sites (see also the caption of Fig. 1). This is inspired by the growth experiments of Escherichia coli [16,17], and enables us a clearer tracking of the localized nucleus formation and subsequent spreading of the periodically injected species across the spatial grid. In the following, we will show that indeed three routes to the emergence of coherent target waves can be distinguished.  Fig. 3(b), for M = 10 −3 the oscillation period of the overall density of species A on the spatial grid is much larger than that of the density within the paced area. Thus, the two cases constitute two different routes towards target waves [see also Figs. 6(b) and (d) for the corresponding snapshots] in the examined system. The third route emerges if M = 10 −4 , where the synchronization between local and overall oscillations occurs only intermittently. As a consequence, the resulting target waves do not have a well-defined wavelength, as can be observed in Fig. 6(c). Figure 5(k) depicts succinctly the intervals of M (for the considered T 0 = 600) resulting in a given route to coherent target formation. Notably, the gray region of M (aperiodic regime) results only in disordered turbulent spiral-like waves [see also Fig. 6(a) for the corresponding snapshot].
From the results presented in Fig. 5, we thus conclude that there exist three different routes to coherent target waves in the examined cyclic predator-prey model incorporating a periodic current.
First, target waves can emerge due to the synchronization between the periodic current and oscillations of the density of the three species on the spatial grid, as presented in Fig. 5(i). The resulting spatial portrait is shown in Fig. 6(b). Due to the regularity of the oscillations the target waves also have an easily inferable wavelength λ = D/L, as noted on the corresponding snapshot. The second route is similar to the first, the difference being that the synchronization sets in only intermittently. The resulting spatial portrait is shown in Fig. 6(c), whereby the seeming lack of a well-defined wavelength is a direct consequence of the intermittent temporal outlay of the variations of the species' densities, as shown in Figs. 5(c) and (g). Note that in Fig. 6(c) the circular ring of the blue species, for example, is sometimes thinner and sometimes thicker, i.e. not always the same. This is because the oscillation period of the underlying temporal trace is not exact but varies intermittently. Finally, the third route towards target waves is realized when the frequency of the pacemaker is much higher than that characterizing the oscillations of the overall density of the three species, as presented in Fig. 5(j). The resulting spatial portrait is shown in Fig. 6(d is also worth noting that for sufficiently large M, exceeding 5.0 × 10 −4 , the mobility quickly decimates any given species that is injected within the paced area, and thus the artificially established dominance there is very short-lived [note the sharp peaks in Fig. 5(d)]. In fact, individuals of an injected species always aggregate near the boundary of the paced area, and only then start to spread outwards across the spatial grid. But since this aggregation and propagation takes a long time, the period of overall oscillations is much longer than that of oscillations within the paced area. On the other hand, for extremely small values of M (< 3.0 × 10 −6 ), individuals spreading out of the area that is subjected to the periodic current distribute almost homogeneously across the spatial grid due to the relatively high rate of predation and reproduction if compared to the exchange rate (mobility). Thus, the injected individuals become fragmented, in turn failing to evoke ordered target waves, as depicted in Fig. 6(a). Accordingly also, the oscillations of the species' densities, both within the paced area and overall, are aperiodic [see Figs. 5(a) and (e)].
Focusing further on the wavelength λ of the target waves depicted in Fig. 6 (where applicable), one can thus observe that the mobility very much affects the typical spatial distance between neighboring wave fronts with the same species, and while the first [panel (b)] and the third route [panel (d)] yield target waves with an easily inferable λ, this is not necessarily the case for the second route [panel (c)]. In fact, results presented in Fig. 7(a) show that the wavelength of target waves displays non-monotonous behavior in dependence on the mobility of individuals for T 0 = 600 and T 0 = 400, while increases with M according to λ ∼ log M 1/2 for T 0 = 300. Notably, this is different from what has been reported for spiral waves [10], where λ ∼ M 1/2 has been found applicable. Furthermore, results in Fig. 7(a) suggest that indeed both, M and T 0 affect the route towards the emergence of target waves. Figure 7(b) presents results of an extensive analysis performed over the M − T 0 parameter space, in particular indicating the properties of the oscillations of species' densities, similarly as in Figure 5(k). Gray color denotes that, for the particular combination of M and T 0 values, the emergence of coherent target waves is impossible due to an aperiodic outlay of the oscillations, both within the paced area as well as overall [see Figs. 5(a) and (e) for an example]. Yellow color indicates the emergence of coherent target waves via the first route, which is due to the synchronization between the periodic current and oscillations of the overall density of the three species on the spatial grid [see  Fig. 7(b) it is also clear that the non-monotonous dependence of λ on mobility for T 0 = 400 and T 0 = 600 depicted in Fig. 7(a) is due to the emergence of target waves via the second route for intermediate M (red stripes). Finally, we note that when the route towards the emergence of target waves is switched, narrow parameter regions of spiral wave formation may appear [see also Fig. 4(b)], which, interestingly, were also observed in the Dictyostelium discoideum studies [6,37].

Summary
In sum, we have studied the emergence of target waves in a cyclic predator-prey model incorporating a periodic current of the three competing species in a small area situated at the center of the square lattice. We have shown that a pacemaker-like periodic current may evoke coherent target waves, provided the mobility of individuals on the spatial grid and the frequency of the forcing are adequately adjusted. Furthermore, we have identified three possible routes towards the emergence of coherent target waves, depending on the properties of the oscillations of species' densities within the paced area and across the whole spatial grid. In particular, target waves can emerge due to the complete or intermittent synchronization between the periodic current and oscillations of the overall densities of the three species, or they can emerge so that the frequency of the pacemaker is much higher than the frequency characterizing the induced global oscillations of the three densities. Irrespective of the scenario, however, the global oscillations and the associated target waves are a direct consequence of the introduced pacemaker. Which route is chosen depends both on the mobility of the individuals as well as on the selected oscillation period of the localized current. We thus provide insights into the mechanisms of pattern formation resulting due to the interplay between local and global dynamics in a system that is governed by cyclically competing species. Since the identified mechanisms affect the spatiotemporal dynamics in a predictable way, and moreover, can be tuned effectively via accessible system parameters like the frequency of the pacemaker, our work indicates possibilities of pattern control and selection in systems governed by cyclical interactions. Cyclical competitions evidently play an important role by the maintenance of ecological biodiversity and emergence of cooperation in structured populations [38,39,40,41,42,43]. Related to this, complex spatiotemporal formations have often been observed, as for example spiral wave structures in the realm of the rock-paper-scissors game [44]. Coherent target waves, however, are rarely reported in this context. Although target waves were observed in excitable systems, the main difference of the present model is that the diffusion induced by individual mobility is the same for all species. Therefore, the emergence of target waves is due to a different mechanism. Notably, the periodic injection method has recently also been used in a complex Ginzburg-Landau system for the purpose of spatiotemporal chaos control [29,45]. According to our study, one can control the route towards target wave formation by adjusting the frequency of the periodic current. In this way also the wavelength of target waves can be controlled, which may find useful applications in a variety of realistic systems, ranging from targeted drug delivery to sustenance of biodiversity. We hope that our study will promote the understanding of social dynamics and strengthen the established importance of physics [46] when striving towards this goal.