Environment driven oscillation in an off-lattice May–Leonard model

Cyclic dominance of competing species is an intensively used working hypothesis to explain biodiversity in certain living systems, where the evolutionary selection principle would dictate a single victor otherwise. Technically the May–Leonard models offer a mathematical framework to describe the mentioned non-transitive interaction of competing species when individual movement is also considered in a spatial system. Emerging rotating spirals composed by the competing species are frequently observed character of the resulting patterns. But how do these spiraling patterns change when we vary the external environment which affects the general vitality of individuals? Motivated by this question we suggest an off-lattice version of the tradition May–Leonard model which allows us to change the actual state of the environment gradually. This can be done by introducing a local carrying capacity parameter which value can be varied gently in an off-lattice environment. Our results support a previous analysis obtained in a more intricate metapopulation model and we show that the well-known rotating spirals become evident in a benign environment when the general density of the population is high. The accompanying time-dependent oscillation of competing species can also be detected where the amplitude and the frequency show a scaling law of the parameter that characterizes the state of the environment. These observations highlight that the assumed non-transitive interaction alone is insufficient condition to maintain biodiversity safely, but the actual state of the environment, which characterizes the general living conditions, also plays a decisive role on the evolution of related systems.

According to the Darwinian selection hypothesis only the most viable competitor should survive as a result of a selection process. But we witness an amazing diversity of species in nature, which begs alternative explanations in ecology and in other complex competitive systems. The presence of a cyclic dominance among competitors is an elegant and very simple clue to resolve this contradiction. Indeed, scientists have observed several cases in living systems where the mentioned type of interaction can be observed. Examples can be given from microbial and plant communities, to coral reef, lizards, salmons and human interactions [1][2][3][4][5][6] . We should stress, however, that similar cyclic dominance can also be detected in so-called social systems, where different strategies may dominate each other in a non-transitive way 7,8 .
The basic model describing such kind of interactions between system members is based on the well-known rock-scissors-paper game where every member is a predator of another member and a prey for the third one simultaneously. Naturally, the strength of the dominance could be different between some predator-prey pairs and this asymmetry can provide some counter-intuitive system behaviors. One of them is the so-called "survival of the weakest" effect where the species having the lowest invasion rate develops the highest fraction in the population 9,10 . This phenomenon was confirmed in several modified models during the years 8,[11][12][13] and in general the related dynamical behavior of cyclically dominant systems has collected significant research interest in the last decade [14][15][16][17][18][19][20][21] . For example, the n number of cyclically competing species can be extended to an arbitrary number where n plays a central role on the resulting dynamics. This generalized version has a very rich structure leading to the formation of multi-domains of one or more species, which are separated by interfaces 22,23 . Also, the increase of the number of species usually leads to the development of more complex dynamical patterns. By focusing on the interplay between competition and partnership in spatial environments it can be observed that the development of neutral associations between individuals belonging to enemy partnership seems to affect www.nature.com/scientificreports/ the development of the dynamical structure along interfaces separating competing domains. For further details and a general overview of the present state of this research avenue we refer the interested reader to recent review papers [24][25][26] .
Technically the related problems can be studied in Lotka-Volterra and May-Leonard models where the spatial distribution of species have a decisive role on the final evolutionary outcomes [27][28][29][30] . In Lotka-Volterra models the application of 3 ≤ n species offers the simplest extension of a rock-scissors-papers-like cyclic dominance where predation and reproduction may occur in an elementary process [31][32][33] . In a May-Leonard model the cyclic invasion is split into a "selection" and a probabilistic reproduction step which makes the sum of all individuals a nonconserved quantity. Strongly related to the scope of our present study, it turned out that mobility has a decisive role on the evolving pattern [34][35][36] . More precisely, when the typical length of rotating spirals become comparable to the system size due to strong diffusion then the system can easily evolves into a trapped, or absorbing state where only a single species survives.
The above mentioned unequal invasion rates could be the result of an environmental factor, which is in general a parallel research avenue in complex systems. Heterogeneous environment can modify dynamical process directly [37][38][39][40][41][42][43] , which could be a local or a seasonal, or time-dependent change 44,45 . But we may control the state of the environment to modify the fractions of competing agents intentionally [46][47][48] . Furthermore, the actual state of the environment can determine the vitality of the whole population fundamentally because adverse conditions may prevent individuals to survive while beneficial environment with unlimited resources can offer optimal living conditions, hence supports species. A natural question is how general environmental conditions influence the established cyclic dominance in the whole population. Is there any consequence on the evolving patterns when the environment makes easy or difficult for species to reproduce? An extreme case could be when a possible death of the individuals due to starvation is considered, which happens when a certain individual fails a given number of times when attempting predation 49 . In this case it was observed that the death of these individuals provide a crucial contribution to preservation of coexistence.
As a first step toward a more comprehensive understanding, in the following we study a model where the general state of the environment is modeled via a single parameter which determines the local carrying capacity of the system. In this way we can vary the living conditions of all competitors uniformly and monitor how such changes influence the resulting evolutionary outcome. We note, however, that similar question was raised by other scientists previously who studied a well-mixed system or a spatially structured metapopulation 50,51 . For our present study it is important to stress that emerging rotating spirals and spiral waves are frequently observed accompanying patterns of cyclic dominance in spatial systems [51][52][53][54][55] . Therefore it is a fundamental question to study these arrangement when external conditions are varied. From this viewpoint there is a crucial technical circumstance that need to be mentioned. Generally, the application of lattice-type interaction topology makes simulations significantly easier, while the most important behaviors are still observable in these systems [56][57][58][59] . However, there is a drawback of the mentioned modeling technique in our present case which has a paramount importance. In particular, a lattice topology allows to change external conditions by discrete steps only. An alternative technique could be the so-called off-lattice simulations where the positions of individuals, hence their neighborhood may change continuously [60][61][62][63] . The latter makes us possible to tune system parameters almost continuously, hence the control parameter which characterizes of the status of the environment can be varied finely. Indeed, the latter technique is more demanding and requires larger numerical efforts, but in certain cases we cannot avoid this difficulty. For example this is the proper way to study certain phenomena, like clustering 61,64 .
As we will show, the emerging spatio-temporal pattern depends sensitively on the actual state of the external environment which directly determines the general living condition of the whole population. In what follows, we first present the suggested off-lattice version of the May-Leonard model and its mathematical details. We then proceed with the presentation of the main results and followed by a discussion of their wider implications.

Model specification
In this paper we shall consider a square box of linear size L = 1 with periodic boundary conditions, which is the stage of our off-lattice simulations within the framework of a May-Leonard model. However, we should note that because of off-lattice character of the simulation the actual shape has no particular significance on the final outcomes. According to the traditional setup three different species, A, B and C, are fighting for space where they dominate each others cyclically. In particular, the species A preys the species B, that preys the species C and the C preys the species A, hence closing the cycle. The predation only occurs if there is a prey inside a circle of radius ℓ p (predation length) centered around the predator and in this case the closest prey dies out. If there is no prey within the circle, then nothing happens.
An alternative elementary process is a prey-independent reproduction of the focal individual. In this case a successor emerges within a radius ℓ m , but only if the total number of all individuals within the reproduction range ( ℓ r ) is smaller than M. The latter parameter, called local carrying capacity, characterizes the actual state of the environment. In a harsh environment the value of M is low because limited resources cannot keep more individuals alive, while in a benign environment this value is higher. The third microscopic process is individual movement. In this case the chosen individual moves by a distance ℓ m (movement length) in a randomly chosen direction. Notably, this step is always executed, while the success of predation and reproduction may depend on other circumstances, like the presence of a prey, or the total number of individuals in the reproduction neighborhood.
Initially n A = n B = n C = 10 4 individuals of species are distributed randomly where the horizontal and vertical coordinates of every individuals are continuous variables. The total number of the whole population is N = n A + n B + n C , which may change in time due to the above mentioned stochastic processes. It is worth www.nature.com/scientificreports/ noting that when we calculate the proper distance of two individuals we consider the mentioned periodic boundary conditions. During an elementary Monte Carlo (MC) step an active individual is chosen randomly, which may move, predate or reproduce. The related probabilities are m, p and r, respectively. For a full MC step we repeat the elementary steps N times. In the following we have chosen m = 0.5 and p = r = 0.25 unless otherwise stated. Furthermore, for the characteristic lengths we used ℓ p = ℓ r = 0.02 , and ℓ m = 0.01 . These values allow us to observe proper behavior of the spatial system. But we stress that similar observations can be made if we use other values of our model parameters therefore in the following we present characteristic and typical system reactions in dependence of environmental changes. We should also emphasize that by choosing too large length values, when ℓ ≈ 1 , the scales of microscopic steps become comparable to the system size. In this extreme case we would terminate onto a well-mixed system where the actual spatial distribution of individuals has no particular importance. To obtain the expected accuracy for the above mentioned parameter values we have repeated every run 10 3 times by using independent initial conditions and averaged the individual results. Further details of our numerical experiments are given in the next section.

Results
We first present a general overview about the impact of environmental change on the emerging spatial pattern of competing species. Figure 1 depicts six representative snapshots of our off-lattice May-Leonard-model for different values of local carrying capacity M. The competing A, B, and C species are marked by red, blue, and yellow colors respectively, while white marks empty space. The comparison illustrates very clearly that the rotating spirals become more evident as we increase the value of M. In parallel, the total number of individuals also increases by increasing M because the portion of white color becomes gradually smaller.
Our last observation is summarized in Fig. 2 in a more quantitative way. Here we plot the N average number of individuals in the whole population in dependence of M. We stress that this is an average value because the temporary number of individuals may change in time. To illustrate it, in the inset we show a representative distribution of N in the stationary state for M = 15 . Here the position of the calculated average is marked by a vertical red line. Note that the error bars are also marked in the main plot and to obtain the requested accuracy we calculated the time average over 4 × 10 4 MC steps after 10 8 MC steps of relaxation for each M value. We note www.nature.com/scientificreports/ that the actual number of iterations depends on the size of population, therefore to obtain the related data for higher M requires significantly larger numerical efforts. As the main plot suggests, the average value follows a linear dependence on M where the value of the slope is quite robust and does not depend on microscopic details, like the ℓ m value. In particular, if we allow more intensive individual movement, for instance, then the increment of total individuals will change in a similar way, having the same slope, as we increase M.
As we already noted, we have also explored what happens in the extreme cases when a characteristic length of movement becomes comparable to the system size. For example when ℓ m is increased too long, say ℓ m = 0.8 then the typical domain sizes grow which shows similar behavior can be observed in a well-mixed system. In this case the fluctuations can be so large that one of the competing species goes extinct, which breaks the symmetry and the population eventually terminates to a homogeneous state. This behavior is similar to those previously observed on a lattice structure 34 . Another interesting behavior can be observed when the predation and reproduction lengths become comparable to the system size. In this case empty space behaves like an "additional" species and the real species become seemingly isolated from each other. As a result, the deserted areas occupy a significant portion of available space. We will discuss this phenomenon later but first we focus on other aspects of emerging spirals.
The emergence of rotating spirals has a detectable consequence not only on spatial, but also on temporal patterns. The latter fact can be illustrated properly if we monitor the time dependence of the fraction of a certain species. This phenomenon is shown in Fig. 3 where we plot the temporary number of individuals for species A for different values of M. As the plot suggests the oscillation becomes more intensive for larger M values. Obviously, similar patterns can be obtained for the remaining two species because the non-transitive interaction establishes a symmetry among the competing species. While for small M values, when the environment is harsh, the time course seems to be noisy, but for high values of M the environment becomes rich of resources hence it is capable to maintain a large populations stably. This can also be read out from the plot because the average level of n A increases gradually as we increase the value of M. There is, however, an important feature of the time dependence which underlines the main conclusion of our study. In a stochastic simulation it is a generally expected behavior that for larger population the system behavior becomes less noisy. Indeed, this also happens in our model and the oscillation becomes more regular as we reach higher M values, hence indirectly enlarge the size of the whole  www.nature.com/scientificreports/ population. The amplitude of the oscillation, however remains significant hence indicating the emergence of spirals and waves we already have shown in Fig. 1.
To give a more quantitative description about the oscillations we apply the Fourier analysis to the ρ(t) function which determines the fraction of a species in time. The temporal discrete Fourier transform can be given as where the coefficient of f is calculated from N G = 10 4 components. To collect reliable data for the stationary states we always used 1000 MC steps of relaxation. The resulting power spectrum �|ρ A (f )| 2 � of species A is shown in Fig. 4 where we plotted the curves obtained for different M values simultaneously. These values are indicated in the legend. We note that for an appropriate scale all power spectrum values are multiplied by a 10 5 constant factor. For the requested accuracy we averaged the data over 250 independent simulations where the system evolution was launched from different initial states. Evidently, similar curves can be obtained for the remaining B and C species. The comparison of curves indicates that the location of the peak shifts toward higher frequency values as we increase M. For example for M = 30 we have 107 oscillations during 10 4 Monte Carlo steps. In parallel, the height of the peak also increases by enlarging M, which suggests that the oscillation becomes more characteristic as the living conditions are improved. Notably, the position of the peaks show a nice logarithmic scaling as it is shown in the inset of Fig. 4. This quantitative analysis confirms what we already observed in Fig. 3. Namely, even if there is a clear non-transitive cyclic interaction among competing species, the well-known rotating spirals and accompanying oscillations are hardly detectable if the environment is poor of resources and can only maintain a stunted population. We stress that the biodiversity is still maintained, but not in the presence of rotating spirals we frequently expect from a spatial system having cyclic dominant microscopic dynamic. However, if the living conditions are improved then the anticipated rotating spirals of spatial distribution and the time-dependent oscillation of species become evident again.
Before concluding we should highlight that alone the introduced carrying capacity parameter is just an initial step to model the general conditions of the environment. More precisely, the single value of M does not determine the state of the environment accurately, because this parameter is linked to the reproduction length ℓ r . When the latter is large then even a relatively high M value still represents a poor environment. This is demonstrated in Fig. 5 where we applied five times larger reproduction and predation lengths as previously. Therefore, even if we used quite large M values, the total sum of individuals remain low in the population. As an accompanying effect, the rotating spirals diminish from the pattern. But they can be recovered again if we increase M, as it is done in panels (b) and (c) of Fig. 5. Interestingly, the portion of empty space remains high, and seems to behave as an organic additional member of the spirals. This is a phenomenon that cannot be observed in a cyclic system where the size of the population is strictly constant in time 51 .

Discussion
It is a well-known fact that cyclic dominance among competing species not only maintains diversity but frequently generates spiral waves in a spatial system where agents have limited access to other interacting partners. There are, however, some circumstances when this general picture is broken. An example could be when the high mobility of species destroys the above described patterns and jeopardizes the coexistence of all competitors 34 . Similarly, the breaking of unidirectional invasions, or heterogeneity in site-specific invasion rates could also www.nature.com/scientificreports/ terminate the coexistence of competing species 12,65 . External factors, like the proper state of environment, which generally determines the viability of the population, also seems to be a crucial ingredient to this problem. Motivated by this argument, in this work we explored how the actual state of environment influences directly the competition of equally strong opponents. In our simple model the mentioned condition can be controlled via a single parameter that determines the local carrying capacity of the system. The other central feature of our model was the off-lattice topology which made possible to vary this parameter and the related external condition gently. Our first observation is the stable coexistence of competing species for low value of M when the general maintaining capacity of the environment is moderate. But in this limit, spiral patterns, which characterizes the traditional three-species model, cannot be detected. In parallel, the time dependence of a certain species shows an irregular sequence. This anomaly, however, can be restored if we enlarge M, which involves a more supporting environment hence can maintain a more crowded population. It is worth noting that conceptually similar observation was analyzed by Szczesny, Mobilia, and Rucklidge, who considered a metapopulation model where every patches, which are organized in a square grid, have a limited carrying capacity 51 . In their cases when the sizes of the well-mixed sub-populations were small the rotating spirals disappeared, but they were restored again by increasing the size of the mentioned sub-populations. The agreement between their partial differential equation approach and our present off-lattice simulations underlines the universality and broadens the validity of the presented observation.
We also demonstrated that there is a linear dependence of the average size of the population on the M parameter and this slope is largely independent on microscopic details of the used May-Leonard model. We have also studied how a regular oscillation of time dependent fraction of species emerges as we improve the general quality of the environment. Furthermore, we demonstrated that the characteristic frequency of the oscillation shows a logarithmic scaling with the introduced M parameter.
Summing up, our work pointed out that the quality of the environment, which provides a stage for competition of species and determines the general living conditions, can be a decisive factor for what kind of spatial and temporal patterns emerge. When the general living conditions are poor, because of the lack of resources or for other reasons, then the well-known rotating spirals characterizing such kind of spatial system are missing, but they become detectable again if the living conditions are improved. Similar behavior was already observed in a case when intensive diffusion resulted in the disappearance of spirals. This last effect, however, can be easily explained because fast individual movement helps more intensive mixing hence it drifts the spatial system toward the well-mixed behavior. But in our case the poor environment has not similar consequence therefore the lack of rotating spirals seems to be surprising at first sight. On the other hand the reported behavior is robust because other modeling approach, by assuming metapopulation setup, also confirmed it. Nevertheless, this phenomenon can also explain why we can observe spiral patterns rarely in field studies where cyclic dominance is present otherwise. . This figure illustrates that alone the value of M is insufficient to characterize the state of environment, because its combination with the length scales and the probabilities determining dynamical processes could tell us the proper living conditions. The increase of M, however, when all other parameters are fixed, can restore the rotating spirals we expect from a cyclically dominant spatial system. Interestingly, the large predation and reproduction ranges result in large portion of empty space which behaves as an additional inseparable part of the spirals.