QUASI-STATIONARY OSCILLATIONS IN GAME-DRIVEN EVOLUTIONARY DYNAMICS

Nonlinear dynamics makes use of attactors to describe asymptotic behavior of complex systems. However, in real life can be unattainable, as achieving them might require time much exceeding all relevant timescales of a system. Therefore, there is an increasing interest in quasi-stationary states, where the system rapidly converges to and remains for a long time, before getting into an absorbing (asymptotic) state. Exemplifying in the fa-mous Dawkins‘ Battle of the Sexes game, we demon-strate that quasi-stationary distributions can produce not simply different, but a much more complex behavior, then the asymptotic ones, that is transient self-sustained oscillations of player numbers and the corresponding non-unimodal probability distribution. We ﬁnd that parameters of the quasi-stationary limit cycle depend on the population size.


Introduction
Nonlinear dynamics offers powerful tools to describe behavior of complex physical, biological, chemical or socio-ecological systems. If a system is dissipative, typically, one or another asymptotic regime, a stationary state, limit cycle or chaotic attractor is reached after a transient process. At the same time, in certain cases the transient process itself can take a substantial time. Sometimes, an asymptotic state may not be reached during a relevant timescale, for example, the coherence time of a quantum system, or a population lifetime [Biroli and Kurchan, 2001;Rose et al., 2016;Assaf and Mobilia, 2012]. In this case, one speaks about metastable or quasistationary states, rapidly converged on by a trajectory from generic initial conditions, and resided at for a long time [Rabinovich et al., 2008;Macieszczak et al., 2016]. Whether such states can be even more complex than an attractor, remains a challenging problem.
As a particular complex system we take an evolutionary "Battle of Sexes" game of populations that models the gender conflict over parental care [Dawkins, 1976]. In this game males and females have two alternative behavioral strategies, which stem from various expectations for courtship and different assistance in raising an offspring. On each step of the game a random player dies and a random player replicates itself, so that the population size remains constant. When all offsprings inherit a strategy from a parent of the same gender, an extinct strategy cannot return. We find that although such a state is absorbing, it can be unattainable over a very long time, if the population size is large. Employing the concept of the quasi-stationary distribution in absorbing Markov chains [Darroch et al., 1965;Collet et al., 2012], we find that the metastable dynamics is oscillatory. The paper is structured as follows. In Section 2 we introduce the model. In Section 3 and 4 its Moran process formulation and Markov chain model. In Section 5 we introduce the nonlinear dynamical model for the mean-field approximation. Section 6 outlines different numerical approaches to find quasistationary states. Section 7 makes the statement of the results and Section 8 concludes the paper.

Model
We consider two-player game "Battle of Sexes" (BoS) with finite size populations. Players A (males) have two available strategies s ∈ {1, 2} for reproduction. Similarly players B (females) have strategies s ∈ {1, 2}. Interaction between males and females is determined by some payoffs. The payoff matrix quantifies the reward received by a player after it has played against a member of opposite-sex population. The payoffs differ for males and females, as well as for combination of strategies, see Table 1.
Here a ss is the payoff for the male with strategy s ∈ {1, 2}, who plays against the female, who use strategy s ∈ {1, 2}. Similarly payoffs for the female b s s are defined. After the m-th game round the state of the populations is fully specified by the number of players using the first strategy: i (males) and j (females), 0 ≤ i, j ≤ N , where N is the size of populations.
The population split can be described by the number of players that stick to the first strategy: i (males) and j (females). Thus the state of the system can be expressed as an N ×N matrix p, which elements p i,j are the probabilities to find the system in the states i and j, respectively: The four states i ∈ {0, N }, j ∈ {0, N } are impossible to leave, and are called absorbing. Getting into one of it means that a populational strategy has degenerated and respective players have got extinct. The absorbing states are attractors of the evolutionary dynamics for any finite N, and the finite-size fluctuations will eventually lead a population to one of them. Ultimately, only one strategy survive in each of populations.
The states i ∈ {0, N }, j ∈ {1, . . . , N − 1} and i ∈ {1, . . . , N − 1}, j ∈ {0, N } are absorbing boundary states, a step to absorbing points. As soon as the population gets to the boundary, it will move in the direction of one of the two nearest absorbing states.
Other states are transient.

Moran Process
We use the Moran process to describe evolution of finite male and female populations in the BoS game [Moran, 1962;Nowak et al., 2004;Nowak, 2006]. The Moran process is a microscopic description of a birthdeath process. The average payoffs for males playing strategies s ∈ {1, 2} are Simirlarly, for females playing strategies s it follows Payoffs define the probabilities for members to be chosen for reproduction. For example, for the male populations these probabilities are /N is the average payoff for males, and ω ∈ [0, 1] is a baseline fitness. When ω = 0, the probability P A s (i, j) = 1/N is uniform across the population and does not depend on players. For the female population P B s (i, j) is determined in the similar way.
Essentially, the dynamics described by the Moran process consists of three steps: (i) an individual is selected for reproduction in accordance with probability (4), (ii) the selected player produces one identical offspring; (iii) the offspring replaces a randomly chosen member. This evolutionary mechanism acts in both populations, males and females, at the same time. An individual can not change its strategy and an offspring inherits the strategy of its parent. Therefore, the population size N is constant over time.

Markov Chain
At each time step the number i of males or j of females with first strategy can increase or decrease only by one or stay constant. Transition rates from the state i to states i + 1, i − 1 and i for male population A are expressed as [Traulsen et al., 2005] All other transitions are forbidden, therefore, the corresponding transition probabilities are equal to zero. Transition rates for female population B are defined analogously. The BoS game can be described by the Markov chain, that is stochastic transitions between states of the system [Behrends, 2000]. In our case the game dynamics can be evaluated by multiplying the state p with the transition fourth-order tensor S, with elements S(i, j, i , j ), where (i, j) is the current state and (i , j ) is a new one after a round of the Markov process. By using the bijections k = (N +1)i+j and l = (N +1)i +j , we can unfold the probability matrix p(i, j) into the vector p(k), and the fourth-order tensor S(i, j, i , j ) into the matrix S(k, l). This reduces the problem to a Markov chain p m+1 = S p m , where m is the number of the game round. The elements of the transition matrix S are determined using

Mean Field Model
The dynamics of infinite size populations can be described by a deterministic model that is written in terms of frequencies of the two (or more) strategies. In the continuous (mean-field) limit N → ∞ the dynamics of the variables x = i N and y = j N is defined by the adjusted replicator equations [Smith, 1982] x where In the asymptotic regime populations reach Nash equilibrium [Nash, 1950] of the game (x = 0.5, y = 0.5), that is, (i = N 2 , j = N 2 ), see Fig. 1. For finite size populations the mean-field model is only an approximation. The most drastic distinction is that absorbing states become asymptotically stable, and Nash equilibrium gets unstable. To take into account stochastic effects, one can intorduce noice and study Langevin Equations (10).

Finding the Quasi-stationary State
According to the definition in the above, by the quasistationary state q we will call the state, the system rapidly converges to and remains in for a long time, before getting into an absorbing state.
Below we describe three methods that we used to find and independently confirm the quasi-stationary state in the BoS game.

Darroch and Seneta approach
Following the method proposed in [Darroch et al., 1965], the quasi-stationary state q (9) can be found as the normalized right eigenvector that corresponds to the maximum eigenvalue λ of the reduced transition matrix R, that is obtained by striking out the rows and columns corresponding to absorbing boundaries: By inverse bijection, the vector q can be transformed into a two-dimensional probability density function (pdf) p QS . The resulting distribution p QS can be represented as a color 3D histogram. The advange of this method is that it returns the "numerically exact" solution (potentially, up to machine precision). However, it becomes unfeasible for large populations, because the algorithm implies the diagonalization of (N − 1) 2 × (N − 1) 2 matrix, and a computer quickly runs out of memory.
The two following methods allow to circumvent the problem and to approximate the quasi-stationary distribution by sampling it from random ensembles.

Numerical Simulation of the Moran Process
The simulation of the Moran process was performed by launching a random process many times from a random initial point (i, j), i, j ∈ {1, N − 1} and following the changes in the population structure over many rounds M . If the state remains unabsorbed after M rounds, it is sampled for the probability density function.

Langevin Equations
Another way to approximate a quasistationary state is to numerically solve the Langevin equations [Lemons et al., 1997;Gardiner, 2004]. The Langevin equation (10) is a stochastic differential equatioṅ Here x = i N , y = j N are continuous variables, ξ(t) represents uncorrelated Gaussian white noise.
The coefficients g can be represented by the 2 × 2 matrix The matrix G can be expressed in terms of the diffusion matrix Variables in (10) and (12) are defined using Kampen, 2007;Risken et al., 1996] where K, L ∈ {A, B}.
In order to obtain the matrix G we use the method outlined in [Traulsen et al., 2012]. We diagonalize the matrix D = U ΛU T and obtain the matrix G as where U is an eigenvectors matrix, Λ is an eigenvalues matrix of the matrix D.
Note that in the infinite population size limit, N → ∞, the stochastic terms in Eq.(10) vanish, and Langevin equations reduce to the mean-field equations Eq.(7).
The choice of symmetric payoff matrixes is convenient to keep all transition probabilities symmetric. We will see that this results in symmetric quasi-stationary distribution. The results for these specific payoff matrices are applicable to other payoff matrices if they have a similar structure.
We make use of the Darroch and Seneta approach to find the quasi-stationary distribution. The main result of the paper is that the found distributions are not unimodal, but they have a minimum in the center and a maximum around the it. An example of the quasistationary state for N = 200 is shown in Fig. 2. It suggests that the actual stochastic dynamics is governed by the metastable limit cycle about the Nash equilibrium (i * , j * ) = (i = N 2 , j = N 2 ). We perfomed the  Moran process and Langeven simulations for an independent validation, and found that approximate solutions give a good agreement with the numerically exact quasistationary state. Fig. 3 and Fig. 4 show the results of numerical simulations of the Moran process, unraveling dynamics on the quasistationary solution. It demonstrates oscillations in the number of players in both populations, the Nash equilibrium state being the center. In constast to the nonlinear mean-field solution (Fig.1), the trajectories diverge from the center. Notably, the trajectory makes many rotations on the metastable "limit cycle" before getting into the absorbing state, illustrating the nature of quasi-stationary dynamics.
Langevin simulations for a range of system sizes, N , reveal the size dependence of the metastable cycle radius R M C . It agrees with the theoretical expectation that as N → ∞, the mean-field model must be recovered, so that the distribution should become unimodal with the maximum in the Nash equilibrium, (i * , j * ), that is, (x = 1 2 , y = 1 2 ). Precisely, we find that the radius decreases with the increase of N , and scales to zero at R M C ∼ N −1/2 . In terms of the Langevin-oriented approach to the dynamics of finite populations, the emergence of the metastable limit cycle and its disappearance with N → ∞ can be viewed as a stochastic Hopf bifurcation [Arnold, 2003].

Conclusion
We studied quasi-stationary stats in the "Battle of Sexes" evolutionary game. We showed that under certain parameters the quasi-stationary distribution is nonunimodal, that corresponds to a metastable limit cycle, whereas the nonlinear mean-field approximation displays a stable equilibrium point. The radius of this cycle depends on the size of the population N . With increasing N , the cycle radius decreases and metastable state degenerates to a point, that is Nash equilibrium. This transition can be interpreted as a stochastic Hopf bifurcation. The results provide an insight into the quasi-stationary dynamics of finite-size populations, much more complex than either absorbing states or the mean field attractors.