Mean-field interactions in evolutionary spatial games

We introduce a mean-field term to an evolutionary spatial game model. Namely, we consider the game of Nowak and May, based on the Prisoner's dilemma, and augment the game rules by a self-consistent mean-field term. This way, an agent operates based on local information from its neighbors and nonlocal information via the mean-field coupling. We simulate the model and construct the steady-state phase diagram, which shows significant new features due to the mean-field term: while for the game of Nowak and May, steady states are characterized by a constant mean density of cooperators, the mean-field game contains steady states with a continuous dependence of the density on the payoff parameter. Moreover, the mean-field term changes the nature of transitions from discontinuous jumps in the steady-state density to jumps in the first derivative. The main effects are observed for stationary steady states, which are parametrically close to chaotic states: the mean-field coupling drives such stationary states into spatial chaos. Our approach can be readily generalized to a broad class of spatial evolutionary games with deterministic and stochastic decision rules.


Introduction-
The mean-field approximation, initially devised in 1907 in the context of (ferro)magnetic properties of materials [1], is since playing a truly central role in a variety of branches of physics, ranging from superconductivity [2] to ferromagnetism of disordered metals [3]; from ultracold quantum gases [4] to spin glasses [5], to name just a few. The main idea-that a behavior of a collection of a macroscopic number of strongly interacting components can be reduced to a single-particle problem in a self-consistent external potential-proves fruitful for multiple research fields outside of traditional physics: for instance, describing spatially extended networks [6], droplet traffic [7], traffic flow [8], and social dilemmas [9].
In this Letter, we demonstrate how mean-field ideas can be applied to a system that does not have an explicit representation in the language of statistical mechanics. Instead, it is a deterministic dynamical system with a steady-state, i.e., a long-lived state with a well-defined "thermodynamic limit" of the infinite system size. We start with the completely deterministic evolutionary game model of Nowak and May [14] and add a novel ingredient, the coupling to the instantaneous mean density of cooperators. The modification does not change the deterministic nature of the process in contrast with the previous research (see, f.e., [31][32][33][34][35]) in which the stochastic nature of modifications obviously opens the possibility to using methods of statistical physics. Once again, our model is deterministic (specifically, there is no temperature, no fluctuations, and no random forces). The introduction of the mean-field interaction leads to a nontrivial diagram of steady states and opens new directions for possible applications of structured multi-agent systems. Our approach can be generalized to any deterministic or stochastic spatial evolutionary game.
Studying the dynamics of ensembles of agents is a common theme in biology, economics, and social sciences. Collective behavior emerges through repeated pairwise interactions between individuals, where each agent acts to maximize its fitness. The evolutionary game theory approach model elementary interactions as gametheoretic contests between agents having varying strategies [10]. For well-mixed populations-in physics language, this corresponds to a mean-field approximationthe macroscopic description proceeds via the master equation for the time evolution of populations of strategies, the so-called replicator equation [11,12].
One limitation of the replicator equation approach is that it does not allow for structured populations, where the interaction range is restricted to a local neighborhood. Such populations are known to exhibit long-time behavior beyond standard replicator dynamics [13]. Various local arrangements have been considered: following the pioneering work [14], noisy dynamics were studied on regular, decorated, and frustrated two-dimensional lattices [13,15,16] and coevolving random networks [17]. Deterministic imitation dynamics have been investigated on regular square [18,19] and triangular lattices [20], on a simple cubic lattice 3D [21], and diluted 2D lattices [19,22].
In the Letter, we construct an evolutionary game model which features both kinds of effects: purely local interactions with nearest neighbors, and a selfconsistent, mean-field-type coupling to the order parameter. We start with the classic game of Nowak and May (NM) [14,26], where the order parameter can be chosen as the mean density of cooperators f c in the steadystate regime of evolution. The mean-field term is the coupling to the instantaneous density of cooperators, f c (t), which fluctuates due to the game dynamics and has a well-defined mean value in the steady-state regime. The addition of the MF term to the decision-making rule is quite realistic. It can be interpreted as the influence of the "averaged strategy of society" or a mass-media influence that reflects and drives society's opinions. In the arXiv:2107.11088v2 [cond-mat.stat-mech] 24 Sep 2021 limit of vanishing mean-field (MF) term, the phase diagram reproduces known results.
We note that our approach is very different from what is known as mean-field games in mathematical economics [23]: the latter is a class of stochastic differential games (which, in turn, can be mapped onto the nonlinear Schroedinger equation [24]). In our work, we only consider discrete deterministic games based on the Prisoner's dilemma.
The Prisoner's dilemma.-Consider two agents, α and β, which interact via the rules of the Prisoner's dilemma game. An agent is characterized by a strategy, which takes one of two values: cooperate, C, or defect, D. In an "interaction", agents receive payoffs which depend on their strategies [25]. Encoding strategies of an agent by length-two vectors, s, such that C=(1, 0) T and D=(0, 1) T , the payoff of an agent α is P α = s T α H s β , and the payoff of the agent β is P β = s T β H s α , where H is the payoff matrix, and s α and s β are the states of agents α and β, respectively. Following Ref. [14], we take Essentially, Eq. (1) means that in a C-C interaction, both agents receive a unit payoff (we set it to unity to fix the payoff scale), and in a D-C interaction, the D receives a payoff of b. We set all other payoffs to zero so that a single parameter, b, governs the game. In general, one may vary the two additional parameters, which correspond to the payoffs in the D-D and C-D interactions. We keep both these parameters equal zero for simplicity and to allow a direct comparison with Ref. [14].
The rules of the Nowak-May spatial game [14,26]. -Consider a population of L 2 players arranged at the vertices of a regular L×L square lattice with periodic boundary conditions. The game is played at discrete time steps, t=0, 1, 2, . . . , T . The transition of the system from time t to the next time step t+1 consists of the update of strategies of all players in parallel (we call this a game round). At time step t, an agent located at the lattice site x "interacts" with its neighbors and receives a total payoff, where s x ≡ s x (t) is the strategy of an agent x at a time step t, and the summation runs over the eight neighbors of the agent x (similar to the chess king moves) [36]. After all, payoffs are calculated, each player changes its state, s x (t)→ s x (t+1), to accept the strategy of their neighbor with the largest payoff, This concludes the game round, and the game repeats with the new strategies of players.
To fully specify the time dynamics of the game, we need to define the initial state at t = 0. We follow Ref. [14] and start the evolution with a random initial state where each player is C with a probability f 0 and a D with a probability 1 − f 0 . In the numerical simulations, we average over several realizations of initial conditions at fixed f 0 .
Note that the initial conditions provide the only source of randomness in the game: the dynamics is fully deterministic given the initial state of all strategies at t=0. The strategy of a player x at the next time step, s x (t+1), depends on the current state of itself, s x (t), its eight neighbors, and their neighbors-i.e., on states of 25 players in total. In the cellular automata language, the spatial evolutionary game of Nowak-May (NM game) is a cellular automaton with the single-player's transition matrix containing 2 25 rules.
The steady-state of the spatial evolutionary game.-The game starts with some initial distribution of strategies at t=0 and proceeds repeatedly. After a certain transition time, the system reaches a steady-state, with the nature of the steady-state strongly dependent on the value of the game parameter b [14]. Following Ref. [14], we characterize the instantaneous state of the game by the density of cooperators, f c (t), and the steady-state is characterized by the mean density of cooperators f c : The summation in Eq. (5) runs over the states of all agents at time t, and δ ( s x (t)=C) is an indicator variable taking the value of 1 if s x (t)=C and zero otherwise. In other words, Eq. (5) is nothing but the number of cooperators divided by the total number of agents. In Eq. (4), T is the total number of game rounds, and τ is the burn-in time such that | f c −f c (t)| f c for t>τ . The mean-field modification.-In this Letter, we include the mean-field term as follows: instead of Eq (2) we consider payoffs of the form where In other words, the interaction of a player with its neighbors is identical to the the NM game, Eq. (2). Additionally, we introduce coupling to the instantaneous density of cooperators, f c (t). This way, transition rules become non-local in space due to the dependence on the density of cooperators at the current time step f c (t). Here, we stress that f c (t) is the instantaneous density so that the game is memory-less. In Eqs. (6)- (7), λ is the relative strength of the MF coupling. The limit λ→0 is simply the NM game, Eq. (2). In this work, we only consider λ=1, so that b is the only parameter that governs the game dynamics.
The transition from time step t to t+1 is identical to the NM game: After all payoffs are calculated, each agent changes its strategy according to Eq. (3), and the game repeats with the new strategies of players.
Dynamical regimes.-The discrete structure of the payoffs in the NM game (1)-(2) leads to a very specific dependence of the game dynamics on the payoff parameter b. Comparing the payoffs of C and D in the neighborhood of an agent, one finds for the NM game [18,26] that the dynamics of the game changes at the values of b=m/n with 1 m, n 8. Here m and n are integers and are related to the numbers of C in a local neighborhood. The steady-state density of cooperators, f c takes discontinuities at these special values of b, as shown in Fig. 1(a) with the open circles.
For the mean-field game (6)-(7) (MF game), a similar comparison of the payoffs of neighboring agents shows that an agent changes state depending on the ratio of m+f c (t) and n+f c (t) -which self-consistently depends on the density of cooperators at time t. Since in the steady state we expect f c (t)≈ f c , we conjecture that the steady state density is related to the payoff parameter via where m and n are integers and 1 m, n 8. Note that Eq. (8) cannot be directly interpreted as giving the boundaries of dynamical regimes. To clarify the status of Eq. (8) we compare its predictions to direct numerical simulations. Numerical simulations.-We simulate the dynamics of the game on the square lattice with linear size L=400. We use periodic boundary conditions to minimize finitesize effects. Additional simulations with system size up to L=800 demonstrate no significant finite-size effects. Therefore we believe our simulations correspond to the thermodynamic limit of the model [37]. For all values of b, we simulate 100 runs with random initial distributions of C and D using the initial density of cooperator states f 0 . Each configuration converges to a particular steady-state, and we compute averages for both time fluctuations in the steady-state and realizations of initial conditions. The steady-state behavior is largely insensitive to the precise value of f 0 , and we fix f 0 =0.9 in most simulation runs [38], which is the value used in Refs [18,26]. The number of game rounds is typically a few thousand for the burn-in time towards the steadystate, τ , cf Eq. (4), which is larger than the typical relax-ation time for the initial state. After burn-in, we collect statistics for a few thousand further game rounds, T −τ .  Another remarkable feature of the MF game, Fig. 1(a), is the flat segment between the (m, n)= (8,5) and (5,3) states at f c ≈0.3-which corresponds to the value of f c in the chaotic regime of NM game. The switch between two regimes is reminiscent of the first-order phase transitions, driven by the competition between two stable states. In our case, it is the switch between two competing states-with (m, n)=(8, 5) and (5, 3)-driven by the mean-field coupling. The mean-field coupling, Eqs. (6)- (7), changes the nature of the transitions between steady state regimes: for the NM game the transitions are discontinuous; for the MF game, the density f c is continuous, but its derivative w.r.t. b jumps. The locations of transitions, b * , are defined by the crossings of constant values f = f * and hyperbolas, Eq. (8). In the vicinity of a transition, expanding Eq. (8) for small |b − b * |, we find where the amplitudes A m,n ∝ 1/(m − n) and the pseudocritical exponent β = 1. Table I lists the numeric values of amplitudes for transitions shown in Fig. 1(a).   Fig. 1(a). See text for discussion.
Spatial chaos.-A salient feature of the NM game is the appearance of the chaotic regime: for 8/5<b<5/3 all agents change their strategy with typical time scales of several game rounds while maintaining the mean density constant, f c ≈0.3 [26]. (For the NM game with selfinteractions, the chaotic regime occurs for 9/5<b<2 with a slightly larger density, f c ≈12 ln 2−8≈0.318 [14]).
The standard way of probing the chaotic behavior is via the so-called persistence [27][28][29] P (t 1 , t 2 ), which in the present context is nothing but the fraction of players who never change their strategy between time steps t 1 and t 2 . Note that we are primarily interested in the asymptotic steady state persistence, P (t 1 , t 2 →∞). To this end, we compute P (τ, T ) with τ and T having the same meaning as in Eq. (4). Averaging this value over the realizations of initial conditions, we obtain an estimate of the asymptotic value, P (0, ∞). We stress that P (τ, ∞) is qualitatively different from P (0, ∞): The latter is sensitive to the overlap between the steady state and the initial distribution of strategies at t=0. Fig. 1(b) compares the asymptotic persistence for the NM and MF games. Two kinds of regimes are seen for both games: those with P (τ, ∞) =0 -where the steadystate features finite clusters (or "cliques") of agents who cooperatively maintain a fixed strategy, and chaotic regimes where P (τ, ∞)=0 -where everyone changes their strategy at least once.
The dependence of P (τ, ∞) for the NM and MF games is very similar, and the main difference is observed for 3/2<b<5/3: The NM game displays chaotic steady states for 8/5<b<5/3 as expected; for the MF game the chaotic regime is shifted to 1.53 b 1.6, which is exactly the range where the mean density of cooperators, f c , follows Eq. (8).
In other words, the effect of the mean-field coupling, Eqs. (6)- (7), is drastic for steady states of the NM game, which are stable but is parametrically close in b to chaotic states: the MF coupling drives such stationary states into spatial chaos.
We also show in Fig. 1(b) the behavior of P (0, ∞) for the MF game. For b<1, about 40% of agents keep their initial strategy. Upon increasing b, P (0, ∞) jumps to below ∼20% at b=1 and then further decreases before dropping off to zero at the onset of the chaotic regime. Note that P (0, ∞) keeps being essentially zero for b>1.63 where P (τ, ∞) is non-zero. This means that while the finite fraction of agents stabilize, the steady-state of an agent is not related to the initial state at t=0, which is forgotten over the first <τ time steps.. Blue (light grey) color corresponds to the cooperator state C, red (dark grey) color is the defector state D, green color is a C which was a D at the previous time step, and yellow is a D which was a C.
Game field configurations. -Fig. 2 shows representative snapshots of steady states of the game field of the MF game depending on the payoff parameter b. For b<1.53 typical configurations feature web-like structures of D players. These webs are random, i.e., precise positions of the branches vary depending on an initial state of the field(see Fig. 2(a)). In the steady state, these webs remain mostly static, with possible "blinking" players at the interfaces. The widths of the branches of defectors increases with increasing the payoff parameter b, and at around b≈1.53-which is the onset of the chaotic regime-the D webs melt into collections of smaller clusters of C and D, which are no longer static and grow, shrink and collide instead.
At b>1.53 field is dominated with clusters of C and D of varying shapes and sizes ( Fig. 2(b)). Upon further increase of the payoff parameter, for 1.6<b<1.62 the field is dominated by square-shaped clusters of C(see Fig. 2(c)). Isolated clusters grow in all directions, and disintegrate upon colliding with neighboring clusters. For b>1.63 the fields feature static collections of small squareshaped clusters of C embedded into the D background(see Fig. 2(d)).
The rectangular shapes of spatial patterns that emerge in Fig. 2 can be readily understood from the synchronous evolution rules and the payoff structure of the NM game. Consider an interface between regions of C and D agents. The interface can be either straight (i.e., oriented along the square lattice directions), or "angled" (e.g., oriented along the diagonal of the lattice). For the NM game (λ = 0 in Eq. (7)), a square cluster of C is stable for 5/3 < b < 8/3; however the "angled" interface is only stable for 7/3 < b < 8/3. For b < 7/3, the angled interface evolves into a straight one in several time steps-and the straight interface is stable. Introducing the meanfield term, Eq. (7) changes switch-over ratios, m/n, into density-dependent hyperbolas, Eq. (8).
It is instructive to examine the distribution of the sizes of clusters of C and D in different steady-state regimes. Figure 3 shows the exponential decay in the distribution of the cluster sizes in chaotic regimes, b=1.55 and 1.59, while it is slower than exponential for b=1.2, the regime with the long memory of initial states (see Fig. 1(b)).
Asymptotically, the fractal dimension of the clusters as well as of the interfaces between clusters of cooperators C and defectors D tends to the dimension of the plane, d=2, which also coincides with the scaling dimension of masses of clusters [18]. These interfaces form a special type of random fractals, quite different from those at the second-order or first-order phase transitions [30].
Discussion.-We propose a way of introducing selfconsistent mean-field-like couplings into evolutionary games with spatially structured populations. We perform direct numerical simulations of a simplest yet non-trivial model-where decision rules include non-local information due to the mean-field coupling-and find a variety of steady states: chaotic states and self-organized stationary states where a finite fraction of agents form stable clusters (or "cliques"). These states are not critical because cluster sizes are distributed exponentially.
We find that the mean-field coupling induces qualitative changes for stationary steady states, which are parametrically close to chaotic states.
In this work, we only consider a restricted version of the Prisoner's Dilemma, Eq. (1). Our approach however, can be readily generalized to a wide class of games, including PD games with non-zero C-D and D-D payoffs and versions of the Public Goods game.
The introduction of the mean-field type couplings opens several intriguing avenues for possible future work. It would be interesting to clarify the effect of the strength of the MF coupling, λ =1. Non-uniform, site-dependent couplings (e.g., λ x drawn from some random distribution) would reflect natural variations between individuals. The role of the local connectivity of the population needs to be clarified: this includes both higherdimensional regular lattices and disordered and hierarchical random lattices and graphs. One may introduce an external field. It would be interesting to see the interplay between the mean-field couplings and noisy decision rules (e.g., where the decision rules contain pseudotemperature [13,15,33]). L.S. is supported within the framework of State Assignment of Russian Ministry of Science and Higher Education. E.B. and A.D. acknowledge support within the Project Teams framework of MIEM HSE. We thank J.J. Arenzon and A. Szolnoki for valuable suggestions. Simulations were carried out in part through computational resources of HPC facilities at HSE University.