The Influence of Mobility Rate on Spiral Waves in Spatial Rock-Paper-Scissors Games

We consider a two-dimensional model of three species in rock-paper-scissors competition and study the self-organisation of the population into fascinating spiraling patterns. Within our individual-based metapopulation formulation, the population composition changes due to cyclic dominance (dominance-removal and dominance-replacement), mutations, and pair-exchange of neighboring individuals. Here, we study the influence of mobility on the emerging patterns and investigate when the pair-exchange rate is responsible for spiral waves to become elusive in stochastic lattice simulations. In particular, we show that the spiral waves predicted by the system's deterministic partial equations are found in lattice simulations only within a finite range of the mobility rate. We also report that in the absence of mutations and dominance-replacement, the resulting spiraling patterns are subject to convective instability and far-field breakup at low mobility rate. Possible applications of these resolution and far-field breakup phenomena are discussed.


I. INTRODUCTION
Understanding the mechanisms allowing the maintenance of species coexistence is an issue of paramount importance [1]. Evolutionary game theory [2][3][4], where the success of one species depends on what the others are doing, provides a fruitful framework to investigate this question by means of paradigmatic schematic models. In this context, cyclic dominance is considered as a possible motif enhancing the maintenance of biodiversity, and models of populations in cyclic competition have recently received significant attention.
The rock-paper-scissors (RPS) game -in which rock crushes scissors, scissors cut paper, and paper wraps rock -and its variants are paradigmatic models for the cyclic competition between three species. Examples of RPS-like dynamical systems can be Uta stansburiana lizards, and communities of E.coli [5,6,8,9], as well as coral reef invertebrates [10]. In the absence of spatial degrees of freedom and mutations, the presence of demographic fluctuations in finite populations leads to the loss of biodiversity with the extinction of two species in a finite time, see e.g. [11][12][13][14][15][16]. However, in nature, organisms typically interact with a finite number of individuals in their neighborhood and are able to migrate. It is by now well established both theoretically and experimentally that space and mobility greatly influence how species evolve and how ecosystems self-organize, see e.g. [17][18][19][20][21][22]. Of particular relevance are the in vitro experiments with Escherichia coli of Refs. [5][6][7][8] showing that, when arranged on a Petri dish, three strains of bacteria in cyclic competition coexist for a long time while two of the species go extinct when the interactions take place in well-shaken flasks. On the other hand, in vivo experiments of Ref. [23] with bacterial colonies in the intestines of co-caged mice have shown that mobility allows the bacteria to migrate between mice and to maintain their coexistence. These observations illustrate that mobility can both promote and jeopardize biodiversity in RPS games, as argued in Refs. [24][25][26]: In the experiments of Ref. [23] biodiversity is maintained only when there is migration, whereas in Ref. [6] species coexistence is lost in well-shaken flasks corresponding to a setting with a high mobility rate.
These considerations have motivated a series of studies aiming at investigating the relevance of spatial structure and individual's mobility on the properties of RPS-like systems. For instance, various two-dimensional versions of the model introduced by May and Leonard [27], characterized by cyclic "dominance removal" in which each species "removes" another in turn (see below), have received much attention [24][25][26][28][29][30][31]. It has been shown that species coexist for a long time in these models with pair-exchange among neighboring individuals: below a certain mobility threshold species coexist by forming intriguing spiraling patterns below a certain mobility threshold, whereas there is loss of biodiversity above that threshold [24][25][26]. Another popular class of RPS models are those characterized by a zero-sum cyclic interactions (via "dominance replacement") with a conservation law at mean-field level ("zero-sum" games) [32][33][34][35][36][37][38][39][40] and whose dynamics in two-dimensions also leads to a long-lasting coexistence of the species but not to the formation of spiraling patterns, see e.g. [28,31,37]. Recent studies, see e.g. [41][42][43][44][45][46][47], have investigated the dynamics of two-dimensional RPS models combining cyclic dominance removal and replacement, while various generalization to the case of more than three species have also been considered, see e.g. [48][49][50]. In Refs. [43][44][45][46][47] we studied the spatio-temporal properties of a generic twodimensional RPS-like model accounting for cyclic competition with dominance-removal and dominance-replacement, along with other evolutionary processes such as reproduction, mutation and mobility via hopping and pair-exchange between nearest neighbors. By adopting a metapopulation formulation and using a multiscale and size-expansion analysis, combined with numerical simulations, we analyzed the properties of the emerging spatio-temporal dynamics. In particular, we derived the system's phase diagram and characterized the spiraling patterns in each of the phases and showed how non-linear mobility can cause the far-field breakup of spiral waves. In spite of the predictions of the theoretical models, it is still unclear under which circumstances microbial communities in cyclic competition would self-arrange into spiraling patterns as those observed in other systems such as myxobacteria and slime molds [51,52].
Here, we continue our investigation of the generic two-dimensional RPS-like model of Refs. [43][44][45] by focusing on the influence of pair-exchange between nearest-neighbors, as simplest form of migration, on the formation of spiraling patterns in two-dimensional lattice simulations. In particular, we demonstrate a resolution issue: on a finite grid spiral waves can be observed only when the migration is within a certain range. We also show that in the absence of mutations and dominance-replacement, e.g. as in Refs. [24][25][26]41], the spiraling patterns emerging from the dynamics are subject to convective instability and far-field breakup at low mobility rate. This paper is structured as follows: The generic metapopulation model [44,53] is introduced in Sec. 2 and, building on Refs. [43,45], the main features of its description at The composition of a patch evolves in time according to the processes (1) and (2). Furthermore, migration from the focal patch (dark gray) to its four nearest-neighbors (light gray) occurs according to the processes (3), see text. Adapted from [45]. mean-field level and in terms of partial differential equations (PDEs) are outlined. The section 3 is dedicated to a summary of the characterization of the system's spiraling patterns in terms of the underlying complex Ginzburg-Landau equation (CGLE). Section 4 is dedicated to our novel results concerning the resolution issues in finite lattices and the far-field breakup and convective instability under low mobility. Finally, we conclude by summarizing and discussing our findings.

II. MODEL
As in Refs. [43,45], we consider the generic model of cyclic dominance between three competing species defined on an L × L periodic square lattice of patches, L being the linear size of the grid, where each node of the grid is labelled by a vector = ( 1 , 2 ). As illustrated in Fig. 1, each patch consists of a well-mixed population of species S 1 , S 2 , S 3 and empty spaces ∅ and has a limited carrying capacity N : In each patch there are therefore at most N individuals N S i ( ) of species S i (i = 1, 2, 3), and there are also N ∅ ( ) = N −N S 1 ( )−N S 2 ( )−N S 3 ( ) empty spaces. Within each patch , the population composition evolves according to the most generic form of cyclic RPS according to the following schematic reactions: where the index i ∈ {1, 2, 3} is ordered cyclically such that S 3+1 ≡ S 1 and S 1−1 ≡ S 3 . The reactions (1) describe the generic form of cyclic competition, that comprises dominanceremoval with rate σ and dominance-replacement with rate ζ [43,45,46]. The processes (2) allow for the reproduction of each species (with rate β) independently of the cyclic interaction provided that free space (∅) is available within the patch. The biological interpretation of the mutations S i − → S i±1 (with rate µ) is, e.g., that they mimic the fact that sideblotched lizards Uta stansburiana undergo throat-color transformations [9], while from a mathematical perspective they yield a supercritical Hopf bifurcation at mean-field level, see Refs. [54,55], about which a multiscale expansion is feasible, see below and Refs. [43,45,46].
Since we are interested in analyzing the spatio-temporal arrangement of the populations, in addition to the intra-patch reactions (1)-(2), we also allow individuals to migrate between neighboring patches and via pair exchange, according to where At an individual-based level, the model is defined by the Markov processes associated with the reactions (1)-(3). The model's dynamics is thus described by the underlying master equation, and the stochastic lattice simulations performed using the Gillespie algorithm [56], as explained in Refs. [43,45,47]. The metapopulation formulation of the model makes it well-suited for a size expansion of the master equation in 1/N [45,47,[57][58][59][60]. Such an expansion in the inverse of the carrying capacity has been detailed in Ref. [45] where we showed that to lowest order in the continuum limit (L 1) on a square domain of size S ×S the system evolves according to the partial differential equations (with periodic boundary conditions) for the species densities and ρ = s 1 + s 2 + s 3 . Here and henceforth, without loss of generality we have rescaled time by setting β = 1, and on average there are N microscopic interactions during a unit of time [45]. As usual, the diffusion coefficient D and the migration rate δ are simply related It is worth noting that contrary to the models of Refs. [43,45], here we do not consider nonlinear diffusion: Movement in (4) appears only through the linear diffusive terms D∆s i .

III. SPIRALING PATTERNS AND THE COMPLEX GINZBURG-LANDAU EQUATION
In Refs. [43,45] we showed that the dynamics in terms of the PDEs (4) yield spiraling patterns whose spatio-temporal properties can be analyzed in terms of the system underlying complex Ginzburg-Landau equation (CGLE) [62]. The latter is derived by introducing the "slow variables" (X, T ) = ( x, 2 t), where = 3(µ H − µ) is the system's small parameter in terms of which a multiscale expansion is performed about the HB [61]. Details of the derivation can be found in Ref. [45] and brief accounts in Refs. [43,46]. Here we quote the system's CGLE for the complex modulated amplitude A(X, T ) which is a linear combination of the rescaled species densities [43,45,46]: where ∆ X = ∂ 2 and after having rescaled A by a constant, we find the parameter As explained in Refs. [43,45], the CGLE (5) allows us to accurately characterize the spatio-temporal spiraling patterns in the vicinity of the HB (for 1 i.e. µ µ H ) by using the well-known phase diagram of the two-dimensional CGLE, see e.g. [62], and to gain significant insight into the system's spatio-temporal behavior away from the HB (we here In all panels, the initial condition is a random perturbation of the homogeneous state s * , see text and [44]. Adapted from [47]. even away from the HB whose boundaries are essentially the same as in the vicinity of the HB, see Fig. 3. At low mutation rate, there is no spiral annihilation and the SA phase is generally replaced by an extended BS phase (with far-field breakup of the spiral waves when σ ζ).
Away from the cores of the spiral waves, the solution to the CGLE (5) can be approximated by the travelling-wave ansatz A(X, T ) = Re i(k·X−ωT ) of amplitude R(c), angular frequency ω and wave number k [45]. As detailed in Ref. [45], the wavelength λ = 2π/( k) of the spiraling patterns in the BS and EI phases near the HB in the physical space (in lattice units) can thus be estimated by λ ≈ λ H , where where the amplitude R 2 (c) is a decreasing function of c in the BS and EI phases (see Fig. 6 in [45]) and has been found numerically to be in the range 0.84-0.95 in the BS phase [45].
In Ref. [45], we showed that the description in terms of the CGLE is not only valid and accurate near the HB (µ µ H ), but it is also insightful in the limit µ µ H of low mutation rate. In fact, we showed that lowering the mutation rate µ below µ H results in three regimes: the AI, EI and BS phases, see Fig. 3. We also found that reducing µ below µ H results into shortening the wavelength λ of the spiraling patterns in the BS and EI phases: The wavelength at low mutation rate satisfies a linear relationship (see Fig. 14 in Ref. [45]): where the wavelength λ 0 at µ = 0 is inferred from the numerical solution of the PDEs (4) and is shorter than the wavelength λ µ H near µ H ; typically λ 0 ∈ [0.3λ µ H , 0.5λ µ H ] and it scales as λ 0 ∼ 2δ(3 + σ)/σ [63]. For instance, when (β, σ, ζ, µ, δ) = (1, 1, 0.6, 0.01, 0.64), we have ≈ 0.308, R 2 ≈ 0.9 and (7) yields λ H ≈ 52 while we found λ 0 ≈ 26 and therefore at µ = 0.01, the wavelength is λ µ=0.01 ≈ 32 which is in good agreement with lattice simulations, see Fig. 4 (bottom, left).

RAL WAVES ON A GRID?
We have seen that near the HB the RPS dynamics is generally well described in terms of the PDEs (4) and CGLE (5) when N 1. Accordingly, the effect of space and individuals' mobility is accounted by linear diffusion. In this setting, rescaling the mobility rate δ → αδ (α > 0) boils down to rescale the diffusion coefficient and space according to D → αD and x → x/ √ α. Hence, the size of the spatial patterns increases when the individuals' mobility is increased, and it decreases when the mobility is reduced [24][25][26]45]. Furthermore, according to the description in terms of the CGLE, the mobility rate and diffusion coefficient do not affect the stability of the spiraling patterns near the HB but only change their size. Below, we discuss the cases of very small or large mobility rate and show that this may result in spiral waves being elusive and/or unstable even in regimes where the PDEs (4) and CGLE (5) predict that they would exist and be stable.

A. Resolution issues
In this section, we focus on resolution issues and show that while the description of the dynamics in terms of the PDEs (4) predict the formation of spiral waves these cannot be observed on finite lattice due to resolution issues. In other words, we report that, even when the PDEs (4) and CGLE (5) predict that the dynamics leads to stable spiraling patterns (BS phase), these may be elusive when the mobility rate is too low or too high as illustrated in Figure 4. In order to determine the range of the mobility rate δ within which stable spiral waves can be observed on a two-dimensional grid, we distinguish three regimes, see Table I: In regime (i), the PDEs (4) predict a myriad of tiny spirals of wavelength of the order of one unit lattice space. Clearly, the resolution of any finite lattice is insufficient to allow to observe spiraling patterns of such a tiny size (of order of one pixel) on the grid, see Figure 4 (top). In this situation, instead of spiral waves lattice simulations lead to apparently clumps of activity. As shown in Figure 4, this phenomenon does not stem from demographic noise since it is present even when N is very large, as shown in Figure 4 (where N = 1024). Yet, due to their small size, these emerging incoherent spatio-temporal structures are prone to be affected by demographic fluctuations and result in being noisy. In regime (ii), the spiral waves' wavelengths, λ H (7) near the HB and λ µ µ H (8) at low mutation rate, are much larger than the inter-patch space and smaller than the domain size. Hence, stable spiraling patterns fit within the lattice and are similar to those predicted by the PDEs (4), see Figure 4 (bottom left). In regime (iii), the predicted λ H and λ µ µ H outgrow the lattice and the arms of the resulting large spirals appear like planar waves, see Figure 4 (bottom right). Interestingly, planar waves have been found in the model of Ref. [42] without mutations (µ = 0) at sufficiently high pair-exchange rate.
In order to estimate the boundaries between these regimes, we can use the relations (7) and (8)   Wavelength λ Mobility Rate δ Diffusion Coefficient D (S = 1) Patterns on Grid

Planar waves
This means that stable spirals of wavelength given by (7) or (8) can be observed in the BS phase in the range of mobility rate 2 δ (L /κ) 2 that grows with L. Hence, when L is sufficiently large lattice simulations will lead to the formation of stable spiral waves almost for any finite mobility rate. However, as the size of plates used in most microbial experiments rarely exceeds L = 100 [6], it is interesting to consider the case where L is not too large. In particular, when the ratio L/κ = O(1), the range 2 δ (L /κ) 2 is finite and spiral waves outgrow the lattice even for a finite mobility rate δ (L /κ) 2 . For instance, in Figure 4, we have L = 64 and (β, σ, ζ, µ, δ, N, L) = (1, 1, 0.6, 0.01, 1024) which yields ≈ 0.3, c ≈ 1.0, R 2 ≈ 0.9 [45,47], and L /κ ≈ 0.92. In this example, with (7) and (8) we find λ(δ = 0.64) = L/2 and λ(δ = 1.44) = 3L/4. According to the above discussion we expect to find visible spiral waves for δ = 0.64 and patterns resembling planar waves when δ = 1.44 which is confirmed by the bottom row of Figure 4. For the example of Figure 4, the PDEs (4) predict the formation of small spiral waves of wavelength λ = 1 and λ = 3 for δ = 0.000625 and δ = 0.005625 respectively, which result in the noisy clumps of activities on the lattice of the top panel of Figure 4.

B. Far-field breakup of spiral waves under weak pair-exchange rate
Variants of the two-dimensional RPS model (1)-(3) without mutation (µ = 0) have received significant interest and many authors have studied under which circumstances the dynamics leads to the formation of stable spiraling patterns, see, e.g., [24-26, 28, 30, 32, 35-38, 41]. In Ref. [24][25][26], where only the dominance-removal was considered, it was found that the cyclically competing populations moving under pair-exchange always form persisting spiraling patterns under a critical mobility threshold whereas no such coherent patterns were found in a similar system where the cyclic competition was implemented according to the dominance-replacement process, see e.g. [28,31,37]. By means of an approximate mapping onto a CGLE, the authors of Ref. [28] concluded that while the model with dominanceremoval (σ > 0, ζ = 0) can sustain spiral waves, this is not the case of models with dominance-replacement (σ = 0, ζ > 0). In Ref. [41], it is found that the combination dominance-removal and dominance-replacement processes can lead to stable spiraling patterns as well as to convectively unstable spiral waves. This picture was complemented and unified in our recent work [43,[45][46][47] where these questions were considered in the presence/absence of a small mutation rate and nonlinear mobility (pair-exchange and hopping processes were divorced). In particular, we showed that when the mutation rate is low or vanishes, nonlinear mobility alters the stability of the spiral waves and is responsible for their far-field breakup.
In this section we report that a similar intriguing phenomenon also occurs in the case where the mobility of the individuals is implemented by the simple nearest-neighbor pairexchange (3) which result in linear diffusive terms in the corresponding PDEs (4). To  [24][25][26]28] (where N = 1). Based on these previous works, see, e.g. [24-26, 28, 43, 45], we would anticipate that the dynamics of such a variant of the model would be characterized by the formation of stable spiral waves. As shown in Fig. 5 (top, left), this is indeed the case when the pair-exchange is sufficiently high (δ = 0.4). However, when δ is lowered the spiral waves become far-field unstable, see Fig. 5, after the shortening of their wavelength according to the scaling λ ∼ √ δ. Hence, as the wavelength is reduced under low mobility, it appears that the core of the spirals can sustain its arms only for a few rounds before a convective instability starts growing, very much like in the EI phase, and eventually cause the far-field breakup of the spiral waves. While the detailed mechanism of this far-field breakup has still to be elucidated, we believe that it stems from the nonlinear nature of the problem since demographic noise is not at its origin (for N = 1024 fluctuations are here negligible). We also think that a spiral far-field breakup always arises when ζ ≈ 0 and µ = 0, but depending on the mobility rate δ it occurs outside the lattice (high mobility rate) or within the grid (low mobility rate). In fact, a careful analysis of the PDEs (4) explains this phenomenon of far-field breakup with µ = 0: For fixed σ > 0, spiral waves exhibit far-field breakup when ζ σ/2 or when ζ is close to zero.
However, this analysis is beyond the scope of this paper and will be given elsewhere [63]. .

V. SUMMARY & CONCLUSION
We have studied the influence of a simple form of mobility, modeled by a pair-exchange between nearest neighbors, on the spatio-temporal patterns emerging in the generic twodimensional model of Refs. [43,45] for the cyclic rock-paper-scissors competition between three species. The underlying evolutionary processes are cyclic dominance-removal and cyclic dominance-replacement interactions, reproduction, migration (pair-exchange), and mutation. While various properties of this system, formulated as a metapopulation model, were investigated in Refs. [43,45] by combining multiscale and size expansions with numerical simulations, here we have analyzed the influence of the simple pair-exchange process on the properties of the spiraling patterns characterizing the dynamics of this system.
First, we have highlighted resolution issues that arise on finite lattices: While the description of the dynamics in terms of the underlying partial differential equations predict the formation of spiraling patterns in the so-called "bound-state phase", spiral waves may simply be elusive in the simulations on a finite lattice if the mobility rate is too low. More precisely, when the size of the lattice is finite and comparable to the size of experimental plates, we show that well-defined spiral waves can be observed only when the mobility rate is within a finite range: When the mobility rate is too low, the PDEs predict the emergence of spiral waves of wavelength of order of the lattice spacing which cannot be resolved, whereas when the mobility rate is sufficiently high the resulting spiraling patterns have a wavelength of the order of the lattice size and appear like planar waves. Spiral waves can be observed in lattice simulations when the mobility rate is between these two values. In fact, building on the analysis carried out in Refs. [43,45] in terms of the system's complex Ginzburg-Landau equation, we have estimated the critical value of the mobility rate. While the range within which spiral waves in bound-state phase grows with the size of the system, we have found that such a range may be finite and therefore spiraling patterns too small to be resolved and observable on lattices of a size comparable to the plates used in most experiments (96-well plates, see e.g. [8]). We believe that these "resolution issues" may therefore be particularly relevant when one tries to interpret experimental results and can explain why spiral waves appear to be elusive in microbial experiments as those of [6,8].
Second, we have focused on the version of the model with cyclic dominance-removal and without mutations that has received significant attention in the recent years, see e.g. [24][25][26]. While previous works reported that in this case the underlying rock-paper-scissors dynamics (with ζ = µ = 0) leads to well-defined spiraling patterns, here we show that spiraling patterns become convectively unstable and that a far-field breakup occurs when the mobility rate is lowered (with the other rates kept constant). We have verified that the mechanism underlying this phenomenon does not originate from demographic fluctuations and refer to the future work for a detailed analysis of its mechanism in terms of the system's partial differential equations [63].
While we have specifically discussed the biologically-relevant case of the two-dimensional metapopulation model, this analysis can be readily extended to one-dimensional and threedimensional lattices on which we would respectively expect traveling and scroll waves instead of spiral waves. It would also be interesting to study whether similar phenomena would arise to the oscillating patterns characterizing some RPS games on small-world networks [64][65][66].

VI. ACKNOWELDGMENTS
The support to BS via an EPSRC PhD studentship (Grant No. EP/P505593/1) is gratefully acknowledged.
The authors declare no conflict of interest.