Universal scaling of extinction time in stochastic evolutionary dynamics

Evolutionary dynamics is well captured by the replicator equations when the population is infinite and well-mixed. However, the extinction dynamics is modified with finite and structured populations. Experiments on the non-transitive ecosystem containing three populations of bacteria found that the ecological stability sensitively depends on the spatial structure of the populations. Based on the Reference–Gamble–Birth algorithm, we use agent-based Monte Carlo simulations to investigate the extinction dynamics in the rock-paper-scissors ecosystem with finite and structured populations. On the fully-connected network, the extinction time in stable and unstable regimes falls into two universal functions when plotted with the rescaled variables. On the two dimensional grid, the spatial structure changes the transition boundary between stable and unstable regimes but doesn’t change its extinction trend. The finding of universal scaling in extinction dynamics is unexpected, and may provide a powerful method to classify different evolutionary dynamics into universal classes.

www.nature.com/scientificreports/ community in the high-Arctic tundra in Greenland 61 and the plant communities 62,63 . These findings in nature make rock-paper-scissors game not only a theoretically interesting model but also a practical system to study.
In this Report, we use agent-based Monte Carlo simulations to investigate the extinction dynamics in the rock-paper-scissors ecosystem with finite and structured populations. Making use of numerical simulations, one can extract the extinction time in different regimes. However, for the finite and structured populations, the usual Moran algorithm 64,65 and its generalized versions 24,66,67 may not be appropriate. One shall turn to the evolutionary theory on graphs 21,[68][69][70] to include the local mutual competitions. Following the same spirit, we introduce a local three-party Reference-Gamble-Birth (RGB) algorithm and perform the agent-based Monte Carlo simulations for the evolutionary dynamics. The RGB algorithm only requires local information in the stochastic processes and captures the network structure into account. And, it correctly reproduces the replicator dynamics when the populations are infinite and well-mixed.
Two major results are found in our numerical simulations. The first one is the universal scaling: the extinction time is always finite and falls onto the universal function when plotted with the rescaled variables. On the fullyconnected network, two universal functions for stable and unstable regimes are found and the phase transition between them corresponds to neutral (critical) dynamics where the extinction time shows power-law behavior. The unexpected universal scaling may provide a powerful method to classify ecosystems with different dynamical behaviors into universal classes. The second major result is the notion of global payoff versus local payoff. On the two dimensional grid, the spatial structure changes the transition boundary between stable and unstable regimes but the universal function remains the same. This observation motivates us to introduce the notion of global payoff which turns out to be consistent with our numerical simulations. Our numerical simulations not only resolve the conflicting results in the previous studies 33,34 but also suggest a potential tool to classify different evolutionary dynamics by universal scaling.

Methods
Rock-paper-scissors game. Here we consider the symmetric cyclic competitions between species A, B and C described by the payoff matrix where the positive-definite parameter g denotes the gain in the cyclic competitions. The replicative dynamics of the ecosystem with populations N A , N B , N C for the species A, B, C is captured by the following differential equations, where a(t) ≡ N A /N is the frequency (relative population) of species A in the total population N = N A + N B + N C . The other frequencies b(t) ≡ N B /N and c(t) ≡ N C /N are defined in similar fashion. The parameter φ = (g − 1)(ab + bc + ca) is the average fitness of the ecosystem.
The symmetric replicator equations 11,34 host a non-trivial fixed point (a * , b * , c * ) = ( 1 3 , 1 3 , 1 3 ) . If the gain is greater than unity, g > 1 , the flows around the fixed point is stable. On the other hand, if the gain is less than unity, g < 1 , the fixed point becomes unstable. At the critical gain g c = 1 , it corresponds to the zero-sum game and the flow around the fixed point is neutral. Therefore, according to the replicator equations, the ecosystem goes through a phase transition from the stable regime to the unstable one when the gain g is tuned through the critical point g c = 1. RGB algorithm. The above scenario is well known for the rock-paper-scissors game when the population is infinite and well-mixed. But, what happens when the population is finite and the competitions among different species occur on the network with specific spatial structure? To investigate these important issues, we utilized the agent-based Monte Carlo simulations. To incorporate the locality of mutual competitions on the reasonable ground, we introduce a local three-party Reference-Gamble-Birth (RGB) algorithm. For each RGB update, it consists four steps as shown in Fig. 1. In the RGB algorithm, the third-party Reference provides local survival standard and determines the death probability of Gamble.
The death matrix D plays a central role in RGB algorithm. Similar to payoff matrix, its matrix elements D ij defining the death probability for all possible competitions among species. The definition of death matrix indicates that the death matrix share the same matrix dimension with the payoff matrix and all its matrix elements must sit inside the range 0 ≤ D ij ≤ 1 . In this section, we are going to show (i) the equivalent dynamics between replicator equations and the deterministic dynamics of stochastic updates in RGB algorithm in well-mixed and infinite populations; (ii) a simple scaling relation, i.e. D = P/s , to construct the death matrix D for RGB algorithm from an arbitrary negative payoff matrix P; (iii) extracting the extinction time (i.e. a simple scaling T ex = T RGB /s ) for a payoff matrix P by using the RGB algorithm.
In a single update, there are two possible ways to change the population size N i of species i: (1) N i = 1 with probability D kj x i x j x k coming from the process of the decease of a species-k Gamble when contacting with a The RGB algorithm introduces selection among species by a death matrix, and therefore its deterministic dynamics is captured by the replicator equations without the birth terms. Now, we can derive the relation between the death matrix in RGB algorithm and the payoff matrix in replicator equations when frequency-dependent selection is taken into account. The gauge redundancy of the payoff matrix indicates the replicator equations remain invariant when adding a arbitrary constant to all elements in the payoff matrix 11 . Based on the gauge redundancy, we consider a negative payoff matrix P with negative values for all matrix elements (i.e. P ij ≤ 0 for all i and j) without loss of generality. Because 0 ≤ D ij ≤ 1 , a proper rescaling factor is required to establish the equivalence between the payoff matrix P ij and the death matrix D ij . Introducing a proper rescaling factor s, which has a lower bound s m = max(|P ij |) defined by the payoff matrix, we define the death matrix from the negative payoff matrix as The RGB algorithm has several advantages over the Moran process. First of all, only local information is involved in the stochastic process, more realistic than the need of global information for the Moran process. Secondly, the RGB algorithm reproduces the dynamics of frequency-dependent selection described by the replicator equations in infinite and well-mixed populations when the rescaling of time is restored properly. Finally, the the RGB algorithm takes the network structure into account, capturing the sensitive dependence of evolutionary dynamics in structured populations.

Results
Universal scaling in well-mixed populations. We perform agent-based Monte Carlo simulations to extract the mean extinction time. For each run, when one of the species goes extinct, we stop the simulated evolution and mark the time of extinction. The mean extinction time T ex (N, g) , depending on the population size N and the gain g defined in the payoff matrix, is computed by taking ensemble average of 10,000 runs in numerical simulations. We first focus on the well-mixed population, i.e. the fully connected network where every individual is connected to each other. It is rather remarkable that all numerical data can be collapsed onto the universal functions, where g ≡ |g − g c | denotes the deviation from the critical gain g c = 1 . As shown in Fig. 2, the numerical data for g > g c (stable regime) and g < g c (unstable regime) collapse onto the universal scaling functions F ± (x) respectively. The exponent α = 1 can be extracted from extinction time of different population sizes at the critical point g = 0 and the other exponent β = 1 is the optimal fitting to collapse all data with different N and g. In the following paragraphs, we would like to show that these scaling variables can be understood within the framework of the renormalization group (RG) 71 . But, it is worth mentioning that the scaling variable x = N(g − 1) agrees with the our expectation from population genetics 72 that the extinction time in finite populations will scale as the product of population size and selection strength.
Let us derive the scaling form of the extinction time 71 by RG analysis. Assuming the scaling dimensions for T ex , N and g are d T , d N and d g respectively,  Fig. 2). Our numerical simulations show that the extinction time is always finite but the ecological stability can still be defined through the trend of extinction time. Plotted with the rescaled variables, stable and unstable regimes fall onto two universal functions F ± (x) and can be differentiated without ambiguity. In general, the exact forms of the scaling functions F ± (x) are difficult to obtain. However, its asymptotic form can be understood by the following arguments. Consider the stable regime first. In the asymptotic limit ( |x| = N| g| ≫ 1 ), the extinction of the ecosystem arises from a chain of successive rare event against the flow back to the stable equilibrium. In the stable regime, the clusters within the correlation length are locally stable with a small probability p to go extinct. Suppose the typical size of the cluster is N c so that there are roughly N/N c clusters. The probability for the whole ecosystem to go extinct is thus p N/N c = q N , where q = p 1/N c . In consequence, the extinction time grows exponentially with the population size N and the asymptotic form of F + (x) should be Here A is some positive constant and N + = 1/(A�g) is the characteristic constant for the exponential growth. In the unstable regime, the population decays exponentially, N(t) = Ne −γ t . The extinction time is usually estimated by the criterion N(T ex ) ≈ O(1) , leading to the logarithmic trend T ex ∼ ln N . The asymptotic form of the scaling function F − (x) thus takes the form, The asymptotic forms in Eqs. 16, 17 has been derived in the previous studies 34 and agree with the numerical simulations presented in Fig. 2 rather well. Even though the asymptotic forms of the scaling functions can be understood rather well, the exact forms of the scaling functions F ± (x) are difficult to obtain. Global payoff matrix in structured populations. Now we turn to the extinction dynamics on the 2D grid. As explained in the previous paragraphs, only local neighbors participate in each RGB updates. It is striking that the biodiversity is enhanced due to the spatial structure of the network: the extinction time always follows the exponential trend as shown in Fig. 3. Even in the extreme case g = 0 , the stability of the cyclic-competing ecosystem is still protected by the grid structure. Is the extinction time described by the same scaling function as in the well-mixed case? The answer, after finding appropriate rescaled variable, turns out to be positive.
Even though the phase boundary between stable and unstable regimes are modified by the spatial structure of the network, the universal functions remain the same. We are thus inspired to introduce the notion of global payoff matrix as explained below. Local population dynamics generated by local interactions depends on the payoff matrix parametrized by the local gain g. But, the extinction time is a global property: even with the same local gain g, the 2D grid and the fully-connected network (well-mixed population) give different results. Because the fully-connected network looks the same either locally or globally, it is invariant under renormalization-group transformation and serves as a good reference. In RG analysis, the partition function is kept invariant to extract the renormalized coupling constants 71 . Here we adopt the same method and introduce the global gain G on an arbitrary network via the relation where the superscript FC stands for "fully-connected". That is to say, the extinction time for the local gain g on the considered network is set equal to that for the global gain G on the fully-connected network. According to the definition, g = G on the fully-connected network, reflecting the fact that local and global properties are the same in the well-mixed population. In our numerical simulations, the asymptotic limit is approximated by the largest population size in available simulations. To accurate the mapping, the characteristic constant in universal scaling function ( F + ) is used to defined the global payoff matrix on the 2D grid. The mapping between G and g on the 2D grid is carried out numerically (shown in Fig. 4). The red line ( G = g ) represents the trivial relation for (12) N =b d N N, www.nature.com/scientificreports/ the fully-connected network. All data sit in the stable regime (blue background) and confirm that the dynamics on the 2D grid is always ecologically stable with G > 1 . The mapping between the local and the global gains in the 2D grid reveals strong renormalization due to the spatial structure of the network. The global gain provides a convenient method to compare the effects of population structures on the extinction time.
The global gain G turns out to be the appropriate rescaled variable. When plotted with the global gain, the extinction time on the 2D grid collapse onto the same scaling function F + (x) as shown in Fig. 5. Because the instability of the ecosystem is wiped out by the grid structure, the other branch F − (x) does not show up. The reappearance of universal scaling indicates that the population structure only renormalizes the phase boundary but the underlying extinction behavior remains the same. Spatial and temporal correlations on the grid. To explore the different dynamics arisen from the network structure, it is helpful to study the spatial and temporal correlation functions, where the summation is over all species, X = A, B, C . The random variable X(r, t) = 1 when the species X is spotted at the distance r and in the later time t. Otherwise, its value is zero. The ensemble average X(r, t)X(0, 0) measures the probability of finding the same species with the space-time separation (r, t) and equals X(r, t) X(0, 0) if no correlation is present. The cross-species correlation functions are neglected here because they reflect the same trend as observed in C(r, t).
The spatial and temporal dependences of the correlation function is shown in Fig. 6. The spatial correlation is short-ranged with correlation length ξ ≈ 3.3 (lattice unit) for g = 1 . The data fits the usual exponential decay (19) C(r, t) ≡ X �X(r, t)X(0, 0)� − �X(r, t)��X(0, 0)�, Figure 3. Enhanced ecological stability on 2D grid. Differnt from the phase transition observed in the fullyconnected network (well-mixed population), the extinction time on the 2D grid always follows the exponential trend even for g 1 where the average local payoff g − 1 is negative. The spatial structure gives rise to nontrivial normalization on the local payoff matrix and enhances the biodiversity on the 2D grid. . Global payoff on the 2D grid. In spatial games, the local evolutionary dynamics is generated by the local payoff matrix. However, the global ecological stability is not directly reflected by the local payoff matrix due to the renormalization from the network structure. Referencing with the fully-connected network, the global gain G can be extracted by equating the characteristic constant of the extinction time. The mapping between the local gain g and the global gain G not only highlights the renormalization of the payoff matrix but also provides clear evidence for robust biodiversity ( G > 1 ) in 2D grid. www.nature.com/scientificreports/ rather well and implies no long-ranged pattern formation in the 2D grid. The temporal correlation is more subtle -it shows non-trivial oscillations ω ≈ 0.18 with a relatively long correlation time τ ≈ 37.8 . There exists clear oscillatory signature before the temporal correlation damps out. This oscillation in the steady state is not predicted in the replicative dynamics. Because the initial configuration is randomly sampled, the population size of each species is the same at the beginning. According to the replicator equations, the initial condition sits right on the fixed point and the dynamics is trivially static -no oscillation at all. One may cautiously argue that, due to statistical fluctuations, the initial condition is not exactly on the fixed point. Straightforward analysis leads to an intrinsic oscillation frequency ω 0 = (g + 1)/2 √ 3 near the fixed point. The sensitive dependence on the local gain g is not found in the numerics. Thus, the temporal oscillation must come from a different origin.

Scientific Reports
The situation is close to the nematic phase 73 where the correlation in one direction (temporal) is stronger than the other (spatial). One can picture the whole system is roughly divided into many spiraling units of linear size ξ . Each unit generates local spiral waves but never synchronizes into the global spiral due to interactions with the wavefronts of other neighboring units. Because the spiral unit is local, its oscillation cannot be captured by the usual replicative dynamics which works for uniformly mixed populations. This picture purveys a simple explanation why the biodiversity is greatly enhanced in the 2D grid. Suppose the probability for a spiral unit to go out is p. The interactions among the units ignite or suppress the neighboring ones. Thus, the probability for the complete extinction P ex requires all units to go out at the same time, P ex ≈ p N/ξ 2 , where N/ξ 2 is roughly the number of spiraling units. The extinction is proportional to the reciprocal of P ex , where N + = ξ 2 /| ln p| . The exponential trend of the extinction time arises from the concurrent decease of the spiraling units. Although complete understanding of the evolutionary dynamics on the 2D grid is not yet achieved, the snap shots in Fig. 6 supports the qualitative picture proposed here.

Discussion
Two major results are found in the agent-based Monte Carlo simulations. The first one is the emergence of universal scaling. The second one is the renormalized global payoff due to spatial structure of the network. While our numerical results clearly support the above conclusions, a rigorous proof remain open at the point of writing. In fact, some other numerical study 74 on the cooperation in social games also found that the rescaled extinction time follows universal scaling. However, their emphasis is in the neutral regime (critical phenomena) where the power law reigns. But, in the rock-paper-scissors game we studied here, we find the phase transition can be collapsed onto two universal curves when plotted in terms of the rescaled variables. This is quite unexpected and shows that the universal scaling works not just in the neutral regimes (or equivalently, the criticality). Thus, our conclusions are stronger and can serve as an important method to classify different dynamical trends for ecosystems with finite populations.
One can try to understand the extinction trends qualitatively by the dynamics of the biodiversity indicator χ(t) ≡ a(t)b(t)c(t) . When the population is infinite and well-mixed, its dynamics can be derived from the replicator equation, where Ŵ(a, b, c) = 1 − 3(ab + bc + ca) > 0 due to the constraint a + b + c = 1 for the positive-definite frequencies. It is clear that the drift flows toward extinction ( χ = 0 ) for g < 1 while the ecosystem is stable for g > 1 .
At the critical value g = 1 , the drift disappears and the flow is neutral. Now consider the effect of finite population. At the critical point g = 1 , the biodiversity indicator is no longer a constant of motion but suffers an O(1/N) drift toward extinction 34 . The exponential trend for g > 1 can be understood as a chain of successive rare events against the flow while the logarithmic trend for g < 1 is the direct consequence of the exponential decay due to the drift toward extinction. Right at the critical point g = 1 , the neutral drift implies the absence of any characteristic time scale and thus hints for the emergence of the universal scaling. Detail calculations for the fluctuation-induced drift leads to the extinction time T ex ∼ N and explains the exponent α = 1 . However, it remains an open question why the classification of the dynamics can be readily described by the scaling arguments developed for phase transitions in thermostatistics.
The 2D grid posts even more challenging puzzles. The success of applying global gain G to collapse all data onto the same scaling function is encouraging for further studies along the route of renormalization-group techniques. This shall provides an independent check on the relation between G and g, which is obtained numerically in our study here. The plausible coarse-graining in the time domain is further supported from the prolonged correlation time in the temporal correlation. The renormalization group may also provide more information about the exponents and the scaling functions.

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.