Discovering the effect of nonlocal payoff calculation on the stabilty of ESS: Spatial patterns of HawkDove game in metapopulations
Abstract
The classical idea of evolutionarily stable strategy (ESS) modeling animal behavior does not involve any spatial dependence. We considered a spatial HawkDove game played by animals in a patchy environment with wrap around boundaries. We posit that each site contains the same number of individuals. An evolution equation for analyzing the stability of the ESS is found as the mean dynamics of the classical frequency dependent Moran process coupled via migration and nonlocal payoff calculation in 1D and 2D habitats. The linear stability analysis of the model is performed and conditions to observe spatial patterns are investigated. For the nearest neighbor interactions (including von Neumann and Moore neighborhoods in 2D) we concluded that it is possible to destabilize the ESS of the game and observe pattern formation when the dispersal rate is small enough. We numerically investigate the spatial patterns arising from the replicator equations coupled via nearest neighbor payoff calculation and dispersal.
Keywords: Evolutionary game dynamics, metapopulations, pattern formation.
1 Introduction
Evolutionary game theory is a mathematically accessible way of describing different behavioral traits in a population each of which is associated with a pure strategy of the underlying game. The initial focus of evolutionary game theory was on the concept of evolutionarily stable strategies (ESS) which is used to enhance our understanding of the evolution of animal behavior by [32, 31]. A strategy is defined to be an ESS if a small number of individuals playing a different strategy cannot invade a population playing it. An important question regarding ESS is if such a strategy is attainable. The study by [43] extended the realm of evolutionary game theory to include dynamics. In other words, they introduced the replicator equations relating the ESS concept with the equilibria of these equations [21]. Since then replicator equations are at the core of evolutionary game theory. This classical model describes the evolution of behavioral traits in an infinite population assuming that a given individual is equally likely to interact with any other. Key advances were made by relaxing some of the above mentioned assumptions. In particular, the inclusion of finite populations and spatial structure in evolutionary games have accelerated the progress of the theory.
Stochastic processes have been applied to model evolutionary game dynamics in finite populations [42]. There are a variety of microscopic rules describing the game dynamics in finite populations such as birthdeath update, deathbirth update or pairwise comparison rules [38]. These update rules describes a class of Markov chains, transition probabilities of which are assumed to depend on frequencies of phenotypes and game parameters. For many of these processes, replicator equation is not only a limiting deterministic case [46] but also describes the mean dynamics of the underlying Markov chains.
Considering the fact that natural environments possess a spatial dimension, meaning that individuals have limited mobility and interact with their neighbors, led many scholars to incorporate this important property into the study of evolutionary games. To analyse the effect of spatial structure on evolutionary game dynamics different approaches have been taken in to account: numerical simulations of games on grids (see e.g. [35, 34, 37]) or more generaly on graphs [1, 41]; and analytical studies of replicatordiffusion equations (see, for example, [22, 19, 20, 13]). We would like to note that a similar partial differential equation is found as a hydrodynamical limit of the frequency dependent Moran process [9]. In addition, integrodifferential replicator equations taking nonlocal payoff calculation into consideration were obtained as a mesoscopic limit of the structured individual based models (see e.g. [23, 5]).
Here we let individuals play the HawkDove game, originally developed by [32], to describe certain scenarios in animal conflict modeling a contest over a shareable resource. There are two subtypes or morphs of one species with two strategies, Hawk and Dove. Both subtypes first display aggression. The Hawk escalates into a fight until it either wins or is injured (loses), whereas the Dove runs for safety if faced with major escalation. Unless faced with such escalation, the Dove attempts to share the resource.
Rather than taking the continuous space into account as in replicatordiffusion or integro differential equations, we consider a landscape ecology perspective and subdivide the environment into distinct but identical patches with periodic boundaries each of which contains a population of individividuals. We would like to note that stochastic population models considering a collection of patches on which a number of individuals lives were studied by [45, 2]. In evolutionary game theory literature, a similar approach was taken by [18] to study a finite population of individuals subdivided into demes with the assumption of local interactions meaning that individuals interact only with other members of the same deme.
In our dynamical setting we assume that a patch containing two subtypes of the same species is chosen randomly at each time step. Then an individual is drawn from this site with a probability depending on its fitness and replaces a randomly chosen individual from its natal site with probability or one of the neighboring sites. Such a process can be seen as a spatial coupling of frequency dependent Moran processes [42] via dispersal. Hence identifies the dispersal probability of a newborn. The fitness calculation, on the other hand, has two major components. The first one is the payoff calculation and the other is related to determining the effect of payoffs on fitness.
Using replicatordiffusion equations as a modeling tool leads us to the assumption of local interactions (i.e. local payoff calculation) as in [18]. In the context of the HawkDove game, this is to say that individuals living at a patch do not compete for the resources with the residents of other sites. This very same assumption, on the other hand, was criticized by [7, 12] and relaxed by supposing that reproduction and hence population dynamics take place in a habitat patch whose resources are also used by individuals that live and reproduce in neighboring patches through foraging. This relaxation leads us to the fact that individuals from neighboring sites compete for the common resources and hence the payoff calculation for individuals residing in a patch does not depend solely on the local population configuration but also on the weighted average of the frequencies of the morphs in a certain neighborhood of the site. We would like to remark that there is a vast litrature on individual based models of evolutionary spatial games taking nonlocal payoff calculation in a neighborhood of a spatial location into account (see, for example, pioneering studies by [35, 34]). The second major component of fitness calculation is related to the intensity of selection , a parameter determining the strength of payoffs compared to the baseline fitness. As [39] point out simple as well as illuminating results arise in the limit of weak selection, This is also the case for our model, hence we assume that the effect of payoffs is small when compared to baseline fitness.
The process considered here is a coupled system of Markov chains taking the spatial structure of the environment into account. This coupling between these chains is through dispersal and nonlocal payoff calculation. We find that the mean field dynamics of this coupled system of Markov chains is a coupled system of replicator equations (CRE). Here, we hypothesized that the magnitude of dispersal probability is comperable to that of the small selection parameter This hypothesis is shown to be a requirement to destabilize the ESS of the HawkDove game. In particular, we analyse this limiting deterministic case near the ESS of the underlying HawkDove game and find that small dispersal rate gives rise to spatial pattern formation in 1D and 2D spatial regions with periodic boundaries when the magnitude of the dispersal probability is of order We would like to note that spatial pattern formation is also possible in the strong selection regime. For a brief discussion and illustration of these patterns, see C.
The emergence of spatial inhomogeneity is related to the nonlocal (or quasilocal) interactions that are shown to be a mechanism for spatial pattern formation for a number of ecolgical processes [7, 14, 3, 4, 25, 12, 30, 26, 47]. In these earlier works it was shown that the spatial inhomegeniety is due to the fluctuating density of a species. Whereas, pattern formation in our model describes the fluctuations in the frequencies of the morphs due to the fact that replicator equations are used to model evolution in phenotype space.
The structure of this article is as follows: In Section 2, we give a detailed description of our stochastic model and relate it to deterministic meanfield equations. In Section 3, we perform a linear stability analysis of the model and investigate the conditions for pattern formation. In Section 4, we study the pattern formation numerically for the nearest neighbor interactions and investigate the effects of patch sizes and neighborhood types. Lastly, we discuss and summarize our findings in Section 5.
2 From Coupled Moran Process to CRE
In this section we consider a class of two player symmetric games whose payof matrix is given as follows:

Here we consider a landscape ecology perpective and divide the habitat into identical distinct patches each of which contains individuals. Since our aim is to obtain and analyse the mean dynamics of a stochastic model, we assume that is large. We denote the set of these patches by Each site in the set has a dispersal neighborhood denoted by
Suppose that each individual in the population is either type or At each transition time, a site of origin is chosen randomly, and the following actions take place in order:

An individual from the site is chosen to reproduce according to a frequency dependent probability

With probability the offspring migrates to one of the neighbouring sites in equally likely and replaces a randomly chosen individual in this site.

The offspring replaces a randomly chosen individual from the site of origin with probability
Before proceeding to an introduction of our stochastic model, we discuss how to take a nonlocal payof calculation into account. Suppose that, for any site the frequency of type individuals at time is given by and denote the vector of these frequencies by The payoff calculation is directly related to the foraging range of the species denoted by For the sake of simplicity we suppose that an individual from site is able to play the game with any individual in her site of origin or one of its neighboring sites in and collects her payoff. Hence the foraging neighborhood may be taken as Then the frequency of type individuals in this neighborhood is given by
(1) 
where is used to denote the cardinality of set Hence payoffs of phenotypes and can be calculated as follows:
(2) 
We would like to remark that it is common to exclude selfinteractions in calculating payoffs when finite population models of evolutionary game theory are taken into account. Nevertheless, our assumption that is large enables us to construct our model without this exclusion as done by many scholars (see, for instance, [46, 9, 6]).
Using the selection parameter to describe the degree to which payoffs contribute the fitnesses of these phenotypes leads us to the following fitness functions of type and individuals at site
(3) 
and
(4) 
respectively. Hence, the frequency dependent probability to choose an individual with strategy at site is given by
It is clear that the probability of choosing a type individual is given by
Note that the selection parameter affects the contribution of payoffs to fitnesses of each type. In most studies a small selection parameter () is used to simplify the analysis. This assumption is also biologically reasonable since there might be many other factors effecting the choice of the individual who is going to reproduce. Assuming small selection pressure leads us to the following approximation of the frequency dependent probabilty:
Now we are ready to set up our Markov process. Suppose that transition times of coupled Moran process is given by where denotes its transition probabilities as follows:
It is now possible to construct the transition probabilities of the Markov process. At site an increase in the number of type individuals is possible only if an individual of type is chosen from the site and replaced by a type individual. The probability of chosing a type individual from the origin site is This individual replaces a type individual with probability in site Or an individual of type can be chosen from the neighborhood to reproduce and disperse to site with probability and replaces a type individual from site In a similar way, one can construct the probability of an increase in the number of type individuals. Since we are dealing with a population of populations on a grid, we assume that each site has the same number of neighboring sites i.e. for any In the following section, we specify a few neighborhoods to study the stability of the mixed strategy ESS of the HawkDove game.
Here we hypothesize that the magnitude of the dispersal probability is comperable with that of the selection parameter i.e. for some positive constant This weak dispersal hypothesis is required to simplify the equations governing the population dynamics and to destabilize the space homogenous ESS of the game in the weak selection regime. The former simplification is due to the fact that For a justification of the latter, see Remark in Section 3. Destabilization of the ESS of the HawkDove game is also possible (see C) when the magnitude of dispersal probability is comperable to that of the selection parameter which is not assumed to be small. Consequently, we have the following transition probabilities in the weak selection regime:
and
Let be a solution to coupled replicator equations (CRE) that is given as follows:
(5) 
where the discrete dispersal operator is defined by the following equality:
for any
To connect the above given discrete time Markov chain with a continuous time differential equation we need to consider a piecewise affine interpolation of the process at any time as follows:
for any Connection between the above defined interpolated Markov chain and CRE (5) is given by the following lemma.
Lemma 1.
For any finite time there is a scalar such that
for any and all
The above given result can be proved following the proof of [6, lemma 1] which only requires replacing the vector field used in their paper with the Lipchitz continuous vector field associated with equation (5). The result implies that the trajectories of the Markov chain and equation (5) deviate from each other with a probability that exponentially approaches to zero as increases. In particular, the above given bound can be used to obtain a strong law of large numbers result for any finite time as follows:
which mimics the result obtained by [27] for continuous time density dependent Markov chains.
3 Stability of Mixed Strategy ESS of the HawkDove Game
Recall the classic HawkDove game involving two discrete strategies. Individuals playing a Hawk strategy always fight to injure or kill their opponent. Whereas individuals employing the Dove strategy always display and never escalate the contest to serious fighting. If two individuals meet and both adopt the Hawk strategy, at least one will be seriously injured in the contest. Similarly, if two players both adopt the Dove strategy, there is some cost to continued displaying. When one player adopts a Hawk strategy and the other plays Dove, the Hawk wins the contested resource.
The payoff matrix of a HawkDove game satisfies and To simplify the notation we introduce the following parameters
(6) 
In terms of these two parameters above the given condition can be repeated as and the mixed strategy ESS of the game is given by Hence, we denote the spatially homogenous solution of the CRE in vector form
We study the stabiliy of mixed strategy ESS of the HawkDove game under CRE. In Section 2, we did not specify the neighborhoods or the spatial dimension of the habitat. Here we first describe the simple nearest neighbor dispersal and interaction in 1D and 2D spatial habitats used in both ecological studies by [12, 47, 26, 30, 8], and game theoretical studies initiated by [34, 35]. Then we give a general form of the mean field equation and find the conditions under which the ESS can be destabilized.
3.1 Calculating Eigenvalues in 1D Habitat
Consider a one dimensional periodic habitat i.e. which forms a ring. Such a fragmented habitat may be reasonable for a species living around a lake or a pond. The nearest neighborhood interaction (or payoff calculation) makes sense if the foraging behavior of animals is taken into account. Thus, we consider the following dispersal and interaction neighborhoods:
(7) 
The former neighborhood automatically identifies a one dimensional discrete dispersal operator:
(8) 
which is widely used as discretization of the diffusion term in reactiondiffusion equations (see e.g. [8, 30]).
In order to describe the use of the latter neighborhood we introduce the discrete convolution operator of two vectors and and denote it by The definition of dicrete convolution is given as follows:
(9) 
where Hence, the average frequency of type individuals around any site can be written as a discrete convolution operator
(10) 
where is the interaction kernel whose nonzero elements are given as Hence (5) with an appropriate time scaling can be written as
(11) 
where parameters and are as defined in (6).
We linearize (11) by using the first order expansion where is a vector independent of time denoting a spatial perturbation term. Hence this formula can be considered as a separation of variables formula for (11). Substituting this ansatz to the equation gives the following first order relation:
(12) 
To obtain an explicit expression for the eigenvalues of the above given linear equation we need to use discrete Fourier series (DFS). Following [40, 29], DFS of a vector is denoted by and given by
(13) 
for any Note that is used to denote the imaginary unit number.
One can easily see that the DFS of the convolution of two vectors and is given by the Hadamart product of DFS of these vectors i.e. we have:
(14) 
Note that the operator can be written as the sum of a convolution operator and an identity operator (see A.1). Hence taking the DFS of both sides of equation (12) one obtains explicit equalities for the eigenvalues as follows:
(15) 
where for any The explicit calculations are given in A.1. In particular, we obtained the above equality by using formulas (30) and (31).
3.2 Calculating Eigenvalues in 2D Habitat
Extending the 1D neighborhood to 2D may lead us to different neighborhoods. We consider two major neighborhoods used both in game theory literature for cellular automata simulations (see e.g. [33]) and in metapopulation studies (see e.g. [12]). First consider the von Neumann neighborhood of site which can be defined as folows:
(16) 
Another extension of the one dimensional nearest neighbour interaction range is so called Moore neighborhood defined as follows:
(17) 
Hence dispersal neighborhoods are given as for Note that these two neighborhood types are the simplest generalizations of 1D case. Projecting both neighborhoods to 1D leads us to the one dimensional neighborhood
Now we are ready to define the discrete dispersal operators for von Neumann and Moore neighborhoods as follows:
(18) 
.
To write down the CRE in 2D habitat, we need to define the 2D discrete convolution operator of two vectors (or matrices) and which is given as follows:
(19) 
Therefore, (5) in 2D habitat with an appropriate time scaling can be written as follows:
(20) 
for Here and are used to denote discrete uniform probabiliy distributions on the neighborhoods and of any site respectively. The explicit expressions for these discrete probability distributions are given in A.2.
Following [40, 29], discrete Fourier series (DFS) of a matix is denoted by and given by
(21) 
for any Note that the discrete convolution theorem also holds in 2D and states that
(22) 
Similarly using the first order expansion and taking the DFS of resulting linear equation, we obtain the explicit expressions for the eigenvalues.
3.3 Generalization of the CRE and instability of the ESS
In the above subsections we only considered the nearest neighbor interactions and derived the explicit expressions for the eigenvalues of the linearization of CRE. It is possible to write this equation down in a more general form and identify the conditions on the interaction kernel and dispersal rate to destabilize the ESS. Recall the discussion in Section 3.1 stating that the discrete dispersal operator can be written as the sum of a convolution and the identity operator i.e. for a discrete probability distribution This is also true in 2 spatial dimensions (see A.2). Hence we can write the CRE as
(25) 
where and are some discrete probability kernels. Linearizing this equation around as done in previous subsections and taking the DFS leads us to the following eigenvalues:
(26) 
Since both and are discrete probability distributions their DFS (which is fully determined by the characteristic functions of the kernels) takes values inbetween and Hence the real part of the term must be nonpositive. This implies that the dispersal has a stabilizing effect on the ESS. Destabilization of the ESS is only possible if the real part of the term is negative for some and the dispersal rate is small enough. We state this result in the following theorem.
Theorem 1.
The ESS of the game is unstable under the equation (25) if the real part of the DFS of the interaction kernel takes negative values and the dispersal rate is small enough.
Nonlocal interactions for a given number of groups result in the overlapping of foraging areas that can be thaught as a mechanism altering the frequencies of morphs forming clusters. The above result, therefore, may imply the existence of clusters altering around the ESS for some foraging range and frequencies determined by the nonlocal interaction kernel
In Section 1 we mentioned the studies of evolutionary games via replicatordiffusion equations. With the assumption that all diffusion rates are equal; replicatordiffusion equations describe a monotone operator and one cannot seek a mechanism destabilizing the ESS of the game. To study diffusion driven intabilities [48, 10] introduced a special form of the population regulation to allow for different diffusion rates. In particular, they used a nonlinear term providing a local regulation of the frequencies of the morphs to keep solutions to replicatordiffusion equations in an extension of the game simplex. The following result related to generalized CRE ensures that the frequency of each type stays in between 0 and 1 at each patch i.e. it stays in the extended simplex where is the number of sites into which habitat is divided.
Theorem 2.
For any initial datum solution to CRE (5) stays in the extended simplex.
The proof of this result is given in B. Note also that the above theorem is valid not only for the HawkDove game but all symmetric two player games.
Now we verify that the hyphotheses of the Theorem 1 holds for the 1D and 2D cases examined above. To shorten the analysis we write the eigenvalues in a common form. Recall that the cardinality of the neighborhood is denoted by for Hence all three eigenvalues can be written as follows:
(27) 
for Note that
(28) 
This implies that for all Hence we can conclude that the dispersal has a stabilizing effect on the ESS for any
Remark 3.
Since all three nearest neighbour interaction kernels are related to the uniform distribution, the DFS of each kernel take negative values. The DFS of the kernel for example, takes negative values since for some values of near Similarly, it can be easily seen that the DFS of the interaction kernels take negative values for von Neumann and Moore neighborhoods.
Suppose that then the ESS of the game is unstable for nearest neighbour interactions in 1D or 2D habitats. By continuity it is also unstable for small enough. On the other hand it is stable for large enough Hence, we conclude that there exists a critical diffusion rate such that the ESS is stable for and unstable for Such a critical diffusion rate can be calculated by solving the following root finding problem:
(29) 
for The argument of this maximization problem is denoted by and called the most unstable wavenumber.
4 Numerical Simulations
Throughout this section, we consider the following payoff matrix

descrribing a HawkDove game with and
First we consider the equation (11) in a 1D habitat forming a ring. Our model with nearest neighborhood interaction and dispersal gives rise to two types of patterns in phenotype space. The emergence of these two types of patterns are related to the number of patches contained in the habitat In particular the number of patches determines the number of competing eigenvalues and the emergence of different spatial patterns.
Figure 1 panels (a) and (b) illustrates these patterns for and For there exists only one most unstable wave number giving rise to patterns shown in Figure 1 Panel (a). More generally, when the number of patches are even we have only one the most unstable wavenumber which is given by on which the DFS of the interaction kernel takes its smallest value. We refer to these patterns as regular stationarywave type patterns. On the other hand, when the patch number is odd there are two competing wavenumbers causing the emergence of a different type of stationary patterns. An example of such patterns for is shown in Figure 1 Panel (b). We refer to these patterns as degenerate stationary patterns.
For our numerical calculations we used (29) and the eigenvalues (15) to calculate the critical dispersal rates for the parameters given in the above payoff matrix. Having this value enables us to calculate the most unstable wave number(s) giving rise to patterns shown in Figure 1. The critical dispersal rates and associated wave numbers for and are given as follows.

for for which the most unstable wave number is associated with the eigenvalue

for where we get two competing eigenvalues namely and
Next we consider the equation (20) in a 2D habitat with von Neumann neighborhood. Again the types of patterns arising from this equation are related to the number of patches contained by the habitat. For we observed regular stationary patterns (see Figure 2 Panel (a)) for for which we have only one unstable wavenumber. For and we observed the emergence of patterns shown in Figure 2 Panel (b). We would like to call attention to a similarity between these patterns and the ones shown in Figure 1. For a fixed one can easily see that these patterns are similar to the one obtained in Figure 1 Panel (b). On the other hand, the obtained patterns are similar to 1D patterns shown in 1 Panel (a) for fixed Hence we call these patterns semidegenerate stationary patterns. Finally, considering the case leads us to the patterns shown in Figure 2 Panel (c). These patterns can be thought as a 2D degenetarte stationary patterns since fixing either or results in degenerate 1D degenerate patterns shown in Figure 1 Panel (b).
For our numerical calculations we again used (29) and the eigenvalues (23) to calculate the critical dispersal rates for the game parameters given in the above payoff matrix. The critical dispersal rates and associated wave numbers for different values of and are given as follows.

for Here the most unstable wavenumber is associated with the eigenvalue

for and Here we get two competing eigenvalues namely and

for Here we get four competing eigenvalues namely and
We remark that the values of and are chosen only for illustration purposes. The factor determining the number of unstable wavenumbers is related to the evenness or oddness of these values as in the 1D case.
Finally we consider the equation (5) in 2D habitat with Moore neighborhood. Again the patterns arising from this equation is related to the number of patches contained by the environment. We observed patterns shown in Figure 3 Panel (a) for for which there exist two competing unstable wavenumbers. This patters can be considered as a composition of two different 1D regular patterns. If the patch numbers are taken as we observe strip patterns shown in Figure 3 Panel (b). For a fixed one can easily see that these patterns are similar to the regular patterns obtained in Figure 1 Panel (a). When we observe the patterns shown in Figure 2 Panel (c). These patterns can be considered as a composition of two different 1D degenetarte patterns.
For our numerical calculations we again used (29) and the eigenvalues (24) to calculate the critical dispersal rates for the parameters given by the above payoff matrix. The critical dispersal rates and associated wave numbers for different values of and are given as follows.

for Here we get two competing eigenvalues namely and

for and Here the most unstable wavenumber is associated with the eigenvalue

for Here we get four competing eigenvalues namely and
It is easy to observe that Moore neighborhood gives rise to more complex patterns in general. The reason behind this observation is that the eigenvalues associated with Moore neighborhood identifies distinct wavenumbers. Whereas all of the most unstable wavenumbers asociated with the von Neumann neighborhood are near
5 Discussion
In this paper, we have obtained the coupled replicator equations (CRE) for two player symmetric games as mean dynamics of a frequency dependent Moran process coupled via spatial dispersal and nonlocal payoff calculation in a patchy environment. In particular, we studied the stability of the mixed strategy ESS of the HawkDove game in this patchy environment. We have shown that the spatial dispersal operator has a stabilizing effect on the spatially homogenous solution of CRE (i.e. mixed strategy ESS of the game) as in other ecological models (see for example [16, 15, 11]) while nonlocal payoff calculation may destablize the ESS. In this paper we only explicitly calculate the Fourier transforms for kernels modeling nearest neighborhood interactions and dispersal in 1D and 2D periodic habitats to keep the technicality to a minimum. However, we also explicitely stated the form that the CDE takes when one needs to consider more complex dispersal and interaction structures to incorporate movement to further patches, or to include unequal weighting. As a result, we found that the ESS of the game is unstable if the (real part of) interaction kernel takes negative values and the dispersal rate is sufficiently small.
We also numerically studied the patterns arising from CRE with nearest neighbour symmetric dispersal and interaction kernels. We observed that the evenness or oddness of the number of sites affect the shape of the spatial patterns. If the number of patches are even and the habitat is one dimensional, the arising patterns are similar to the regular stationary wave type patterns arising from reaction diffusion equations. Whereas, If the number of patches are odd we observed the formation of degenerate stationary patterns (see Figure 1 Panel (b)) For von Neumann neighborhood, we see that the patterns in 2D habitat as shown in Figure 2 are just 2D generalizations of the ones observed in Figure 1. On the other hand, our numerical study indicates that Moore neighbourhood gives rise to the emergence of more complex patterns when compared to the ones arising from CRE coupled via von Neumann neighborhood (see Figure 3).
The effect of spatial structure in discrete individual based models was first mentioned by [35, 36] for the HawkDove game. Subsequently, it was discussed in more detail by [24] whose numerical simulations indicated that, the consideration of the spatial structure enables doves to become spontaneously organised into clusters so that they can outweigh losses against hawks. The patterns obtained here can be interpreted as a mechanism for spontaneous cluster formation as done by [49]. Hence our results support the numerical findings of [24]. However, [17] considered the HawkDove game on a lattice and remarked that spatial extension generally favours the hawk strategy when compared to the replicator equations. Our model, however, suggests that a strategy is favored when compared to the mean field dynamics only if the ESS corresponding to this strategy is less than (see Figure 4).
It was revealed in the literature that pattern formation can be observed when modified replicatordiffusion models are used [48]. In addition, there are further studies exploring spatially extended public goods games and it was shown that reactiondiffusion equations modeling the public good phenomenon exhibits not only pattern formation but also chaotic dynamics [49]. In both continuous space models it was assumed that the diffusion rates of two morphs are different. Thus instabilities observed from these models are essentially diffusiondriven instabilities. Contrary to these studies our model takes nonlocal interactios or payoff calculation into account. Such interactions can be considered as another mechanism destabilizing the ESS.
Finally, there are several noteworthy issues that should be further explored or extended. In this paper, we focused on twostrategy games with an interior ESS. The derived results may be extended to strategy games with , such as the RockScissorsPaper game. In addition, following [4], a nonlinear analysis of the CRE can be performed to study the phase transitions near stability boundaries. We also would like to note that one can extend the analysis performed in Section 3 from the spatial lattices to more general networks. Such an extension, however, requires a welldefined Fourierlike transformation(s) on the generalized network structure. An apparent example of such graphs is Cayley graphs on which Fourier transform is defined [44]. For wavelet transformation on more general graphs, see [28, Section 8.2].
Appendix A Explicit Calculation of Eigenvalues
a.1 1D Case
Note that the discrete dispersal operator can be written as where is the discrete probability kernel on with and is the identitiy operator. By definition of one dimensional DFS one can calculate the DFS of the kernel as follows:
(30)  
Similarly one can calculate the DFS of the interaction kernel with as
(31) 
a.2 2D Cases
Suppose that if the individuals interact in a von Neumann neighborhood, the discrete dispersal operator can be written as where is the discrete probability kernel on with and is the identitiy operator. By definition of two dimensional DFS one can calculate the DFS of the kernel as follows:
(32)  
The corresponding interaction kernel is given by and its DFS is given as follows:
(33) 
Assuming Moore neighborhood, we get the following equalities:
(34) 
and
(35) 
Appendix B Proof of Theorem 2
Let be two functions inthe space of bounded uniformly continuous functions on equipped with the norm Suppose that for all and let be a solution to equation
(36) 
with the initial data Here the discrete probability kernel is assumed to be a step function with support Since is a discrete probability distribution, for any Notice that if for some and some The trajectory of is therefore always beyond 0. Similarly, if for some and some Hence, for any
Now we will show that mapping is a contraction. Let and be two functions satisfying the above conditions and suppose that and denote the corresponding solutions to equation (36) with the same initial data. Let and , so that
(37) 
where and . Integration up to time yields
which implies
Let for , and take the supremum over and , we get
For sufficiently small , . Hence, for sufficiently small , is a contraction preserving the closed subset and therefore admits a unique fixed point. Since is independent of the size of the data, this argument can be iterated from to any desired time interval.
Appendix C The MeanField Equation in Strong Selection Regime
Here we consider the meanfield equation of the frequency dependent Moran process when the selection parameter is large. In particular, we suppose that yielding the original Markov chain studied by [42]. In this case, the fitness of each morph equals to its frequency dependent payoff (2). This implies the following equality:
Using this probability, one can construct the Markov process to yield its mean dynamics. The meanfield equation of this process is given by the following equation:
where and Here we would like to note that the spatially homogenous solution loses its stability for sufficiently small dispersal probability Figure 5 illustrates degenerate stationary patterns for Q = 21.
References
 Allen, B., Nowak, M. A., 2014. Games on graphs. EMS surveys in mathematical sciences 1 (1), 113–151.
 Arrigoni, F., Pugliese, A., 2002. Limits of a multipatch sis epidemic model. Journal of mathematical biology 45 (5), 419–440.
 Aydogmus, O., 2015. Patterns and transitions to instability in an intraspecific competition model with nonlocal diffusion and interaction. Math. Mod. Nat. Phen. 10 (6), 17–29.
 Aydogmus, O., 2017. Phase transitions in a logistic metapopulation model with nonlocal interactions. Bulletin of mathematical biology.
 Aydogmus, O., Zhou, W., Kang, Y., 2017. On the preservation of cooperation in twostrategy games with nonlocal interactions. Mathematical biosciences 285, 25–42.
 Benaïm, M., Weibull, J. W., 2003. Deterministic approximation of stochastic evolution in games. Econometrica 71 (3), 873–903.
 Britton, N., 1989. Aggregation and the competitive exclusion principle. J. Theor. Biol. 136 (1), 57–66.
 Cantrell, R. S., Cosner, C., Fagan, W. F., 2012. The implications of model formulation when transitioning from spatial to landscape ecology. Math. Biosci. and Eng. 9 (1), 27–60.
 Chalub, F. A., Souza, M. O., 2009. From discrete to continuous evolution models: a unifying approach to driftdiffusion and replicator dynamics. Theoretical Population Biology 76 (4), 268–277.
 Cressman, R., Vickers, G., 1997. Spatial and density effects in evolutionary game theory. Journal of Theoretical Biology 184 (4), 359–369.
 Doebeli, M., 1995. Dispersal and dynamics. Theor. Pop. Biol. 47 (1), 82–106.
 Doebeli, M., Killingback, T., 2003. Metapopulation dynamics with quasilocal competition. Theor. Pop. Biol. 64 (4), 397–416.
 Ferriere, R., Michod, R. E., et al., 2000. Wave patterns in spatial games and the evolution of cooperation. The Geometry of Ecological Interactions: Simplifying Spatial Complexity, 318–339.
 Genieys, S., Volpert, V., Auger, P., 2006. Pattern and waves for a model in population dynamics with nonlocal consumption of resources. Math. Mod. Nat. Phen. 1 (01), 63–80.
 Gyllenberg, M., Söderbacka, G., Ericsson, S., 1993. Does migration stabilize local population dynamics? analysis of a discrete metapopulation model. Math. Biosci. 118 (1), 25–49.
 Hastings, A., 1993. Complex interactions between dispersal and dynamics: lessons from coupled logistic equations. Ecology 74 (5), 1362–1372.
 Hauert, C., 2002. Effects of space in 2 2 games. International Journal of Bifurcation and Chaos 12 (07), 1531–1548.
 Hauert, C., Imhof, L. A., 2012. Evolutionary games in deme structured, finite populations. Journal of theoretical biology 299, 106–112.
 Hofbauer, J., 1998. Equilibrium selection via travelling waves. In: Game theory, experience, rationality. Springer, pp. 245–259.
 Hofbauer, J., Hutson, V., Vickers, G. T., 1997. Travelling waves for games in economics and biology. Nonlinear Analysis: Theory, Methods & Applications 30 (2), 1235–1244.
 Hofbauer, J., Sigmund, K., 1998. Evolutionary games and population dynamics. Cambridge university press.
 Hutson, V., Vickers, G. T., 1992. Travelling waves and dominance of ess’s. Journal of Mathematical Biology 30 (5), 457–471.
 Hwang, S.H., Katsoulakis, M., ReyBellet, L., 2013. Deterministic equations for stochastic spatial evolutionary games. Theoretical Economics 8 (3), 829–874.
 Killingback, T., Doebeli, M., 1996. Spatial evolutionary game theory: Hawks and doves revisited. Proceedings of the Royal Society of London B: Biological Sciences 263 (1374), 1135–1144.
 Killingback, T., Loftus, G., Sundaram, B., 2013. Competitively coupled maps and spatial pattern formation. Physical Review E 87 (2), 022902.
 Kisdi, É., Utz, M., 2005. Does quasilocal competition lead to pattern formation in metapopulations? an explicit resource competition model. Theor. Pop. Biol. 68 (2), 133–145.
 Kurtz, T. G., 1978. Strong approximation theorems for density dependent markov chains. Stochastic Processes and their Applications 6 (3), 223–240.
 Lézoray, O., Grady, L., 2012. Image processing and analysis with graphs: theory and practice. CRC Press.
 Mandal, M., Asif, A., 2007. Continuous and discrete time signals and systems. Cambridge University Press.
 Maruvka, Y. E., Shnerb, N. M., 2006. Nonlocal competition and logistic growth: patterns, defects, and fronts. Phys. Rev. E 73 (1), 011903.
 Maynard Smith, J., 1974. The theory of games and the evolution of animal conflicts. Journal of theoretical biology 47 (1), 209–221.
 Maynard Smith, J., Price, G., 1973. The logic of animal conflict. Nature 246, 15.
 Nowak, M. A., 2006. Evolutionary dynamics. Harvard University Press.
 Nowak, M. A., Bonhoeffer, S., May, R. M., 1994. More spatial games. International Journal of Bifurcation and Chaos 4 (01), 33–56.
 Nowak, M. A., May, R. M., 1992. Evolutionary games and spatial chaos. Nature 359 (6398), 826–829.
 Nowak, M. A., May, R. M., 1993. The spatial dilemmas of evolution. International Journal of bifurcation and chaos 3 (01), 35–78.
 Nowak, M. A., Sigmund, K., 2000. Games on grids. The Geometry of Ecological Interactions: Simplifying Spatial Complexity, 135–150.
 Ohtsuki, H., Hauert, C., Lieberman, E., Nowak, M. A., 2006. A simple rule for the evolution of cooperation on graphs and social networks. Nature 441 (7092), 502–505.
 Ohtsuki, H., Nowak, M. A., 2006. Evolutionary games on cycles. Proceedings of the Royal Society of London B: Biological Sciences 273 (1598), 2249–2256.
 Smith, J., 2007. Mathematics of the discrete Fourier transform (DFT): with audio applicaitons. W3K Publishing.
 Szabó, G., Fath, G., 2007. Evolutionary games on graphs. Physics reports 446 (4), 97–216.
 Taylor, C., Fudenberg, D., Sasaki, A., Nowak, M. A., 2004. Evolutionary game dynamics in finite populations. Bulletin of mathematical biology 66 (6), 1621–1644.
 Taylor, P. D., Jonker, L. B., 1978. Evolutionary stable strategies and game dynamics. Mathematical biosciences 40 (1), 145–156.
 Terras, A., 1999. Fourier analysis on finite groups and applications. Vol. 43. Cambridge University Press.
 Tilman, D., Kareiva, P. M., 1997. Spatial ecology: the role of space in population dynamics and interspecific interactions. Vol. 30. Princeton University Press.
 Traulsen, A., Claussen, J. C., Hauert, C., 2005. Coevolutionary dynamics: from finite to infinite populations. Physical review letters 95 (23), 238701.
 Utz, M., Kisdi, É., Doebeli, M., 2007. Quasilocal competition in stagestructured metapopulations: a new mechanism of pattern formation. Bull. Math. Biol. 69 (5), 1649–1672.
 Vickers, G., 1989. Spatial patterns and ess’s. Journal of theoretical biology 140 (1), 129–135.
 Wakano, J. Y., Hauert, C., 2011. Pattern formation and chaos in spatial ecological public goodsgames. Journal of theoretical biology 268 (1), 30–38.