In this work, we performed stochastic simulations of a cyclic nonhierarchical system composed of 5 species. To this purpose, we implemented a standard numerical algorithm largely used to study spatial biological systems11,13,41. We considered a generalisation of the rock-paper-scissors game for 5 species, whose rules are illustrated in Fig. 1a. The arrows indicate a cyclic dominance among the species. Accordingly, individuals of species i beat individuals of species (i+1), with (i=1,2,3,4,5).
The dynamics of individuals’ spatial organisation occurs in a square lattice with periodic boundary conditions, following the rules: selection, reproduction, and mobility. We assumed the May-Leonard implementation so that the total number of individuals is not conserved43. Each grid point contains at most one individual, which means that the maximum number of individuals is ({mathcal {N}}), the total number of grid points.
Initially, the number of individuals is the same for all species, i.e., (I_i,=,{mathcal {N}}/5), with (i=1,2,3,4,5) (there are no empty spaces in the initial state). We prepared the initial conditions by distributing each individual at a random grid point. At each timestep, one interaction occurs, changing the spatial configuration of individuals. The possible interactions are:
Selection: (i j rightarrow i otimes ,), with (j = i+1), where (otimes) means an empty space; every time one selection interaction occurs, the grid point occupied by the individual of species (i+1) vanishes.
Reproduction: (i otimes rightarrow i i,); when one reproduction is realised an individual of species i fills the empty space.
Mobility: (i odot rightarrow odot i,), where (odot) means either an individual of any species or an empty site; an individual of species i switches positions with another individual of any species or with an empty space.
In our stochastic simulations, selection, reproduction, and mobilities interactions occur with the following probabilities: s, r and m, respectively. We assumed that individuals of all species have the same probabilities of selecting, reproducing and moving. The interactions were implemented by assuming the von Neumann neighbourhood, i.e., individuals may interact with one of their four nearest neighbours. The simulation algorithm follows three steps: i) sorting an active individual; ii) raffling one interaction to be executed; iii) drawing one of the four nearest neighbours to suffer the sorted interaction (the only exception is the directional mobility, where the neighbour is chosen according to the movement tactic). If the interaction is executed, one timestep is counted. Otherwise, the three steps are redone. Our time unit is called generation, defined as the necessary time to ({mathcal {N}}) timesteps to occur.
In our model, individuals of one out of the species can move into the direction with more individuals of a target species. The choice is based on the strategy assumed by species. We assumed three sorts of directional movement tactics:
Attack tactic: an individual of species i moves into the direction with more individuals of species (i+1);
Anticipation tactic: an individual of species i goes towards the direction with a larger number of individuals of species (i+2);
Safeguard tactic: an individual of species i walk into the direction with a larger concentration of individuals of species (i-2).
In the standard model, individuals of all species move randomly.
We considered that only individuals of species 1 perform the directional movement tactics, as illustrated in Fig. 1b. The solid, dashed, and dashed-dotted lines represent the Attack, Anticipation, and Safeguard tactics, respectively. The concentric circumference arcs show that individuals of species 2, 3, 4, and 5 always move randomly. For implementing a directional movement, the algorithm follows the steps: i) it is assumed a disc of radius R (the perception radius), in the active individual’s neighbourhood; ii) it is defined four circular sectors in the directions of the four nearest neighbours; iii) according to the movement tactic, the target species is defined: species 2, 3, and 4, for Attack, Anticipation, and Safeguard tactics, respectively; iv) it is counted the number of individuals of the target species within each circular sector. Individuals on the borders are assumed to be part of both circular sectors; v) the circular sector that contains more individuals of the target species is chosen. In the event of a tie, a draw between the tied directions is made; vi) the active individual switches positions with the immediate neighbour in the chosen direction. The swap is also executed in case of the neighbour grid point is empty.
To observe the spatial patterns, we first performed a single simulation for the standard model, Attack, Anticipation, and Safeguard tactics. The realisations run in square lattices with (500^2) grid points, for a timespan of 5000 generations. We captured 500 snapshots of the lattice (in intervals of 10 generations), that were used to make the videos of the dynamics of the spatial patterns showed in https://youtu.be/Ndvk6Rg57m4 (standard), https://youtu.be/JGhkDAHSo74 (Attack), https://youtu.be/ZZp9QlOfv2Q (Anticipation), and https://youtu.be/eFxWdLhIOuQ (Safeguard). The final snapshots were depicted in Fig. 2a–d. Individuals of species 1, 2, 3, 4, and 5 are identified with the colours ruby, blue, pink, green, and yellow, respectively; while white dots represent empty spaces. The simulations were performed assuming selection, reproduction, and mobility probabilities: (s = r = m = 1/3). The perception radius was assumed to be (R=3).
The population dynamics were studied by means of the spatial density (rho _i), defined as the fraction of the grid occupied by individuals of species i at time t, i.e., (rho _i = I_i/{mathcal {N}}), where (i=0) stands for empty spaces and (i=1,…,5) represent the species 1, 2, 3, 4, and 5. The temporal changes in spatial densities of the simulations showed in Fig. 2 were depicted in Fig. 3, where the grey, ruby, blue, pink, green, and yellow lines represent the densities of empty spaces and species 1, 2, 3, 4, and 5, respectively. We also computed how the selection risk of individuals of species i changes in time. To this purpose, the algorithm counts the total number of individuals of species i at the beginning of each generation. It is then counted the number of times that individuals of species i are killed during the generation. The ratio between the number of selected individuals and the initial amount is defined as the selection risk of species i, (zeta _i). The results were averaged for every 50 generations. Figure 4 shows (zeta _i,(%)) as a function of the time for the simulations presented in Fig. 2. The ruby, blue, pink, green, and yellow lines show the selection risks of individuals of species 1, 2, 3, 4, and 5, respectively.
To quantify the spatial organisation of the species, we studied the spatial autocorrelation function. This quantity measures how individuals of a same species are spatially correlated, indicating spatial domain sizes. Following the procedure carried out in literature41,42,44,45,46, we first calculated the Fourier transform of the spectral density as (C({{vec{r^{prime}}}}) = {mathcal{F}}^{{ – 1}} { S({{vec{k}}})} /C(0)), where the spectral density (S({{vec{k}}})) is given by (S({{vec{k}}}) = sumlimits_{{k_{x} ,k_{y} }} {mkern 1mu} varphi ({{vec{kappa }}})), with (varphi ({{vec{kappa }}}) = {mathcal{F}}{mkern 1mu} { phi ({{vec{r}}}) – langle phi rangle }). The function (phi ({{vec{r}}})) represents the species in the position ({{vec{r}}}) in the lattice (we assumed 0, 1, 2, 3, 4, and 5, for empty sites, and individuals of species 1, 2, 3, 4, and 5, respectively). We then computed the spatial autocorrelation function as
$$C(vec{r^{prime}}) = sumlimits_{{|{{vec{r^{prime}}}}| = x + y}} {frac{{C({{vec{r^{prime}}}})}}{{min (2N – (x + y + 1),(x + y + 1))}}}.$$
Subsequently, we found the scale of the spatial domains of species i, defined for (C(l_i)=0.15), where (l_i) is the characteristic length for species i.
We calculated the autocorrelation function by running 100 simulations using lattices with (500^2) grid points, assuming (s = r = m = 1/3) and (R=3). Each simulation started from different random initial conditions. We then captured each species spatial configuration after 5000 generations to calculate the autocorrelation functions. Finally, we averaged the autocorrelation function in terms of the radial coordinate r and calculated the characteristic length for each species. We also calculated the standard deviation for the autocorrelation functions and the characteristic lengths. Figure 4 shows the comparison of the results for Attack, Anticipation, and Safeguard strategies with the standard model. The ruby, blue, pink, green, and yellow circles indicate the mean values for species 1, 2, 3, 4, and 5, respectively. In the case of standard model, the mean values are represented by grey circles, which are the same for all species. The error bars show that standard deviation. The horizontal black line represents (C(l_i), =, 0.15).
To further explore the numerical results, we studied how the perception radius R influences species spatial densities and selection risks. We calculated the mean value of the spatial species densities, (langle , rho _i,rangle) and the mean value of selection risks, (langle , zeta _i,rangle) from a set of 100 simulations in lattices with (500^2) grid points, starting from different initial conditions for (R=1,2,3,4,5). We used (s,=r,=,m,=1/3) and a timespan of (t=5000) generations. The mean values and standard deviation were calculated using the second half of the simulations, thus eliminating the density fluctuations inherent in the pattern formation process. The results were shown in Fig. 6, where the circles represent the mean values and error bars indicate the standard deviation. The colours are the same as in Figs. 3 and 4. Furthermore, to verify the precision of the statistical results, we calculated the variation coefficient – the ratio between the standard deviation and the mean value. Supplementary Tables S1 and S2 show statistical outcomes for species densities and selections risks, respectively.
We studied a more realistic scenario where not all individuals of species 1 can perform the directional movement tactics. For this reason, we defined the conditioning factor (alpha), with (0,le ,alpha ,le ,1), representing the proportion of individuals of species 1 that moves directionally. For (alpha =0) all individuals move randomly while for (alpha =1) all individuals move directionally. This means that every time an individual of species 1 is sorted to move, there is a probability (alpha) of the algorithm implementing the directional movement tactic, instead of randomly choosing one of its four immediate neighbours to switch positions. To understand the effects of the conditioning factor, we observed how the density of species 1 changes for the entire range of (alpha), with intervals of (Delta alpha = 0.1). The simulations were implemented for (R=3) and (s,=r,=,m,=1/3). It was computed the mean value of the spatial density of species 1, (langle , rho _1,rangle), and its standard deviation from a set of 100 different random initial conditions. The results were depicted in Fig. 7, where the green, red, and blue dashed lines show (langle , rho _1,rangle) as a function of (alpha). The error bars indicate the standard deviation.
Finally, we aimed to investigate how the directional movement tactics jeopardise species coexistence for a wide mobility probability range. Because of this, we run 2000 simulations in lattices with (100^2) grid points for (0.05,<,m,<,0.95) in intervals of (Delta , m, =,0.05), with (R=2) and (R=4). The simulations started from different random initial conditions and run for a timespan of 10000 generations. Coexistence happens if at least one individual of all species is present at the end of the simulation, (I_i (t=5000) ne 0) with (i=1,2,3,4,5). Otherwise, the simulation results in extinction. The coexistence probability is the fraction of implementations which results in coexistence. The simulations were performed for two values of perception radius, (R=2) and (R=4); the selection and reproduction probabilities were assumed to be (s,=,r,=,(1-m)/2). The results were depicted in Fig. 8, where yellow, green, red, and blue lines show the coexistence probability as a function of m for the standard model, Attack, Anticipation, and Safeguard tactics, respectively. The solid and dashed lines show the results for (R=2) and (R=4), respectively.
Source: Ecology - nature.com
