A Finite Population Destroys a Traveling Wave in Spatial Replicator Dynamics

Christopher Griffin ¹¹1Applied Research Laboratory, Penn State University, University Park, PA 16802 ³³footnotemark: 3 Riley Mummah²²2Department of Ecology and Evolutionary Biology, UCLA, Los Angeles, CA 90095 Russ deForest ³³3Department of Mathematics, Penn State University, University Park, PA 16802

Abstract

We derive both the finite and infinite population spatial replicator dynamics as the fluid limit of a stochastic cellular automaton. The infinite population spatial replicator is identical to the model used by Vickers and our derivation justifies the addition of a diffusion to the replicator. The finite population form generalizes the results by Durett and Levin on finite spatial replicator games. We study the differences in the two equations as they pertain to a one-dimensional rock-paper-scissors game.

1 Introduction

Evolutionary games using the replicator dynamic have been studied extensively and are now well documented [1, 2, 3, 4, 5, 6]. Variations on the classical replicator dynamic include discrete time dynamics [7] and mutations [8, 9]. Additional evolutionary dynamics, such as imitation [1, 10, 11, 4] and exchange models [12] have been studied. Alternatively evolutionary games have been extended to include spatial models by a number of authors [13, 14, 15, 16, 17, 18, 19, 20]. Most of these papers append a spatial component to the classical replicator dynamics (see e.g., [18]) or discuss finite population replicator dynamics in which total population counts are used (see e.g., [13]). In the latter case, a spatial term is again appended to the classical replicator structure.

In [21], Durrett and Levin study discrete and spatial evolutionary game models and compare them to their continuous, aspatial analogs. For their study the authors focus on a specific class of two-player two-strategy games using a hawk-dove payoff matrix. Because their payoff matrix is $2\times 2$ , a single (spatial) variable $p(x,t)$ can be used to denote the proportion of the population playing hawk (Strategy 1) while a second variable $s(x,t)$ denotes the total population. Using the payoff matrix:

\mathbf{A}=\begin{bmatrix}a&b\\ c&d\end{bmatrix},

the remarkable reaction-diffusion equation is analyzed:

	$\displaystyle\frac{\partial p}{\partial t}$	$\displaystyle=\Delta p+\frac{2}{s}\nabla p\cdot\nabla s+pq((a-c)p+(b-d)(1-p))$		(1)
	$\displaystyle\frac{\partial s}{\partial t}$	$\displaystyle=\Delta s+s\left(ap^{2}+(b+c)p(1-p)+(1-p)^{2}d\right)-\kappa s^{2}.$		(2)

Here $\kappa$ is a death rate due to overcrowding and $q=(1-p)$ . Let $\mathbf{e}_{i}\in\mathbb{R}^{n}$ be the $i^{\text{th}}$ unit vector. Let $\mathbf{u}=\langle{p,q}\rangle$ . When $\kappa=0$ , we can rewrite these equation as:

	$\displaystyle\frac{\partial p}{\partial t}$	$\displaystyle=\Delta p+\frac{2}{s}\nabla p\cdot\nabla s+p\left(\mathbf{e}_{1}^{T}-\mathbf{u}^{T}\right)\mathbf{A}\mathbf{u}$		(3)
	$\displaystyle\frac{\partial s}{\partial t}$	$\displaystyle=\Delta s+s\cdot\mathbf{u}^{T}\mathbf{A}\mathbf{u}.$		(4)

That is, Durrett and Levin have encoded the replicator dynamic into a finite population spatial partial differential equation, which differs from the one used by Vickers [18] because it assumes a finite population in its derivation. The second term in Eq. 3 is an advection of species $p$ along a velocity vector given by the gradient of the whole population. However, this advection speed is inversely proportional to the population size. Thus small populations can have a dramatic effect on this term. As we show in Section 3, this behavior carries through for general games. To illustrate the advection clearly, we can define the potential function:

\phi(x)=2\log\left[s(x,t)\right].

We then have:

\mathbf{v}=\nabla\phi=\frac{2}{s}\nabla s.

Eq. 3 can then be written as:

\frac{\partial p}{\partial t}=\Delta p+\mathbf{v}\cdot\nabla p+p\left(\mathbf{e}_{1}^{T}-\mathbf{u}^{T}\right)\mathbf{A}\mathbf{u},

giving a standard advection-diffusion equation. The logarithm in the potential function conveniently yields the per-capita advection rate.

In this paper we show that Durrett and Levin’s finite population spatial replicator is generalizable to an arbitrary payoff matrix. We then focus our attention on the one-dimensional rock-paper-scissors game under the replicator dynamics in both finite and infinite populations. This system has interesting properties in both the finite and infinite population cases. In particular, we show: (i) the model used by Vicker’s [18] arises naturally as the infinite population limit of the generalization of Durrett and Levin’s model, which in turn can be derived from a stochastic cellular automaton (particle) model as a fluid limit, and (ii) the one dimensional infinite population spatial rock-paper-scissors dynamic has a constant amplitude traveling wave solutions in biased RPS. However, the finite population version does not exhibit such solutions, but does seem to exhibit attracting stationary solutions. We illustrate the latter result numerically.

2 Model

Let $\mathbf{A}\in\mathbb{R}^{n\times n}$ be a payoff matrix for a symmetric game [4]. All vectors are column vectors unless otherwise noted. Below we construct a stochastic cellular automaton and show that the fluid limit of this system yields a generalization of Durrett and Levin’s specific finite population spatial replicator.

The state of cell $i$ at time index $k$ of the cellular automaton is a tuple $\langle{U_{1}(i,k),\dots,U_{n}(i,k)}\rangle$ where $U_{j}(i,k)$ provides the size of the population of species $j$ at position $i$ at time $k$ . For simplicity, we assume that species interaction may only happen between cells and not within cells; i.e., the $U_{1}(i,k)$ members of species $1$ will not play against the $U_{n}(i,k)$ members of species $n$ . This assumption will become irrelevant in the limit.

During state update, an agent $A$ at cell $i$ chooses a random direction (cell $i^{\prime}$ ) and a random member of the population (agent $A^{\prime}$ ) within that cell. Assume Agent $A$ uses strategy $r$ while Agent $A^{\prime}$ uses strategy $s$ . After play, there are $\alpha\cdot\mathbf{A}_{rs}$ additional agents at cell $i$ using strategy $r$ and $\beta\cdot\mathbf{A}_{rs}$ additional agents playing strategy $r$ at cell $i^{\prime}$ , where $\alpha+\beta=1$ is the probability of motion from cell $i$ to cell $i^{\prime}$ . If $\mathbf{A}_{rs}<0$ , then agents are removed from their respective cells. To avoid computational issues when more members of a species die than are present, $\mathbf{A}$ can be modified so that $\mathbf{A}_{rs}\geq 0$ for all $r,s\in\{1,\dots,n\}$ , without altering the evolutionary dynamics in proportion [2]. The update rule for a single agent is illustrated in Fig. 1.

Refer to caption — Figure 1: A single interaction between two agent in a spatial game on a one-dimensional lattice is illustrated. The direction of interaction is chosen at random. Motion is random and governed by $\alpha$ and $\beta$ where $\alpha+\beta=1$ .

The process described above and illustrated in Fig. 1 is assumed to be happening simultaneously for each agent and we assume that the replication/death as a result of game play along with the migration are happening on (roughly) the same time scale. Using these assumptions, we can construct mean-field equations for the population $U_{r}$ at position $i$ and time $k+1$ in a 1-D cellular automaton, assuming an equal likelihood that agents diffuse left or right. The mean number of agents playing strategy $r$ present at position $i$ at time $k+1$ is determined by:

1.

The expected number of agents who remain at cell $i$ : $\alpha U_{r}(i,k)$
2.

The expected number of new agents created at cell $i$ who remain at cell $i$ :

$\frac{\alpha}{2}U_{r}(i,k)\left(\sum_{s}\mathbf{A}_{rs}\left(u_{s}(i+1,k)+u_{s}(i-1,k)\right)\right).$
3.

The expected number of agents who migrate to position $i$ from neighboring cells:

$\frac{\beta}{2}\left(U_{r}(i+1,k)+U_{r}(i-1,k)\right).$
4.

The expected number of agents created in a neighboring cell who migrate to cell $i$ :

$\frac{\beta}{2}\left(U_{r}(i+1,k)+U_{r}(i-1,k)\right)\sum_{s}\mathbf{A}_{rs}u_{s}(i,k).$

Here $u_{s}(i,k)$ is the proportion of the population at cell $i$ playing strategy $s$ at time step $k$ .

Let $\mathbf{u}(i,k)=\langle{u_{1}(i,k),\dots,u_{n}(i,k)}\rangle$ . Re-writing sums as matrix products, the expected number of agents playing strategy $r$ at cell $i$ at time $k+1$ is:

U_{r}(i,k+1)=\alpha U_{r}(i,k)+\frac{\beta}{2}\left(U_{r}(i+1,k)+U_{r}(i-1,k)\right)+\\ \frac{\alpha}{2}U_{r}(i,k)\left(\mathbf{e}_{r}^{T}\mathbf{A}\left(\mathbf{u}(i+1,k)+\mathbf{u}(i-1,k)\right)\right)+\\ \frac{\beta}{2}\left(U_{r}(i+1,k)+U_{r}(i-1,k)\right)\mathbf{e}_{r}^{T}\mathbf{A}\mathbf{u}(i,k).

(5)

Assume the cellular grid has lattice spacing $\Delta x$ . Following [22] and using a Taylor approximation, we can write:

	$\displaystyle U_{r}(x,t+\Delta t)$	$\displaystyle\approx U_{r}(x,t)+\Delta t\frac{\partial U_{r}(x,t)}{\partial t}+O(\Delta t^{2})$
	$\displaystyle U_{r}(x+\Delta x,t)$	$\displaystyle\approx\sum_{j=0}^{2}\frac{\Delta x^{j}}{j!}\frac{\partial^{j}U_{r}(x,t)}{\partial x^{j}}+O(\Delta x^{3}).$

3 Derivation of Fluid Limits

We proceed to derive the mean-field approximation. Passing to the continuous case and assuming that interaction rates decrease linearly with $\Delta t$ , we can write a second order approximation of Eq. 5 as:

\Delta t\frac{\partial U_{r}(x,t)}{\partial t}=\alpha U_{r}(x,t)+\beta\left(U_{r}(x,t)+\frac{1}{2}\Delta x^{2}\frac{\partial^{2}U_{r}}{\partial x^{2}}\right)+\\ \alpha\Delta tU_{r}(x,t)\cdot\mathbf{e}_{r}^{T}\mathbf{A}\left(\mathbf{u}(x,t)+\frac{1}{2}\Delta x^{2}\frac{\partial^{2}\mathbf{u}}{\partial x^{2}}\right)+\\ \beta\Delta t\left(U_{r}(x,t)+\frac{1}{2}\Delta x^{2}\frac{\partial^{2}U_{r}}{\partial x^{2}}\right)\mathbf{e}_{r}^{T}\mathbf{A}\mathbf{u}(x,t)-U_{r}(x,t).

(6)

Where the $-U_{r}(x,t)$ on the right-hand-side arises from the formation of the Newton quotient on the left-hand-side. Expanding and simplifying yields:

\Delta t\frac{\partial U_{r}(x,t)}{\partial t}=\frac{\beta}{2}\Delta x^{2}\frac{\partial^{2}U_{r}}{\partial x^{2}}+\Delta tU_{r}(x,t)\mathbf{e}_{r}^{T}\mathbf{A}\mathbf{u}(x,t)+\\ \frac{\alpha\Delta t}{2}\Delta x^{2}U_{r}(x,t)\mathbf{e}_{r}^{T}\mathbf{A}\frac{\partial^{2}\mathbf{u}}{\partial x^{2}}+\\ \frac{\beta\Delta t}{2}\Delta x^{2}\frac{\partial^{2}U_{r}}{\partial x^{2}}\mathbf{e}_{r}^{T}\mathbf{A}\mathbf{u}(x,t).

(7)

Assume $\beta\in(0,1]$ . Dividing through by $\Delta t$ and assuming that $\lim_{\Delta t\rightarrow 0}\Delta x^{2}/\Delta t=2D/\beta$ yields:

\frac{\partial U_{r}(x,t)}{\partial t}=U_{r}(x,t)\mathbf{e}_{r}^{T}\mathbf{A}\mathbf{u}(x,t)+D\frac{\partial^{2}U_{r}}{\partial x^{2}}.

(8)

The constant $D$ is the diffusion constant and the assumption that

\lim_{\Delta t\rightarrow 0}\Delta x^{2}/\Delta t=2D/\beta

is a variant of the assumption used to derive Fick’s Law [23] and identical when $\beta=1$ .

These are the spatial dynamics used by Durrett and Levin, (in the first part of their paper), but are derived only by adding a diffusion term to the standard finite population growth equations. In [21], Durrett and Levin note that they derive a set of equations they feel are more appropriate for modeling finite spatial systems. Their derivation at the end of [21] (for a specific hawk-dove system) rests on the assumption that migration happens “on a much faster timescale” than game interactions. Our model assumes that migration and game interactions occur on approximately the same time scale. Under this assumption, Eq. 8 is the correct spatial adaptation for finite populations; i.e., one simply adds a diffusion term. In the case where migration happens more quickly, then the derivation in [21] should be used instead.

4 Spatial Replicator with Finite Population

The derivation of Eqs. 1 and 2 are not given in [21]. They can be generalized for an arbitrary evolutionary game using Eq. 8 as the starting point. Let:

M(x,t)=\sum_{s}U_{s}(x,t).

Differentiating we have:

\frac{\partial}{\partial t}\frac{U_{r}(x,t)}{M(x,t)}=\frac{1}{M(x,t)}\frac{\partial U_{r}(x,t)}{\partial t}-\\ u_{r}(x,t)\sum_{s}\frac{1}{M(x,t)}\frac{\partial U_{s}(x,t)}{\partial t}.

Substituting from Eq. 8 we obtain:

\frac{\partial u_{r}}{\partial t}=u_{r}\cdot\left(\mathbf{e}_{r}^{T}\mathbf{A}\mathbf{u}-\mathbf{u}^{T}\mathbf{A}\mathbf{u}\right)+\\ \frac{D}{M}\left(\frac{\partial^{2}U_{r}}{\partial x^{2}}-u_{r}\cdot\sum_{s}\frac{\partial^{2}U_{s}}{\partial x^{2}}\right).

(9)

Unlike in the derivation of the standard replicator dynamic, the rate of change of the population proportion is not solely a function of the proportions themselves.

We can remove dependence on the individual populations to derive an independent (coupled) system of differential equations that includes only the total population. For arbitrary strategy $r$ , we can apply the quotient rule to obtain:

M\frac{\partial u_{r}}{\partial x}=\frac{\partial U_{r}}{\partial x}-u_{r}\frac{\partial M}{\partial x}

Differentiating again, multiplying by $1/M$ and re-arranging yields the expression:

\frac{1}{M}\frac{\partial^{2}U_{r}}{\partial x^{2}}=\frac{u_{r}}{M}\frac{\partial^{2}M}{\partial x^{2}}+\frac{2}{M}\frac{\partial u_{r}}{\partial x}\frac{\partial M}{\partial x}+\frac{\partial^{2}u_{r}}{\partial x^{2}}.

Using this we can write

\frac{D}{M}\left(\frac{\partial^{2}U_{r}}{\partial x^{2}}-u_{r}\sum_{s}\frac{\partial^{2}U_{s}}{\partial x^{2}}\right)=D\left(\frac{2}{M}\frac{\partial M}{\partial x}\frac{\partial u_{r}}{\partial x}+\frac{\partial^{2}u_{r}}{\partial x^{2}}\right)

Using this we can re-write Eq. 9 as:

\frac{\partial u_{r}}{\partial t}=u_{r}\cdot\left(\mathbf{e}_{r}^{T}\mathbf{A}\mathbf{u}-\mathbf{u}^{T}\mathbf{A}\mathbf{u}\right)+\\ D\left(\frac{2}{M}\frac{\partial M}{\partial x}\frac{\partial u_{r}}{\partial x}+\frac{\partial^{2}u_{r}}{\partial x^{2}}\right).

(10)

The dynamics of $M$ can be derived (by addition) from Eq. 8:

\frac{\partial M}{\partial t}=M\mathbf{u}^{T}\mathbf{A}\mathbf{u}+D\frac{\partial^{2}M}{\partial x^{2}}.

Thus, we have a coupled set of differential equations written entirely in terms of $\mathbf{u}$ and $M$ , rather than $U_{r}$ , $M$ and $\mathbf{u}$ :

\left\{\begin{aligned} &\forall r\left\{\begin{aligned} \frac{\partial u_{r}}{\partial t}=&u_{r}\cdot\left(\mathbf{e}_{r}^{T}\mathbf{A}\mathbf{u}-\mathbf{u}^{T}\mathbf{A}\mathbf{u}\right)+\\ &\hskip 20.00003ptD\left(\frac{2}{M}\frac{\partial M}{\partial x}\frac{\partial u_{r}}{\partial x}+\frac{\partial^{2}u_{r}}{\partial x^{2}}\right)\end{aligned}\right.\\ &\hskip 20.00003pt\frac{\partial M}{\partial t}=M\mathbf{u}^{T}\mathbf{A}\mathbf{u}+D\frac{\partial^{2}M}{\partial x^{2}}\end{aligned}\right..

(11)

This is the spatial replicator equation for finite populations. Letting $M=s$ and $D=1$ , we recover the dynamics of Durrett and Levin. In contrast to the aspatial replicator, the inclusion of dynamics for $M$ yields a linearly independent system of differential equations.

Allowing $M$ to approach infinity uniformly in $x$ , we arrive at the fluid limit in terms of $\mathbf{u}$ alone; this is the 1D nonlinear reaction-diffusion equation used by Vicker’s [18, 19]:

\forall r\left\{\begin{aligned} \frac{\partial u_{r}}{\partial t}&=u_{r}\cdot\left(\mathbf{e}_{r}^{T}\mathbf{A}\mathbf{u}-\mathbf{u}^{T}\mathbf{A}\mathbf{u}\right)+D\frac{\partial^{2}u_{r}}{\partial x^{2}}\end{aligned}\right..

(12)

Generalization to $N$ -dimensions is straightforward by replacing $\partial^{2}_{x}$ with the Laplacian $\Delta$ . The $N$ -dimensional spatial replicator with finite population is given by:

	$\displaystyle\forall r\left\{\begin{aligned} \frac{\partial u_{r}}{\partial t}=&u_{r}\cdot\left(\mathbf{e}_{r}^{T}\mathbf{A}\mathbf{u}-\mathbf{u}^{T}\mathbf{A}\mathbf{u}\right)+\left(\frac{2}{M}\nabla M\cdot\nabla u_{r}+\Delta u_{r}\right)\end{aligned}\right.$
	$\displaystyle\hskip 20.00003pt\frac{\partial M}{\partial t}=M\mathbf{u}^{T}\mathbf{A}\mathbf{u}+D\Delta M.$

Thus, the finite population case adds a nonlinear advection term that forces $u_{r}$ to follow the population gradient. A similar system is studied by deForest and Belmonte in [13], where the payoff gradient is followed instead of the population gradient.

In both Eqs. 11 and 12, we see that the aspatial replicator dynamics appear on the right hand side perturbed by a spatial term. It is well known that the dynamics of the aspatial replicator are confined to the $n$ -dimensional simplex $\Delta_{n}$ . This remains true for the spatial replicator dynamics with finite populations. Moreover, the solution $u_{r}=1$ (i.e., there is only one population) is a fixed point for the spatial replicator dynamic since the spatial derivative of the probability distribution of the population proportions is zero and the time derivative is identically zero as expected. Thus, pure populations are constant stationary solutions for these dynamics. Lastly, if $\tilde{\mathbf{u}}(x,t)$ is a constant solution at a Nash equilibrium for the game defined by $\mathbf{A}$ , then the right-hand-side is again identically zero by the Folk Theorem [4] of evolutionary game theory together with the fact that there is no spatial variation. Thus every Nash equilibrium of the matrix game corresponds to a spatially constant stationary solution of the spatial replicator dynamic in both the finite and infinite population cases.

5 Example: One Dimensional Rock-Paper-Scissors

Durrett and Levin’s analysis of Hawk-Dove was aided by the fact that one strategy can be eliminated, leaving a coupled system of two partial differential equations. In the remainder of this paper, we analyze variations of rock-paper-scissors (RPS), which yield more interesting results because of its cyclic three-strategy nature and because it can be easily parameterized as discussed in [2]. Substantial work has been done on (spatial) RPS without assuming a replicator dynamic or assuming a general replicator dynamic [24, 25, 26, 27, 28, 29, 30, 31, 32, 33]. In an example of a very recent generalization, Kabir & Tanimoto study pairwise evolution in RPS with noise [34].

A major focus of work in non-replicator spatial RPS has been on identifying spatial structures (e.g., spirals) that can emerge from these dynamics. For completeness, we quantify the relationship between the work in [26, 27, 28, 29, 30, 31, 32, 33] and the replicator dynamic considered in this work in Appendix A. There, we show that the equation system used in these references does not subsume results on the replicator dynamic in general and hence one cannot extend these results automatically to this case.

The objective in the remainder of this work is to compare the finite and infinite population replicator for a more interesting three strategy game. For simplicity (and in contrast to much of the previous work in spatial RPS), we focus our attention on the one-dimensional PDE.

We consider a generalized RPS payoff matrix as given by:

\mathbf{A}=\begin{bmatrix}0&-1&1+a\\ 1+a&0&-1\\ -1&1+a&0\end{bmatrix}.

When $a=0$ , this is the standard RPS game which has Nash equilibrium $\langle{\tfrac{1}{3},\tfrac{1}{3},\tfrac{1}{3}}\rangle$ . This is the unique interior fixed point and the aspatial replicator exhibits an elliptic fixed point at this Nash equilibrium. This Nash equilibrium is preserved for $a\neq 0$ ; and corresponds to an asymptotically stable interior fixed point when $a>0$ and unstable fixed point when $a<0$ [35] as a result of a degenerate Hopf bifurcation. This degenerate Hopf bifurcation is well understood and leads to traveling wave solutions, which we discuss in Section 5.1. We note $\mathbf{A}$ is not as fully general as the RPS matrix given in [4] or [35], but it is substantially easier to work with.

Let $\mathbf{u}=\langle{u_{r},u_{p},u_{s}}\rangle$ and note:

\zeta(u_{r},u_{p},u_{s})\overset{\Delta}{=}\mathbf{u}^{T}\mathbf{A}\mathbf{u}=a\left(u_{r}u_{s}+u_{r}u_{p}+u_{s}u_{p}\right).

There are at least two classes of global solutions to Eqs. 11 and 12:

Equilibrium Solution

Here, $u_{r}=u_{p}=u_{s}=\tfrac{1}{3}$ and $M$ solves:

0=\frac{a}{3}M+M^{\prime\prime}

or there is a single population (e.g., $u_{r}=1$ ) and $M$ satisfies $M^{\prime\prime}=0$ .

Oscillating Solution

Here, $u_{*}(x,t)\equiv\upsilon_{*}(t)$ with $*\in\{r,p,s\}$ where $\upsilon_{r}$ , $\upsilon_{p}$ , $\upsilon_{s}$ are solutions to the standard RPS replicator and $M(x,t)=\mu(t)$ satisfies linear equation:

\dot{\mu}=\zeta(\upsilon_{r},\upsilon_{p},\upsilon_{s})\mu.

In either case, we may impose appropriate boundary conditions.

5.1 Traveling Wave Solutions

Consider the infinite population model. Letting $z=x-ct$ , we can re-write Eq. 12 in compact form as:

\forall i\left\{\begin{aligned} Dv_{i}^{\prime}&=-u_{i}\left(\mathbf{e}_{i}^{T}\mathbf{A}\mathbf{u}-\mathbf{u}^{T}\mathbf{A}\mathbf{u}\right)-cv_{i}\\ u_{i}^{\prime}&=v_{i}\end{aligned}\right..

(13)

For RPS we have the six dimensional non-linear system:

\left\{\begin{aligned} Dv_{r}^{\prime}&=au_{r}\zeta(u_{r},u_{s},u_{p})-u_{r}(u_{s}-u_{p})-au_{r}u_{s}-cv_{r}\\ u_{r}^{\prime}&=v_{r}\\ Dv_{p}^{\prime}&=au_{p}\zeta(u_{r},u_{s},u_{p})-u_{p}(u_{r}-u_{s})-au_{p}u_{r}-cv_{p}\\ u_{p}^{\prime}&=v_{p}\\ Dv_{s}^{\prime}&=au_{s}\zeta(u_{r},u_{s},u_{p})-u_{s}(u_{p}-u_{r})-au_{s}u_{p}-cv_{s}\\ u_{s}^{\prime}&=v_{s}.\end{aligned}\right.

(14)

If $\mathbf{u}^{*}$ is a Nash equilibrium of $\mathbf{A}$ , then the pair $\mathbf{u}=\mathbf{u}^{*}$ and $\mathbf{v}=\mathbf{0}$ is a fixed point of Eq. 13. Linearizing about $u_{r}=u_{p}=u_{s}=\tfrac{1}{3}$ , $v_{r}=v_{p}=v_{s}=0$ , we obtain the eigenvalues:

	$\displaystyle\lambda_{1,2}$	$\displaystyle=\frac{-3c\pm\sqrt{9c^{2}+12aD}}{6D}$
	$\displaystyle\lambda_{3,4}$	$\displaystyle=\frac{-3c\pm\sqrt{9c^{2}+6aD+6D\sqrt{3}\sqrt{-(a+2)^{2}}}}{6D}$
	$\displaystyle\lambda_{5,6}$	$\displaystyle=\frac{-3c\pm\sqrt{9c^{2}+6aD-6D\sqrt{3}\sqrt{-(a+2)^{2}}}}{6D}.$

We can simultaneously show that for appropriate choice of wave speed, a Hopf bifurcation and hence a two dimensional center manifold exists and therefore a non-decaying traveling wave solution exists for the PDE. As a by-product, we compute the wave speed for a non-decaying traveling wave in terms of $a$ and $D$ . For some constant $b$ (to be determined), let:

(3c\pm bi)^{2}=9c^{2}+6aD\pm 6D\sqrt{3}\sqrt{-(a+2)^{2}}=\\ 9c^{2}+i6aD\pm 6D(a+2)\sqrt{3}.

Expanding the left hand side and relating real and imaginary parts we have:

	$\displaystyle 9c^{2}-b^{2}=9c^{2}+6aD$
	$\displaystyle 6bc=6D(a+2)\sqrt{3}.$

Solving for $b$ and $c$ yields:

	$\displaystyle b$	$\displaystyle=\sqrt{-6aD}$
	$\displaystyle\pm\tilde{c}$	$\displaystyle=\pm\frac{(a+2)\sqrt{2k}}{2\sqrt{-a}}$

Without loss of generality, assume positive $\tilde{c}$ . Using this information, we can rewrite the eigenvalues as:

	$\displaystyle\tilde{\lambda}_{1,2}$	$\displaystyle=\frac{-3\tilde{c}\pm\sqrt{9\tilde{c}^{2}+12aD}}{6D}$
	$\displaystyle\tilde{\lambda}_{3}$	$\displaystyle=\frac{-6\tilde{c}-bi}{6D}$
	$\displaystyle\tilde{\lambda}_{4}$	$\displaystyle=\frac{bi}{6D}$
	$\displaystyle\tilde{\lambda}_{5}$	$\displaystyle=\frac{-6\tilde{c}+bi}{6D}$
	$\displaystyle\tilde{\lambda}_{6}$	$\displaystyle=\frac{-bi}{6D}.$

For this solution to be physically realized, we must have $a<0$ . We also assume $a>-2$ or the dynamics changes (i.e., winning becomes losing). These requirements and our assumption that $\tilde{c}>0$ implies that $\mathrm{Re}(\lambda_{1,2})<0$ for all choices of $D>0$ and $a\in(-2,0)$ . Therefore, this system has a four dimensional stable manifold and two pure imaginary eigenvalues, which satisfies the first requirement of Hopf’s theorem (see [36], Page 152).

Consider $\lambda_{4,6}$ as functions of $c$ with:

\lambda_{4,6}(c)=\frac{-3c+\sqrt{9c^{2}+6aD\pm 6D\sqrt{3}\sqrt{-(a+2)^{2}}}}{6D}.

For $c=\tilde{c}>0$ , we know that $\mathrm{Re}[\lambda_{4,6}(\tilde{c})]=0$ . Differentiating we have:

\lambda_{4,6}^{\prime}(c)=-\frac{1}{2D}+\frac{3c}{2D\sqrt{9c^{2}+6aD\pm 6D\sqrt{3}\sqrt{-(a+2)^{2}}}}.

Evaluating at $c=\tilde{c}$ and simpilfying yields:

\lambda_{4,6}^{\prime}(\tilde{c})=-\frac{1}{2D}+\frac{3\tilde{c}\left(3\tilde{c}\mp bi\right)}{2D\left(9\tilde{c}^{2}-b^{2}\right)},

by choice of $\tilde{c}$ . We conclude:

\mathrm{Re}\left[\lambda_{4,6}^{\prime}(\tilde{c})\right]=-\frac{1}{2D}+\frac{9\tilde{c}^{2}}{2D\left(9\tilde{c}^{2}-b^{2}\right)}=\frac{b^{2}}{2D\left(9\tilde{c}^{2}-b^{2}\right)}.

This is non-zero since we assume $a<0$ and $D>0$ . Thus the real parts of eigenvalues $\lambda_{4,6}$ cross the imaginary axis with non-zero speed, satisfying the second criterion of Hopf’s theorem. Thus we conclude that the six dimensional traveling wave ODE has a solution and moreover exhibits a Hopf bifurcation, implying the existence of traveling wave solutions. In our numeric simulations, we show that fine tuning the parameters leads to a numerically stable traveling wave over the region of integration when $a=-\tfrac{4}{5}$ and $D=\tfrac{1}{12}$ (see Figs. 4 and 6).

5.2 Illustrative Behavior in the Unbiased Games

Consider Eqs. 11 and 12 and assume the periodic boundary conditions $u_{*}(-\pi,t)=u_{*}(\pi,t)$ for all time. When $a=0$ , then $\zeta(u_{r},u_{p},u_{s})=0$ , so the population $M(x,t)$ satisfies the heat equation. For simplicity, we choose a solution to the heat equation that models the diffusion of a population outward:

M(x,t)=\frac{4000e^{-\frac{x^{2}}{4\beta(t+10)}}}{\sqrt{\pi}\sqrt{\beta(t+10)}}.

The initial conditions:

	$\displaystyle u_{r}(x,0)$	$\displaystyle=\frac{1}{3}\left[1+\sin\left(x-\tfrac{4\pi}{3}\right)\right]$
	$\displaystyle u_{p}(x,0)$	$\displaystyle=\frac{1}{3}\left[1+\sin\left(x-\tfrac{2\pi}{3}\right)\right]$
	$\displaystyle u_{s}(x,0)$	$\displaystyle=\frac{1}{3}\left[1+\sin\left(x\right)\right]$

model three populations that are proportionally spread around a circle. The behavior of the three population proportions in both finite and infinite populations are shown in Fig. 2.

In the infinite population case we observe that the three population proportions converge toward the equilibrium stationary solution. This is further illustrated in Fig. 3 (left) where we show a ternary plot of the three strategies at $x=0$ . By contrast the finite population solution converges toward the oscillating stationary solution, as illustrated in Fig. 3 (right). We show the corresponding solution to the ordinary RPS replicator dynamic to which the system converges at all points in space. The convergence of the infinite population system to the equilibrium at all points in space is illustrated in Fig. 3 (left).

Before converging to an oscillating solution, we can see the one-dimensional proportions proportions are affected by the diffusion of the population at large. In Fig. 2, we note that the finite population plots are stretched with respect to their infinite population counterparts. This is particularly noticeable at time $t=5$ and $t=10$ . This stretching, caused by the advection of the total population, leads to the difference in steady state solutions for the same initial conditions.

We contrast this behavior with the case when $a<0$ . Here, the population will collapse since $\zeta(u_{r},u_{p},u_{s})\leq 0$ at all times. As noted, we set $a=-\tfrac{4}{5}$ and $D=\tfrac{1}{12}$ . From this we compute a constant amplitude traveling wave speed of:

\tilde{c}=\frac{1}{2}\sqrt{\frac{3}{10}}.

This traveling wave can be seen in Fig. 4 in the infinite population plots.

In contrast the population collapse (shown in Fig. 5) causes the finite population proportions to converge to an oscillating stationary strategy with greater and greater amplitude. This is precisely the behavior one expects to see from RPS under the replicator dynamics when $a<0$ . This is also shown Fig. 6(right), in which we show an example RPS trajectory with $a=-\tfrac{4}{5}$ and example trajectories at various points in space.

Interestingly, the population collapse slows dramatically as $t$ increases. This is caused by the fact that for (nearly) pure strategies, $\zeta(u_{r},u_{s},u_{t})\approx 0$ . This is illustrated in Fig. 5, where the nearly exponential decay flattens after $t\approx 30$ . Additionally, Fig. 5 shows that the population becomes asymmetric as the system evolves, as a result of the various species dynamics.

Convergence to the stable traveling wave solution is illustrated in Fig. 6 (left). Using the computed $\tilde{c}$ , numerical analysis provided initial conditions which can be used to find a numerical solution to Eq. 14 that produces the closed curve (on the center manifold) that is guaranteed to exist by the eigenvalue analysis performed in the previous section. This is shown in Fig. 6 (left).

6 Conclusion

In this paper we studied a finite and infinite population spatial replicator. We showed how the finite population spatial replicator can be derived from first principles from a stochastic cellular automaton model and from there how the infinite population replicator used by Vickers [18, 19] follows from this. This result generalizes the work of Durrett and Levin [21] who first derived and studied the finite population spatial replicator for a specific game. We then compared the finite and infinite population spatial replicator for rock-paper-scissors on a circle ( $S^{1}$ ). Our results are consistent with the work in [29, 30, 31, 33, 27, 28, 32, 26], which studies various characteristics of RPS, not necessarily in connection with spatial replicator. Consistent with the work in [32, 26] we show that for a certain rock-paper-scissors variant stable amplitude traveling waves can emerge as solutions to the infinite population spatial replicator by proving the existence of a Hopf bifurcation in the traveling wave ODE. These traveling waves are destroyed by population collapse in the finite population spatial replicator.

The finite population spatial replicator is intriguing because it is a highly non-linear reaction-advection-diffusion equation where advection is governed by the per capita bulk population motion. Identifying cases where complex behaviors emerge in the finite population case is a clear future direction. Additionally, studying stationary solutions may yield insights. For example, in our unbiased RPS, solutions to the stationary population equation are just the harmonic functions. Using this simplification may help identify interesting properties of the population proportion equations.

Acknowledgement

Portions of CG’s work were supported by the National Science Foundation Grant CMMI-1932991.

Appendix A Relation to Other Rock Paper Scissors Games

A substantial amount of work has been done on rock-paper-scissors outside of the context of the replicator dynamic. This work does not map conveniently to results on the replicator dynamic or its spatial variants [26, 25, 27, 28, 29, 24, 30, 31, 32, 33, 34]. The earliest work to study competition among three species with cyclic dominance may be [24], which is contemporary but does predate the earliest work in evolutionary game theory (see e.g., [37, 35, 38, 39]). In the past twenty years, there has been substantial work on the spatial dynamics of RPS that is independent of the spatial replicator dynamic [13, 14, 15, 16, 17, 18, 19, 20]. Work till 2014 is reviewed in [28]. Mobility in cyclic competition (RPS) is studied in [29, 30]. The impact of reaction rates on spatial RPS is considered in [31], while the emergence of spiraling waves is studied in [33, 27]. More recently, traveling waves, spirals and heteroclininc bifurcations have been studied in [32, 26]. A discrete time model displaying chaos is considered in [40].

The results most closely related to those in this paper (specifically Section 5.1) can be found in the pair of papers by Postlethwaite and Rucklidge [32, 26]. Both these papers study a traveling wave solution for a specific spatial model of the rock-paper-scissors, with a more formal treatment given in [26]. The motivating aspatial model is given by the system of equations:

$\displaystyle\dot{a}$	$\displaystyle=a(1-(a+b+c)-(\sigma+\zeta)b+\zeta c)$	(15)
$\displaystyle\dot{b}$	$\displaystyle=b(1-(a+b+c)-(\sigma+\zeta)c+\zeta a)$	(16)
$\displaystyle\dot{c}$	$\displaystyle=c(1-(a+b+c)-(\sigma+\zeta)a+\zeta b)$	(17)

The observation is made that when $\sigma=0$ , this system exhibits the conserved quantities $a+b+c=1$ and $abc=K$ for some constant $K$ . Letting $u_{1}=a$ , $u_{2}=b$ and $u_{3}=c$ , then Eqs. 15 to 17 can be written as the standard replicator:

\dot{u}_{i}=u_{i}\left(\left(\mathbf{e}_{i}-\mathbf{u}\right)^{T}\mathbf{A}\mathbf{u}\right),

(18)

where:

\mathbf{A}=\begin{bmatrix}-1&-\zeta-1&\zeta-1\\ \zeta-1&-1&-\zeta-1\\ -\zeta-1&\zeta-1&-1\end{bmatrix}

(19)

This is trivially diffeomorphic to the replicator dynamic with payoff matrix:

\tilde{\mathbf{A}}=\zeta\begin{bmatrix}0&-1&1\\ 1&0&-1\\ -1&1&0\end{bmatrix}.

(20)

This is a scaled version of the standard (unbiased) RPS payoff matrix. However, when $\sigma\neq 0$ , then Eqs. 15, 16 and 17 can be written as:

\dot{u}_{i}=u_{i}\left(1+\mathbf{e}_{i}^{T}\mathbf{B}\mathbf{u}\right),

(21)

with:

\mathbf{B}=\begin{bmatrix}-1&-\sigma-\zeta-1&\zeta-1\\ \zeta-1&-1&-\sigma-\zeta-1\\ -\sigma-\zeta-1&\zeta-1&-1\end{bmatrix}.

(22)

Population evolution equations of this form are considered in [41]. When $\sigma=0$ , then $\mathbf{u}^{T}\mathbf{B}\mathbf{u}=\mathbf{u}^{T}\mathbf{A}\mathbf{u}=-1$ , from which we derive Eq. 18. For $\sigma\neq 0$ , this is not the case and consequently the dynamics of Eqs. 15, 16 and 17 are not trapped on the unit 2-simplex as will be the case when we study (un)biased spatial replicator dynamics in finite and infinite populations. To be clear, the authors of [32, 26] make no claim to this effect, however it does create an important distinction between the work in these papers and the work in this paper.

Additionally, as illustrated in Fig. 7, the ODE dynamics that give rise to the PDE are not trapped on the unit simplex (as they are in the replicator) thus suggesting distinct behaviors may be observed.

References

[1] D. Friedman. Evolutionary games in economics. Econometrica, 59:637–666, 1991.
[2] J. W. Weibull. Evolutionary Game Theory. MIT Press, 1997.
[3] J. Hofbauer and K. Sigmund. Evolutionary Games and Population Dynamics. Cambridge University Press, 1998.
[4] J. Hofbauer and K. Sigmund. Evolutionary Game Dynamics. Bulletin of the American Mathematical Society, 40(4):479–519, 2003.
[5] Jun Tanimoto. Fundamentals of evolutionary game theory and its applications. Springer, 2015.
[6] Jun Tanimoto. Evolutionary Games With Sociophysics. Springer, 2019.
[7] Daniele Vilone, Alberto Robledo, and Angel Sánchez. Chaos and unpredictability in evolutionary dynamics in discrete time. Phys. Rev. Lett., 107:038101, Jul 2011.
[8] Andrew J. Black, Arne Traulsen, and Tobias Galla. Mixing times in evolutionary game dynamics. Phys. Rev. Lett., 109:028101, Jul 2012.
[9] Danielle F. P. Toupo and Steven H. Strogatz. Nonlinear dynamics of the rock-paper-scissors game with mutations. Phys. Rev. E, 91(052907), 2015.
[10] K. H. Schlag. Why imitate and if so how? A boundedly rational approach to multi-armed bandits. J. Econ. Theory, 78:130–156, 1997.
[11] D. Fudenberg and D. K. Levine. The theory of learning in games. MIT Press, 1998.
[12] G. Bard Ermentrout, Christopher Griffin, and Andrew Belmonte. Transition matrix model for evolutionary game dynamics. Phys. Rev. E, 93:032138, Mar 2016.
[13] R deForest and A Belmonte. Spatial pattern dynamics due to the fitness gradient flux in evolutionary games. Physical Review E, 87, 2013.
[14] B Kerr, MA Riley, MW Feldman, and Bohannan BJM. Local dispersal promotes biodiversity in a real-life game of rock–paper–scissors. Nature, 418.6894:171–174, 2002.
[15] MA Nowak and RM May. Evolutionary games and spatial chaos. Nature, 359.6398:826–829, 1992.
[16] CP Roca, JA Cuesta, and A Sánchez. Evolutionary game theory: Temporal and spatial effects beyond replicator dynamics. Physics of life reviews, 6.4:208–249, 2009.
[17] G Szabó, A Szolnoki, and R Izsák. Rock-scissors-paper game on regular small-world networks. Journal of physics A: Mathematical and General, 37.7:2599, 2004.
[18] GT Vickers. Spatial patterns and ess’s. Journal of Theoretical Biology, 140(1):129–135, 1989.
[19] GT Vickers. Spatial patterns and travelling waves in population genetics. Journal of Theoretical Biology, 150:329–337, 1991.
[20] R Cressman and GT Vickers. Spatial and density effects in evolutionary game theory. Journal of theoretical biology, 184(4):359–369, 1997.
[21] R. Durrett and S. Levin. The Importance of Being Discrete (and Spatial). Theoretical Population Biology, 46:363–394, 1994.
[22] K. J. Davies, J.E. F. Green, N. G. Bean, B. J. Binder, and J. V. Ross. On the derivation of approximations to cellular automata models and the assumption of independence. Mathematical Biosciences, 254:63–71, 20014.
[23] Alexander van Oudenaarden. Systems Biology, Chapter X: Diffusion (MIT Lecture Notes). http://web.mit.edu/biophysics/sbio/PDFs/L15_notes.pdf, Fall 2009.
[24] Robert M May and Warren J Leonard. Nonlinear aspects of competition between three species. SIAM journal on applied mathematics, 29(2):243–253, 1975.
[25] Mauro Mobilia. Oscillatory dynamics in rock–paper–scissors games with mutations. Journal of Theoretical Biology, 264(1):1–10, 2010.
[26] Claire M Postlethwaite and Alastair M Rucklidge. A trio of heteroclinic bifurcations arising from a model of spatially-extended rock–paper–scissors. Nonlinearity, 32(4):1375, 2019.
[27] Bartosz Szczesny, Mauro Mobilia, and Alastair M Rucklidge. Characterization of spiraling patterns in spatial rock-paper-scissors games. Physical Review E, 90(3):032704, 2014.
[28] Attila Szolnoki, Mauro Mobilia, Luo-Luo Jiang, Bartosz Szczesny, Alastair M Rucklidge, and Matjaž Perc. Cyclic dominance in evolutionary games: a review. Journal of the Royal Society Interface, 11(100):20140735, 2014.
[29] Tobias Reichenbach, Mauro Mobilia, and Erwin Frey. Mobility promotes and jeopardizes biodiversity in rock–paper–scissors games. Nature, 448(7157):1046–1049, 2007.
[30] Tobias Reichenbach, Mauro Mobilia, and Erwin Frey. Self-organization of mobile populations in cyclic competition. Journal of Theoretical Biology, 254(2):368–383, 2008.
[31] Qian He, Mauro Mobilia, and Uwe C Täuber. Spatial rock-paper-scissors models with inhomogeneous reaction rates. Physical Review E, 82(5):051909, 2010.
[32] CM Postlethwaite and AM Rucklidge. Spirals and heteroclinic cycles in a spatially extended rock-paper-scissors model of cyclic dominance. EPL (Europhysics Letters), 117(4):48006, 2017.
[33] Bartosz Szczesny, Mauro Mobilia, and Alastair M Rucklidge. When does cyclic dominance lead to stable spiral waves? EPL (Europhysics Letters), 102(2):28012, 2013.
[34] KM Ariful Kabir and Jun Tanimoto. The role of pairwise nonlinear evolutionary dynamics in the rock–paper–scissors game with noise. Applied Mathematics and Computation, 394:125767, 2021.
[35] E. C. Zeeman. Population dynamics from game theory. In Global Theory of Dynamical Systems, number 819 in Springer Lecture Notes in Mathematics. Springer, 1980.
[36] John Guckenheimer and Philip Holmes. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, volume 42. Springer Science & Business Media, 2013.
[37] Peter D Taylor and Leo B Jonker. Evolutionary stable strategies and game dynamics. Mathematical Biosciences, 40(1-2):145–156, 1978.
[38] J. Hofbauer. On the occurrence of limit cycles in the volterra- lotka equation. Nonlinear Analysis, Theory, Methods and Applications, 5(9):1003–1007, 1981.
[39] Immanuel M Bomze. Lotka-volterra equation and replicator dynamics: a two-dimensional classification. Biological cybernetics, 48(3):201–211, 1983.
[40] Yuzuru Sato, Eizo Akiyama, and J. Doyne Farmer. Chaos in learning a simple two-person game. PNAS, 99(7):4748–4751, 2002.
[41] Elisabeth Paulson and Christopher Griffin. Cooperation can emerge in prisoner’s dilemma from a multi-species predator prey replicator dynamic. Mathematical Biosciences, 278:56 – 62, 2016.