Evolutionary dynamics of the prisoner's dilemma with expellers

In pairwise social dilemma games, two players have to simultaneously choose between two strategies, i.e. cooperation (C) and defection (D). Mutual cooperation leads to a reward R, whereas mutual defection results in a punishment value P for each individual (R > P). If one player cooperates and the other one defects, then the cooperator obtains a payoff S, i.e. the sucker's payoff, whereas the defector gets a payoff T, which is often described as the temptation to defect (T > S). In the presence of social dilemmas, mutual cooperation should be preferred over unilateral cooperation (i.e. R > S) and over an equal probability of unilateral cooperation and defection (i.e. 2R > T + S). If T > R > P > S is further satisfied, we obtain the famous prisoner's dilemma [1], which is the most challenging pairwise social dilemma game for the evolution of cooperation [2]. Obviously, it is best for a rational individual to defect regardless of the coplayer's decision when facing such a dilemma. Indeed, defection is the only strict Nash equilibrium [3], and also the evolutionarily stable strategy [2, 4, 5] in the one-short prisoner's dilemma. However, if both players choose defection then their individual income is lower than that obtained for mutual cooperation when maximum overall payoff is shared equally. Meanwhile, we observe that cooperation is ubiquitous in nature as well as human society [6–8]. The evolutionary puzzle here arising is how natural selection can lead to cooperation [9].

Until now, five fundamental rules, i.e. kin selection [10, 11], direct reciprocity [12], indirect reciprocity [13], network reciprocity [14–16] and group selection [17, 18], are reported to be able to enforce the evolution of cooperation in different contexts [19]. Of the five mechanisms for the evolution of cooperation, considerable attention from physics community has been paid to network reciprocity [20–22]. In evolutionary games on graphs, players occupy the vertices of a graph. The edges determine who interacts and competes with whom. Following the pioneering work of Nowak and May [14], evolutionary prisoner's dilemma has been explored on various spatial networks [23–25] and different social graphs [26–33]. Nowak and May show that spatial prisoner's dilemma with best-take-over rule (deterministic evolutionary rule) and synchronous updating pattern (discrete-time evolution) can generate chaotically changing patterns by which cooperators are able to coexist with defectors indefinitely [14]. Szabó and Tőke subsequently extend spatial prisoner's dilemma with Fermi-function-style rule (stochastic evolutionary rule) and asynchronous updating pattern (continuous-time evolution) [23], and observe that such stochasticity just simplifies the evolutionary dynamics, but does not alter the basic finding that cooperators and defectors can persist permanently in spatial networks. Masuda and Aihara find that small-world topology realizes rapid convergence to the equilibrium with slightly suboptimal level of cooperation in comparison with spatial networks [27]. Moreover, Santos and Pacheco reveal that scale-free networks provide the most favorable conditions for the evolution of cooperation [28]. Recently, Gómez-Gardeñes et al show that multiplex networks enlarge the coexistence regions of cooperators and defectors: Cooperation is promoted in parameter regions where it cannot survive in the monoplex scenario at the expense of a moderate decrease of cooperation in those where it can dominate [31]. On the whole, the survival or prevalence of cooperators on graphs is mainly due to the formation of network clusters by which cooperators are able to protect themselves against invasion of defectors. Meanwhile, a large amount of works dedicate to studying cooperation dynamics by the interplay between network reciprocity and other additional mechanisms [34–45]. For instance, Perc et al discuss how the annealed [34] and quenched [35] randomness of payoff affects the evolution of cooperation in spatial prisoner's dilemma. They find that the annealed randomness of payoff can be beneficial for the evolution of cooperation in a way similar to coherence resonance within the framework of noise-driven dynamical systems [46], while the quenched one promotes the evolution of cooperation owing to social diversity: high-ranking wealthy players can form robust cooperative clusters with low-ranking obedient neighbors. Szolnoki and Perc investigate the spatiotemporal dynamics of cooperation in spatial networks with multiple populations, and find the sustenance of cooperation under adverse conditions that could never be bridged by network reciprocity alone as well as the spontaneous emergence of cyclic dominance, where defectors of one population become prey of cooperators in another population, and vice versa [45]. In recent years, there has also been a growing interest in evolutionary games on dynamic networks [47–52]. In this framework, individuals are able to not only change their strategies but also adjust their social ties in response to game outcomes. The main finding from these studies is that the coevolutionary dynamics of individual strategies and network structures can lead to robust cooperation in social networks. Besides the five fundamental mechanisms, voluntary participation [53, 54], costly punishment [55, 56] and destruction effects [57] are also reported to support the evolution of cooperation in different situations.

On the other hand, expulsion behavior is frequently observed in amphibians, reptiles and mammals mainly for the purpose of resource occupation [58]. For example, the red-backed salamanders, Plethodon cinereus, use pheromonal markers, visual aggressive displays, and bites in expulsion of intruders from areas in which the residents occur [59–61]. The mountain spiny lizards, Sceloporus jarrovi, use aggression to defend their territories from conspecifics by expressing a series of behaviors with increasing intensity: pushups, fullshows, shudders, approaches, face-offs and charges which are usually accompanied by bites [62–64]. The Mongolian gerbils (Meriones unguiculatus) trigger the expulsion of their family members with the highest reproductive capacity to compete for reproduction within the groups [65, 66]. Other examples include parasite expulsion (e.g. the nematode Trichinella spiralis are expelled from gastrointestinal tract by immune system due to the contribution of both mucosal mast and goblet cells, and the mediators they produce [67]), social ostracism (e.g. three-spined sticklebacks prevent other sticklebacks if those individuals emit cues, which indicates that they have been infested by parasites, from joining their social groups [68]), and so on. In fact, expulsion behavior is common in not only nonhuman animal groups but also human society [69–76].

Despite the availability of such abundant empirical data of expulsion behavior, we are still unclear its evolutionary origin as well as its significance for the cooperation among individuals in the presence of social dilemmas. In the following, we construct a simple evolutionary game theoretical model of prisoner's dilemma with expellers. The expellers not only cooperate with coplayers in the game, but also 'punish' defectors by excluding them from future reciprocal relationships. This is motivated by the observation that individuals prefer neighborhoods with better quality, and are often willing to expel cheaters who accept benefits but do not reciprocate in order to avoid exploitation. Remarkably, it is found that expulsion and cooperation can coevolve in the spatial prisoner's dilemma with expellers as long as population has a sparse structure, though defection is the only evolutionarily stable strategy in the mean-field prisoner's dilemma with expellers.

The paper is organized as follows. Section 2 describes the evolutionary game theoretical model of prisoner's dilemma with expellers in detail. Section 3 presents the main results, including analysis of mean-field dynamics (section 3.1) as well as spatial dynamics (section 3.2). Section 4 gives summarizations and discussions on our findings.

Consider a population of size N structured by an undirected and unweighted graph $G\left(V,E\right)$ where the set of vertices V represents the sites, and the set of edges E denotes the possible pairwise symmetric interactions between individuals. The vertex v_i represents the ith site which can be either occupied by an individual or just empty, while the edge e_ij denotes the possible pairwise partnership between the individual at site v_i and that at site v_j.

Initially, a fraction p_ϕ ($0\leqslant {p}_{\phi }\lt 1$) of sites are set to be vacant, and the remaining sties are occupied by either an expeller (E, denoted by a three-dimensional unit vector $s={\left[1,0,0\right]}^{T}$) or a cooperator (C, $s={\left[0,1,0\right]}^{T}$) or a defector (D, $s={\left[0,0,1\right]}^{T}$). In each time step (e.g. t), pairwise interaction dynamics happens with probability 1−w, whereas strategy updating dynamics occurs with probability w.

2.1. Pairwise interaction dynamics

A pair of sites linked by an edge, e.g. e_ij, is randomly picked up from the graph. If both sites are occupied by individuals, they play the one-shot prisoner's dilemma with each other; otherwise, the current time step ends. The payoff P_ij of the individual at site v_i and the payoff P_ji of the individual at site v_j resulted from such an interaction can be respectively expressed as

$$\begin{eqnarray}&&\left\{\begin{array}{l}{P}_{{ij}}={{s}_{i}}^{T}{{Ms}}_{j},\\ {P}_{{ji}}={{s}_{j}}^{T}{{Ms}}_{i},\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where the 3 × 3 payoff matrix M is given by

$$\begin{eqnarray}M=\begin{array}{lcccc} & & E & C & D\\ E & & 1-r & 1-r & -r\\ C & & 1-r & 1-r & -r\\ D & & 1 & 1 & 0\end{array},\end{eqnarray} \tag{ 2 }$$

which represents the donor-recipient version of prisoner's dilemma with a single rescaled parameter: the cost-to-benefit ratio r (0 < r < 1). Here the cost-to-benefit ratio r measures the strength of dilemma confronted by individuals [77, 78]. Note that, in our model, expellers act as cooperators to play the prisoner's dilemma with their coplayers. But different from cooperators, we assume that expellers will unilaterally exclude defectors at a cost c from their current sites into any other vacant sites in the population after they are exploited by defectors in the prisoner's dilemma.

The payoff of each individual, e.g. the payoff P_i of the individual at site v_i, is accumulated during the pairwise interaction process. Because of the randomness in selecting pairs of individuals, the rounds that each individual participates in vary from each other. To decouple the effect of this heterogeneity, the accumulative payoff of each individual, e.g. the one at site v_i, is finally normalized according to

$$\begin{eqnarray}&&{\bar{P}}_{i}=\left\{\begin{array}{l}{P}_{i}/z\ \mathrm{if}\ z\gt 0,\\ 0\ \mathrm{if}\ z=0,\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where z represents the frequency with which the individual at site v_i is chosen for interaction with other ones during the pairwise interaction process.

2.2. Strategy updating dynamics

Individuals experience strategy updating phase in a synchronous manner. Specifically, each individual at site v_i imitates a randomly selected neighbor at site v_j, if any, with probability given by the Fermi function [20]:

$$\begin{eqnarray}&&{Tr}\left({\bar{P}}_{j}-{\bar{P}}_{i}\right)=\displaystyle \frac{1}{1+\exp \left[-\left({\bar{P}}_{j}-{\bar{P}}_{i}\right)/K\right]},\end{eqnarray} \tag{ 4 }$$

where the noise K is the measure of stochastic uncertainties allowing irrational choices. Otherwise, the player at site v_i does not change its strategy.

Obviously, the parameter w controls time scale between pairwise interaction process and strategy updating process: When w → 0, pairwise interaction dynamics proceeds much faster than strategy updating dynamics. Therefore, the equilibrium of pairwise interaction dynamics determines the local configurations as well as the average payoffs of individuals on a graph (see an analytical method for regular graphs in appendix B). This means that the consequent strategy evolution proceeds as on a graph with no mobility of individuals induced by expulsion events (see a theoretical method for regular graphs in supplementary material available online at stacks.iop.org/JPCO/3/015011/mmedia); when w → 1, one can expect neutral evolution on a graph depending on the initial strategy configuration; whenever 0 < w < 1, evolution is driven by a detailed interplay between these two dynamical processes.

3.1. Mean-field dynamics

Now let us consider the deterministic dynamics of prisoner's dilemma with expellers in the mean-field setting, where the state is characterized by density p of players (i.e. p_E of expellers, p_C of cooperators and p_D = 1 − p_E − p_C − p_ϕ of defectors) (see appendix A). In the limit w → 0, the expected values ${Ex}\left(\bullet \right)$ of average payoffs ${\bar{P}}_{X}$ (X ∈ E, C, D) for expellers, cooperators and defectors are respectively given by

$$\begin{eqnarray}&&\left\{\begin{array}{l}{Ex}\left({\bar{P}}_{E}\right)=\left(1-r\right)\left({p}_{E}+{p}_{C}\right)-\left(r+c\right){p}_{D},\\ {Ex}\left({\bar{P}}_{C}\right)=\left(1-r\right)\left({p}_{E}+{p}_{C}\right)-{{rp}}_{D},\\ {Ex}\left({\bar{P}}_{D}\right)={p}_{E}+{p}_{C}.\end{array}\right.\end{eqnarray} \tag{ 5 }$$

According to the replicator-like equation set describing motion for the concentration of each type of players (equation (A.5)), p_E and p_C tend to zero for arbitrary value of p_ϕ and K as ${Ex}\left({\bar{P}}_{D}\right)\gt {Ex}\left({\bar{P}}_{C}\right)\geqslant {Ex}\left({\bar{P}}_{E}\right)$. The increment of p_ϕ or K merely leads to the slowdown of evolutionary speed (figure 1), but cannot change evolutionary direction of the population as ${Ex}\left({\bar{P}}_{D}\right)\gt {Ex}\left({\bar{P}}_{C}\right)\geqslant {Ex}\left({\bar{P}}_{E}\right)$ for any composition of the population (equation (5)). In short, expellers and cooperators become extinct in the mean-field limit for w → 0 (figure 2(a)). In the opposite limit w → 1, pairwise interaction dynamics almost never occurs. In this case, one can obtain ${Ex}\left({\bar{P}}_{E}\right)={Ex}\left({\bar{P}}_{C}\right)={Ex}\left({\bar{P}}_{D}\right)=0$. Then natural selection will not change the composition of the population. Any mixture of expellers, cooperators and defectors in the simplex S₃ is an unstable equilibrium (figure 2(b)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the general case of 0 < w < 1, we cannot give analytical predictions on the evolutionary outcome of the population since the expressions of ${Ex}\left({\bar{P}}_{E}\right)$, ${Ex}\left({\bar{P}}_{C}\right)$ and ${Ex}\left({\bar{P}}_{D}\right)$ cannot be further simplified (equations (A.1)–(A.3)). Nevertheless, in the mean-field limit, the prisoner's dilemma with expellers can be described by the following transformed 3 × 3 payoff matrix M':

$$\begin{eqnarray}M^{\prime} =\begin{array}{lcccc} & & E & C & D\\ E & & 1-r & 1-r & -c-r\\ C & & 1-r & 1-r & -r\\ D & & 1 & 1 & 0\end{array},\end{eqnarray} \tag{ 6 }$$

according to equations (A.1)–(A.3) From the above payoff matrix, we can see that the strategy D is completely dominant over C and E, and has a basin of attraction extended to the whole domain of the system. Based on this fact, we predict that the variation of w does not alter the evolutionary outcome that defectors will eventually wipe out expellers and cooperators irrespective of the initial composition of the population. To confirm this prediction, we present phase diagrams of prisoner's dilemma with expellers for w = 0.2, w = 0.4, w = 0.6 and w = 0.8 in figure 3. It can be seen that all evolutionary trajectories end up in the pure state of defectors irrespective of the values of w. In sum, defectors are advantageous over expellers as well as cooperators, and can take over the whole population in the mean-field limit.

Zoom In Zoom Out Reset image size

3.2. Spatial dynamics

In realistic populations, unlike the mean-field model, individuals do not interact with all other ones. Let us therefore consider a structured population with individuals arranged on a square lattice with von Neumann neighborhood (i.e. including connections between nearest neighbors) and periodic boundary conditions. Despite their difference from actual social networks [79], lattices provide a convenient entry point for exploring the impacts of network structures on evolutionary dynamics [14, 20, 80]. Furthermore, the interaction and competition between species in realistic systems, especially in biology and ecology, can be represented adequately by means of lattices [81–84]. Therefore, we aim to study the evolutionary dynamics of prisoner's dilemma with expellers on spatial networks.

In our computer simulations, each efficient time step (t_e) gives exactly one chance for every individual to change its strategy. Initially, an equal number of expellers, cooperators and defectors are randomly placed on a square lattice of size N, such that the fraction of occupied sites is 1 − p_ϕ. The stationary densities of expellers (p_E), cooperators (p_C) and defectors (p_D) on square lattices are evaluated over sufficiently long sampling times after the system enters into the equilibrium state. To ensure proper accuracy, final results (i.e. the mean densities of expellers ${\bar{p}}_{E}$, cooperators ${\bar{p}}_{C}$ and defectors ${\bar{p}}_{D}$) are averaged over up to 50 independent runs for each set of parameter values.

3.2.1. The limiting case $w \rightarrow 0$ and $w \rightarrow 1$

For convenience of comparison, we start by presenting the evolutionary outcomes of spatial prisoner's dilemma when the parameter of time scale w → 0 in figure 4. It is shown that the mean fraction of cooperators ${\bar{p}}_{C}/\left(1-{p}_{\phi }\right)$ decreases with the cost-to-benefit ratio r when ratio of vacant sites p_ϕ is low (figures 4(a) and (b)). In this case, cooperators are able to survive in the population by forming compact clusters as long as the cost-to-benefit ratio r is sufficiently low [23]. However, the mean fraction of cooperators ${\bar{p}}_{C}/\left(1-{p}_{\phi }\right)$ becomes almost invariant with the cost-to-benefit ratio r when ratio of vacant sites p_ϕ is high (figures 4(c) and (d)). In this case, the clusters of cooperators become sparse. Cooperators are vulnerable to the invasion of defectors irrespective of the values of cost-to-benefit ratio r. Therefore, only cooperators isolated from the spatial domain of defectors remain in the population [85].

Zoom In Zoom Out Reset image size

Let us now study the spatial prisoner's dilemma with expellers when the parameter of time scale w → 0. In comparison with those in figure 4, results reported in figure 5 demonstrate that expulsion is able to promote the evolution of cooperation irrespective of the values of cost-to-benefit ratio r, ratio of vacant sites p_ϕ and cost of expulsion c. On saturated square lattices (i.e. p_ϕ = 0), there is no vacant sites available for expellers to banish their defective neighbors. Expellers, like cooperators, persist in the population solely by virtue of network reciprocity. If cost of expulsion c = 0, expellers and cooperators become equivalent. Hence the mean fraction of expellers ${\bar{p}}_{E}/\left(1-{p}_{\phi }\right)$ is identical to that of cooperators ${\bar{p}}_{C}/\left(1-{p}_{\phi }\right)$ for any value of cost-to-benefit ratio r (figure 5(a)). However, if cost of expulsion c > 0, the evolutionary performance of expellers becomes inferior to that of cooperators (figures 5(e), (i) and (m)). In contrast, on diluted square lattices (i.e. p_ϕ > 0), expellers perform better than cooperators and can persist (or even prevail) in the population across the whole applicable range of cost-to-benefit ratio r and cost of expulsion c (figure 5). In short, expulsion can coevolve with cooperation in the spatial prisoner's dilemma with expellers in contrast to the evolutionary outcomes in the mean-field limit.

Zoom In Zoom Out Reset image size

To systematically investigate the effects of cost of expulsion c and ratio of vacant sites p_ϕ, we move to study phase diagrams for two representative values of cost-to-benefit ratio r in the spatial prisoner's dilemma with expellers. In the absence of expellers, cooperators survive only if r < 0.006, and they are able to defeat defectors completely for r < 0 on the saturated square lattice. Taking these thresholds as benchmark values, we focus on r = 0.005 and r = 0.5. For r = 0.005, cooperators are able to coexist with defectors solely on the basis of network reciprocity. This value of r thus yields lenient conditions for the evolution of cooperation. For r = 0.5, on the other hand, cooperators are unable to survive. This value of r thus yields adverse conditions for cooperative behavior to emerge. In close vicinity of transition points, we here enlarge the size of system from N = 10⁴ (figures 4 and 5) to 9 × 10⁴ (figure 7) to determine boundaries and types of phase transitions with a higher accuracy. Meanwhile, the parameter of time scale w should be covariant with the size of system N to ensure invariance of the ratio of time scale between pairwise interaction dynamics and strategy updating dynamics. Namely,

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{1}{w}_{{N}_{1}}}{1-{w}_{{N}_{1}}}=\displaystyle \frac{{N}_{2}{w}_{{N}_{2}}}{1-{w}_{{N}_{2}}},\end{eqnarray} \tag{ 7 }$$

where w_N denotes the parameter of time scale associated with a system of size N. Figure 6 shows the dynamical behaviors of two systems with different sizes (i.e. N = 10⁴ and N = 9 × 10⁴) but the same ratio of time scale (i.e. Nw/(1 − w) = 5 × 10⁻³). Clearly, evolutionary trajectories for both systems are largely coincident with each other. Due to the finite-size effects, the current fractions of expellers ${p}_{E}\left({t}_{e}\right)/\left(1-{p}_{\phi }\right)$, cooperators ${p}_{C}\left({t}_{e}\right)/\left(1-{p}_{\phi }\right)$ and defectors ${p}_{D}\left({t}_{e}\right)/\left(1-{p}_{\phi }\right)$ fluctuate around an equilibrium after initial transient times, wherein fluctuations are smaller in the lattice with larger size (figure 6).

Zoom In Zoom Out Reset image size

The phase diagram for r = 0.005 presented in figure 7(c) together with figure 5 suggests that at such a low value of r, expellers rather than cooperators represent the more effective players to outperform defectors except if cost of expulsion c > 0 and ratio of vacant sites p_ϕ = 0. The mixed and dynamical (C + D)_D phase in the bottom line of the c − p_ϕ plane first gives way to the mixed and dynamical (E + C + D)_D phase, subsequently to the mixed and frozen (E + D)_F phase and finally to the mixed and frozen (E + C + D)_F phase as the ratio of vacant sites p_ϕ increases. Here a mixed state is further classified as a frozen or dynamical sub-state. A sub-state is called to be frozen provided two or more of the second order heterogeneous variables (i.e. p_EC, p_ED or p_CD) of the system are equal to zero in the stationary state. Otherwise, the sub-state is defined to be dynamical. While majority of phase transitions is continuous, the ${\left(E+C+D\right)}_{D}\to {\left(E+D\right)}_{F}$ phase transition is discontinuous because of an indirect territorial competition between expellers and cooperators (for details see main text below) [86, 87]. Only if the ratio of vacant sites p_ϕ is either negligible or large, are cooperators able to survive. When the ratio of vacant sites p_ϕ is sufficiently small (e.g. p_ϕ = 0.005), the clusters of cooperators are compact enough to fight against defectors (figure 8(a)). Accordingly, cooperators are able to coexist alongside expellers and defectors in a dynamical manner (see the mixed and dynamical (E + C + D)_D phase in lower part of figure 7(c)). When the ratio of vacant sites p_ϕ is large (e.g. p_ϕ = 0.8), the clusters of cooperators are sparse and thus are easily invaded by their defective neighbors. In this case, only cooperators that are isolated from defectors can persist in the system (figure 8(b)). Therefore, cooperators coexist with expellers and defectors in a frozen manner (see the mixed and frozen (E + C + D)_F phase in top part of figure 7(c)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Interestingly, when the ratio of vacant sites p_ϕ is moderate (e.g. p_ϕ = 0.1), cooperators become extinct, and expellers and defectors are frozenly coexistent in the system (see the mixed and frozen (E + D)_F phase in middle part of figure 7(c)). To understand why expellers can outperform cooperators despite additional costs of expulsion, it is important to analyze spatiotemporal dynamics of spatial prisoner's dilemma with expellers. Let ${{Tr}}_{{XY}}\left({t}_{e}\right)\left(X,Y\in E,C,D\right)$ denote the transition probability of X's changing into Y's at efficient time step t_e. Defining $p{{\prime} }_{X}\left({t}_{e}\right)={p}_{X}\left({t}_{e}\right)/\left(1-{p}_{\phi }\right)\left(X\in E,C,D\right)$, we get the following equation set [20]:

$$\begin{eqnarray}&&\left\{\begin{array}{l}\tfrac{{dp}{{\prime} }_{E}\left({t}_{e}\right)}{{{dt}}_{e}}={G}_{{CE}}\left({t}_{e}\right)-{G}_{{EC}}\left({t}_{e}\right)+{G}_{{DE}}\left({t}_{e}\right)-{G}_{{ED}}\left({t}_{e}\right),\\ \tfrac{{dp}{{\prime} }_{C}\left({t}_{e}\right)}{{{dt}}_{e}}={G}_{{EC}}\left({t}_{e}\right)-{G}_{{CE}}\left({t}_{e}\right)+{G}_{{DC}}\left({t}_{e}\right)-{G}_{{CD}}\left({t}_{e}\right),\\ \tfrac{{dp}{{\prime} }_{D}\left({t}_{e}\right)}{{{dt}}_{e}}={G}_{{ED}}\left({t}_{e}\right)-{G}_{{DE}}\left({t}_{e}\right)+{G}_{{CD}}\left({t}_{e}\right)-{G}_{{DC}}\left({t}_{e}\right),\end{array}\right.\end{eqnarray} \tag{ 8 }$$

where ${G}_{{XY}}\left({t}_{e}\right)=p{{\prime} }_{X}\left({t}_{e}\right){{Tr}}_{{XY}}\left({t}_{e}\right)\left(X,Y\in E,C,D\right)$ represents the current inflow of Y's transitioning from X's at efficient time step t_e. In our analysis, both ${p}_{X}\left({t}_{e}\right)\left(X\in E,C,D\right)$ and ${{Tr}}_{{XY}}\left({t}_{e}\right)\left(X,Y\in E,C,D\right)$ are obtained from computer simulations. Starting from a random strategy distribution (figure 9(e)), the imitation of better performing neighbors generates clusters of individuals with the same strategies: expellers and cooperators form tiny clusters in a sea of defectors (figure 9(f)). In this most early stage of evolutionary process (i.e. ${t}_{e}\in \left[1,8\right]$), the description of mean-field is applicable, and the current fraction of expellers ${p}_{E}\left(t\right)/\left(1-{p}_{\phi }\right)$ decreases with the efficient time t_e in an exponential way due to the advantage of defectors and cooperators in payoff (figures 9(a) and (i)). In comparison, the current fraction of cooperators ${p}_{C}\left(t\right)/\left(1-{p}_{\phi }\right)$ does not vary sharply with the efficient time t_e: it first increases mildly and then decreases slowly. Although exploited by defectors at this stage, cooperators are compensated at the expense of expellers (e.g. at t_e = 1, the net inflow of cooperators changing from defectors ${G}_{{DC}}\left({t}_{e}\right)-{G}_{{CD}}\left({t}_{e}\right)\approx -0.042\,7;$ the net inflow of cooperators changing from expellers ${G}_{{EC}}\left({t}_{e}\right)-{G}_{{CE}}\left({t}_{e}\right)\approx 0.048\,3;$ the overall net inflow of cooperators ${G}_{{DC}}\left({t}_{e}\right)-{G}_{{CD}}\left({t}_{e}\right)+{G}_{{EC}}\left({t}_{e}\right)-{G}_{{CE}}\left({t}_{e}\right)\,\approx 0.005\,6;$ figure 9(i)). The further development is determined by a local equilibrium where strategy updating only occurs at the interfaces between clusters of different strategies (i.e. ${t}_{e}\in \left(8,+\infty \right);$ figures 9(f)–(h)). In the presence of defectors, the evolutionary fate of expellers is decided not only by a direct competition with cooperators, but also by the success of both cooperative strategies against invasion attempts by defectors (the accumulated net inflow of expellers by the manner of direct domain competition with cooperators ${\int }_{{t}_{e}}{G}_{{CE}}\left(t\right)-{G}_{{EC}}\left(t\right){dt}\approx 0.104\,6;$ the accumulated net inflow of expellers through the manner of indirect domain competition with cooperators ${\int }_{{t}_{e}}{G}_{{CDE}}\left(t\right)-{G}_{{EDC}}\left(t\right){dt}\approx 0.226\,2$, where ${G}_{{XYZ}}\left({t}_{e}\right)\left(X,Y,Z\in E,C,D\right)$ denotes the current net inflow of Z's from X's via Y's at efficient time step t_e; figure 9(d)). If ratio of vacant sites p_ϕ is moderate, the clusters of cooperators are not compact enough to defend against defectors in their neighborhood. Furthermore, the spatial configurations of cooperators are constantly subjected to additional destruction due to the attack of migratory defectors induced by the expulsion of expellers (figures 9(f) and (g)). In contrast, expellers are able to continuously banish defectors in their neighborhood, and thus reduce frequency of interaction with defective players (figures 9(f) and (g)). Consequently, the evolutionary performance of cooperators become inferior to that of expellers, which explains the advantage of expellers over cooperators in the direct domain competition (e.g. at t_e = 8, the net inflow of expellers changing from cooperators ${G}_{{CE}}\left({t}_{e}\right)-{G}_{{EC}}\left({t}_{e}\right)\approx 0.001\,1;$ at t_e = 166, the net inflow of expellers changing from cooperators ${G}_{{CE}}\left({t}_{e}\right)-{G}_{{EC}}\left({t}_{e}\right)\approx 0.001\,7;$ figures 9(j) and (k)). Figures 9(b) and (c) demonstrate the excellent competitiveness of expellers against defectors in their evolutionary race. Except extremely few efficient time steps, the current net inflow of expellers changing from defectors is always no less than that of defectors changing from expellers after the self-organization process of initial clusters (e.g. at t_e = 8, the net inflow of expellers changing from defectors ${G}_{{DE}}\left({t}_{e}\right)-{G}_{{ED}}\left({t}_{e}\right)\approx 0.000\,2;$ at t_e = 166, the net inflow of expellers changing from defectors ${G}_{{DE}}\left({t}_{e}\right)-{G}_{{ED}}\left({t}_{e}\right)\approx 0.001\,5;$ figures 9(j) and (k)). Indeed, expellers are able to accumulate more payoffs than defectors in the pairwise interaction process, due to evolutionary mechanisms of network reciprocity as well as expulsion. Interestingly, expellers can spread successfully in the presence of defectors that spatial domains lost by cooperators are finally occupied by expellers (figures 9(f)–(h)), which elaborates the advantage of expellers over cooperators in the indirect domain competition. The synergetic combination of direct and indirect territorial competitions ultimately leads to the eventual extinction of cooperators and the final downfall of defectors (figure 9(a)). The evolutionary competition between expellers and defectors terminates when they are completely separated from each other (figures 9(h) and (l)).

Zoom In Zoom Out Reset image size

Moreover, if ratio of vacant sites ${p}_{\phi }$ is negligible, the increment of cost of expulsion c can lead to phase transition from the mixed and frozen (E + D)_F phase to the mixed and dynamical (E + C + D)_D phase (figure 7(c)). Here the emergence of cooperation is due to the increasing cost of expulsion c, which lends more support for cooperators to compete with expellers. The more specific nature of the phase diagram for r = 0.005 is illustrated quantitatively in figures 10(a) and (b), which shows a cross section across the ratio of vacant sites p_ϕ for cost of expulsion c = 2 and that across the cost of expulsion c for ratio of vacant sites p_ϕ = 0.02, respectively. Interestingly, it can be observed in figure 10(a) that there exists an optimal ratio of vacant sites ${{p}_{\phi }}^{* }$ leading to the maximal overall fraction of altruistic individuals (i.e. expellers and cooperators) in the system. Note that the optimal ratio of vacant sites ${{p}_{\phi }}^{* }$ in present study does not relate to percolation threshold of the host graph, in contrast to the findings in either prisoner's dilemma [85] or public goods game [88]. Here the evolutionary fate of altruistic individuals is determined by both direct and indirect territorial competition between expellers and cooperators. Neither too large nor too small but moderate ratios of vacant sites ${{p}_{\phi }}^{* }$ maximize the positive impact of territorial competition.

Zoom In Zoom Out Reset image size

If the conditions for the evolution of cooperation become harsh, as is the case for r = 0.5, the phase diagram changes qualitatively (figure 7(d)). It is shown that the mixed and dynamical (C + D)_D phase at ratio of vacant sites p_ϕ = 0 as well as the mixed and dynamical (E + C + D)_D phase at low ratio of vacant sites p_ϕ transforms into the pure D phase and the mixed and frozen (E + D)_F phase, respectively. Here, cooperators vanish because of the increasing cost of cooperation. In addition, the parameter region where expellers, cooperators and defectors can frozenly coexist (i.e. the mixed and frozen (E + C + D)_F phase) extends to lower part of the c − p_ϕ parameter plane (compare figures 7(d) with (c)). When the cost-to-benefit ratio r is large (e.g. r = 0.5), expellers and cooperators can survive only if they are separated from defectors. If the ratio of vacant sites p_ϕ is small (e.g. p_ϕ = 0.1), cooperators would not have enough space to escape from the exploitation by defectors. Thus expellers coexist with defectors in a frozen manner (figure 11(a)). If the ratio of vacant sites p_ϕ is large (e.g. p_ϕ = 0.8), as is the case for r = 0.005, cooperators are able to frozenly coexist with expellers and defectors in the system (figure 11(b)). Figures 10(c) and (d) feature two characteristic cross sections of the phase diagram presented in figure 7(d), which span across the ratio of vacant sites p_ϕ for cost of expulsion c = 2 and the cost of expulsion c for ratio of vacant sites p_ϕ = 0.2, respectively. Interestingly, it can be seen from figure 10(d) that the impacts of cost of expulsion c become saturated when it is increased: the mean fractions of players are invariant with the cost of expulsion c.

Zoom In Zoom Out Reset image size

Generally, the results obtained from Monte Carlo simulations have been confirmed by using the extended technique of pair approximation (see appendix B and supplementary material; figures 4, 5, 7 and 10). In terms of equilibrium frequencies of strategies, the extended pair approximation methods correctly predict the trends of simulation data in variation with relevant parameters (figures 4, 5 and 10). However, predictions by the extended pair approximation approaches are less accurate in the vicinity of phase transition points where long-range correlations and fluctuations play a dominant role (compare bottom row with top row of figure 7) [89]. In principle, more reliable predictions can be made by enlarging the basic cluster from two-point (used by pair approximation or two-point approximation) to n-point (n > 2; used by n-point approximation) [20].

In the opposite limit w → 1, the system will stay in the strategy updating process but rarely enter into the pairwise interaction process. In such a case, the whole population is governed by random drift because the payoffs of expellers, cooperators and defectors are approximately equal to zero. In terms of game results, there is no role difference between expellers, cooperators and defectors. Hence, expulsion plays neither positive nor negative role in the evolution of cooperation.

3.2.2. The general case $0\lt w\lt 1$

To explore impacts of parameter of time scale w on the evolutionary outcome of spatial prisoner's dilemma with expellers in detail, we present the contour plots of mean fraction of expellers ${\bar{p}}_{E}/(1-{p}_{\phi })$, cooperators ${\bar{p}}_{C}/(1-{p}_{\phi })$ and defectors ${\bar{p}}_{D}/(1-{p}_{\phi })$ as a function of cost-to-benefit ratio r and parameter of time scale w for different values of ratio of vacant sites p_ϕ in figure 12. On saturated square lattices (i.e. p_ϕ = 0), cooperators are able to survive only if cost-to-benefit ratio r and parameter of time scale w are sufficient small (figure 12(b)). Otherwise, defectors dominate expellers as well as cooperators completely (see top row of figure 12). However, on diluted square lattices (i.e. p_ϕ > 0), expellers can coexist with defectors on the full r − w parameter plane (see middle and bottom rows of figure 12). For relatively low ratio of vacant sites p_ϕ (e.g. p_ϕ = 0.4), expellers perform best when cost-to-benefit ratio r and parameter of time scale w are negligible (figure 12(e)). Meanwhile, cooperators are able to survive only if the parameter of time scale w is moderate (figure 12(f)). Interestingly, the mean fraction of cooperators ${\bar{p}}_{C}/\left(1-{p}_{\phi }\right)$ is largely invariant with cost-to-benefit ratio r. For relatively high ratio of vacant sites p_ϕ (e.g. p_ϕ = 0.8), the evolutionary fate of expellers is unaffected by cost-to-benefit ratio r as well as parameter of time scale w (figure 12(i)).As for cooperators, both cost-to-benefit ratio r and parameter of time scale w have similar effects as they have in the case of relatively low ratio of vacant sites p_ϕ (compare figures 12(f) with (j)). Intriguingly, it is shown that there exists an optimal as well as a worst value of time scale parameter w⁺ and w⁻that leads to the highest and lowest level of expulsion and cooperation, respectively (figures 12(c), (g) and (k)). In present study, pairs of individuals are stochastically chosen to play prisoner's dilemma with expellers in the pairwise interaction process. The randomness of this process can be adjusted by the parameter of time scale w. In one limit w → 0, the pairwise interaction dynamics is deterministic. In the other limit w → 1, the pairwise interaction dynamics becomes stochastic to the largest extend. Whenever 0 < w < 1, the pairwise interaction dynamics falls somewhere between these limits. Therefore, the existence of such an optimal phenomenon can attribute to the coherence resonance scenario trigged by the parameter of time scale w which may be considered as the constructive noise [90, 91]. In this case, the increment of mean fraction of altruistic players (i.e. expellers and cooperators) $\left({\bar{p}}_{E}+{\bar{p}}_{C}\right)/\left(1-{p}_{\phi }\right)$ is regarded as an constructive effect, as widespread altruistic players yield a higher total population payoff in comparison with dominant defectors, and thus is favourable for the population. The mean fraction of altruistic players $\left({\bar{p}}_{E}+{\bar{p}}_{C}\right)/\left(1-{p}_{\phi }\right)$ determines the constructive effects of noise on the system and thus has the same meaning as the signal-to-noise ratio in classical coherence resonance phenomena observed in dynamical systems [46]. On the other hand, the appearance of the worst value of time scale parameter w⁻near w = 1 is due to the unbalanced interplay between pairwise interaction dynamics and strategy updating dynamics. Figure 13 shows the optimal and worst parameter of time scale w⁺ and w⁻as a function of ratio of vacant sites p_ϕ. Both the optimal and worst parameters of time scale w⁺ and w⁻are increased with the ratio of vacant sites p_ϕ. This result indicates that with increasing vacancy rate the optimal as well as worst trade-offs between pairwise interaction dynamics and strategy updating dynamics translate from deterministic to stochastic manner.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The purpose of the present paper is to study the evolutionary dynamics of prisoner's dilemma with expellers, and to determine whether expulsion can provide evolutionary advantages for altruistic individuals (i.e. expellers and cooperators) to compete with defectors in prisoner's dilemma. In the mean-field limit, we have shown that defection is the only evolutionarily stable strategy in prisoner's dilemma with expellers. The increment of either vacant sites ratio or time scale parameter between pairwise interaction dynamics and strategy updating dynamics merely leads to the extension of mean time to reach the stable equilibrium. However, expulsion is, in general, effective in deterring defection in spatial prisoner's dilemma with expellers as long as the population has a sparse structure. In one limit of time scale parameter w → 0, the introduction of expellers can dramatically change the system behavior, wherein we do not find the dominance of defectors as it is in spatial prisoner's dilemma, but do find the maintenance or even prevalence of expellers (compare figures 4 with 5). In the presence of defectors, expellers can take advantages over cooperators by a combination of direct and indirect domain competition. By the evolutionary mechanism of expulsion, expellers are able to decrease frequency of interaction with defectors while increase frequency of interaction between cooperators and defectors (see appendix B). The former factor lowers down the payoffs of defectors in expellers' neighborhood, which helps expellers to reach a superiority over defectors who in turn diminish cooperators: a manner of indirect territorial competition. On the other hand, both factors assist expellers to defeat cooperators in the direct territorial competition. Although there exists a rather constrained and unrealistic parameter region (i.e. cost-to-benefit ratio r → 0, ratio of vacant sites p_ϕ = 0 and cost of expulsion c > 0; figure 5) where cooperators can outperform expellers, these advantages are highly impossible to play a role in reality. The spatial prisoner's dilemma with expellers gives rise to a rich variety of possible dynamics and a number of continuous or discontinuous transitions between qualitatively different system behaviors (figure 7). For low cost-to-benefit ratio r, we have observed several types of strategy coexistence: (1) the mixed and dynamical state ${\left(C+D\right)}_{D}$, (2) the mixed and dynamical state ${\left(E+C+D\right)}_{D}$, (3) the mixed and frozen state ${\left(E+D\right)}_{F}$ and (4) the mixed and frozen state ${\left(E+C+D\right)}_{F}$. Each kind of these states is sustained in a particular way by means of which a portion of altruistic individuals can survive in the presence of defectors. To be specific, cooperators can dynamically coexist with defectors by spatial aggregation in ${\left(C+D\right)}_{D}$ phase. In ${\left(E+C+D\right)}_{D}$ phase, expellers can dynamically coexist with cooperators and defectors by direct and indirect territorial competition, the latter of which also leads to the discontinuous phase transition between ${\left(E+C+D\right)}_{D}$ and ${\left(E+D\right)}_{F}$ phase. In both ${\left(E+D\right)}_{F}$ and ${\left(E+C+D\right)}_{F}$ phases, altruistic individuals can coexist with defectors by spatial isolation in a frozen manner. For high cost-to-benefit ratio r, we have observed more simplified evolutionary outcomes: (1) the pure state D, (2) the mixed and frozen state ${\left(E+D\right)}_{F}$ and (3) the mixed and frozen state ${\left(E+C+D\right)}_{F}$.

Note that our principal discoveries are not expected to change for relaxing parameter of time scale from the limiting case w → 0 to the general case 0 < w < 1 (figure 12). Interestingly, we have found that there exist not only an optimal but also a worst value of time scale parameter w⁺ and w⁻that results in the highest and lowest level of expulsion and cooperation, respectively. By interpreting parameter of time scale w as constructive noise, we can attribute the optimal phenomenon to coherence resonance reported previously in temporal and spatially extended dynamical systems [46, 90, 91]. On the other hand, the existence of worst time scale parameter w⁻is due to the unbalanced interplay between pairwise interaction dynamics and strategy updating dynamics. Furthermore, it is well known that mutation weakens the benefits of network reciprocity to cooperation in social dilemma games [92]. Therefore it is worth investigating how the coevolution of expulsion and cooperation affected by mutation. The preliminary simulation results show that the beneficial effects of expulsive behavior on the evolution of cooperation are fairly robust against the variation of mutation. This indicates that the evolutionary mechanism of expulsion on spatial networks is different from network reciprocity, which deserves additional research to clarify.

It was previously found that punishment represents a typical behavior that is able to promote cooperation in certain situations [93, 94]. Particularly, costly punishment refers to the special kind of punishment that is manipulated by imposing some fine on defectors but at a cost to punishers [86, 95–101]. In fact, this kind of costly punishment can be classified as an active one. In this sense, the expulsive behavior in our model can be considered as a kind of passive punishment: expellers tend to 'punish' defectors in a way that they terminates future interactions between them. On the other hand, expulsion is very similar to ostracism, in which defectors lose all their present partners [102, 103]. In addition, the evolutionary mechanism that assists expellers to fight against the exploitation of defectors share a similarity with that in optional games of social dilemma: reduction of interaction frequency between cooperative strategies and defective strategies [104]. Nevertheless, different from the loner in optional games, expellers in our model can also win over cooperators via enhancing frequency of interaction between cooperators and defectors though it requires additional supports from network reciprocity as well as dilution of graph. This is the reason that leads to different dynamics between games with loners and games with expellers. Still, spatial prisoners dilemma with expellers can also be considered as a kind of coevolutionary model between individual strategies and network structures as the neighborhoods of individuals are also coevolved with strategies in our model [22]. Recently, the impact of social exclusion on the evolution of cooperation governed by group interaction has been extensively studied in infinite populations, finite populations and spatial networks, respectively [105–107]. In their works, excluders are able to deny the rights of defectors to enjoy the benefits of public goods at some cost. In this case, excluders can defeat defectors as long as the net gain of excluders from the public goods game is greater than the cost of exclusion. Thus social exclusion is a very powerful mechanism for the evolution of cooperation even in infinite populations. In contrast, expellers in our model are able to recognize the types of coplayers only after they play games with each other. In other words, expellers need to pay a cost r, which is the cost-to-benefit ratio, to make identification as well as another cost c to execute expulsion so as to expel defectors from their neighborhoods. This is the reason why expellers perform inferior to both cooperators and defectors in the mean-field limit. Finally, cooperators in our model play the role of second-order free riders. To explain this, suppose that there are no defectors in the population. In this case, natural selection cannot distinguish between the players that cooperate and banish defectors (i.e. expellers) and the players that cooperate but do not expel defectors (i.e. cooperators). There is no requirement to expel, so there is no second-order free rider problem. If defectors are present in the population, however, the players that eject defectors must do so at a personal cost. Natural selection will now favor the players that cooperate but do not expel defectors. As a result, cooperators rise at the expense of expellers and eventually introduce the invasion of defectors. However, we have shown that the coevolution of expulsion and cooperation can be promoted even in such an adverse condition.

We thank two anonymous referees for their constructive and insightful comments. XFW is supported by Donghua University Startup Funds for Young Faculty. XFW and WJK gratefully acknowledge financial support from the Fundamental Research Funds for the Central Universities (Grants No. 2232018G-09). WJK is also supported by the National Natural Science Foundation of China (Grants No. 61603088). This research is supported in part by Hong Kong Research Grant Council (Grants No. 15206915 and Grants No. 15208418).

In the mean-field theory, the state of a system is characterized by the density of players. In this situation, the expected values ${Ex}\left(\bullet \right)$ of average payoffs ${\bar{P}}_{X}$ (X ∈ E, C, D) for expellers, cooperators and defectors in prisoner's dilemma with expellers are respectively given by

$$\begin{eqnarray}&&{Ex}\left({\bar{P}}_{E}\right)=w\displaystyle \sum _{I=1}^{\infty }{\left(1-w\right)}^{I}\displaystyle \sum _{{i}_{{EE}}+{i}_{{EC}}+{i}_{{ED}}=1}^{I}\begin{array}{l}\tfrac{I!}{{i}_{{EE}}!\left(I-{i}_{{EE}}\right)!}{{p}_{E}}^{2{i}_{{EE}}}\tfrac{\left(I-{i}_{{EE}}\right)!}{{i}_{{EC}}!\left(I-{i}_{{EE}}-{i}_{{EC}}\right)!}{\left(2{p}_{E}{p}_{C}\right)}^{{i}_{{EC}}}\\ \times \tfrac{\left(I-{i}_{{EE}}-{i}_{{EC}}\right)!}{{i}_{{ED}}!\left(I-{i}_{{EE}}-{i}_{{EC}}-{i}_{{ED}}\right)!}{\left(2{p}_{E}{p}_{D}\right)}^{{i}_{{ED}}}\\ \times {\left(1-{{p}_{E}}^{2}-2{p}_{E}{p}_{C}-2{p}_{E}{p}_{D}\right)}^{I-{i}_{{EE}}-{i}_{{EC}}-{i}_{{ED}}}\\ \times \tfrac{2\left(1-r\right){i}_{{EE}}+\left(1-r\right){i}_{{EC}}+\left(-r-c\right){i}_{{ED}}}{2{i}_{{EE}}+{i}_{{EC}}+{i}_{{ED}}}\end{array},\end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}&&{Ex}\left({\bar{P}}_{C}\right)=w\displaystyle \sum _{I=1}^{\infty }{\left(1-w\right)}^{I}\displaystyle \sum _{{i}_{{CC}}+{i}_{{EC}}+{i}_{{CD}}=1}^{I}\begin{array}{l}\tfrac{I!}{{i}_{{CC}}!\left(I-{i}_{{CC}}\right)!}{{p}_{C}}^{2{i}_{{CC}}}\tfrac{\left(I-{i}_{{CC}}\right)!}{{i}_{{EC}}!\left(I-{i}_{{CC}}-{i}_{{EC}}\right)!}{\left(2{p}_{E}{p}_{C}\right)}^{{i}_{{EC}}}\\ \times \tfrac{\left(I-{i}_{{CC}}-{i}_{{EC}}\right)!}{{i}_{{CD}}!\left(I-{i}_{{CC}}-{i}_{{EC}}-{i}_{{CD}}\right)!}{\left(2{p}_{C}{p}_{D}\right)}^{{i}_{{CD}}}\\ \times {\left(1-{{p}_{C}}^{2}-2{p}_{E}{p}_{C}-2{p}_{C}{p}_{D}\right)}^{I-{i}_{{CC}}-{i}_{{EC}}-{i}_{{CD}}}\\ \times \tfrac{2\left(1-r\right){i}_{{CC}}+\left(1-r\right){i}_{{EC}}+\left(-r\right){i}_{{CD}}}{2{i}_{{CC}}+{i}_{{EC}}+{i}_{{CD}}}\end{array},\end{eqnarray} \tag{ A.2 }$$

and

$$\begin{eqnarray}&&{Ex}\left({\bar{P}}_{D}\right)=w\displaystyle \sum _{I=1}^{\infty }{\left(1-w\right)}^{I}\displaystyle \sum _{{i}_{{DD}}+{i}_{{ED}}+{i}_{{CD}}=1}^{I}\begin{array}{l}\tfrac{I!}{{i}_{{DD}}!\left(I-{i}_{{DD}}\right)!}{{p}_{D}}^{2{i}_{{DD}}}\tfrac{\left(I-{i}_{{DD}}\right)!}{{i}_{{ED}}!\left(I-{i}_{{DD}}-{i}_{{ED}}\right)!}{\left(2{p}_{E}{p}_{D}\right)}^{{i}_{{ED}}}\\ \times \tfrac{\left(I-{i}_{{DD}}-{i}_{{ED}}\right)!}{{i}_{{CD}}!\left(I-{i}_{{DD}}-{i}_{{ED}}-{i}_{{CD}}\right)!}{\left(2{p}_{C}{p}_{D}\right)}^{{i}_{{CD}}}\\ \times {\left(1-{{p}_{D}}^{2}-2{p}_{E}{p}_{D}-2{p}_{C}{p}_{D}\right)}^{I-{i}_{{DD}}-{i}_{{ED}}-{i}_{{CD}}}\\ \times \tfrac{{i}_{{ED}}+{i}_{{CD}}}{2{i}_{{DD}}+{i}_{{ED}}+{i}_{{CD}}}\end{array},\end{eqnarray} \tag{ A.3 }$$

where I represents the frequency of randomly chosen pairs of individuals to endure pairwise interaction dynamics between strategy updating events. i_XY denotes the number of XY pairs satisfying ${\sum }_{{XY}}{i}_{{XY}}=I$. Herein, X and Y stand for the state of a site, which is occupied by either an expeller or a cooperator, or a defector.

Since the present dynamical rule has the form of Fermi function (equation (4)), the motion for concentrations of expellers, cooperators and defectors can be given by the following equation set:

$$\begin{eqnarray}&&\left\{\begin{array}{l}\tfrac{{{dp}}_{E}}{{dt}}=-{p}_{E}{p}_{C}\tanh \left(\tfrac{{Ex}\left({\bar{P}}_{C}\right)-{Ex}\left({\bar{P}}_{E}\right)}{2K}\right)-{p}_{E}{p}_{D}\tanh \left(\tfrac{{Ex}\left({\bar{P}}_{D}\right)-{Ex}\left({\bar{P}}_{E}\right)}{2K}\right),\\ \tfrac{{{dp}}_{C}}{{dt}}=-{p}_{C}{p}_{E}\tanh \left(\tfrac{{Ex}\left({\bar{P}}_{E}\right)-{Ex}\left({\bar{P}}_{C}\right)}{2K}\right)-{p}_{C}{p}_{D}\tanh \left(\tfrac{{Ex}\left({\bar{P}}_{D}\right)-{Ex}\left({\bar{P}}_{C}\right)}{2K}\right),\\ \tfrac{{{dp}}_{D}}{{dt}}=-{p}_{D}{p}_{E}\tanh \left(\tfrac{{Ex}\left({\bar{P}}_{E}\right)-{Ex}\left({\bar{P}}_{D}\right)}{2K}\right)-{p}_{D}{p}_{C}\tanh \left(\tfrac{{Ex}\left({\bar{P}}_{C}\right)-{Ex}\left({\bar{P}}_{D}\right)}{2K}\right).\end{array}\right.\end{eqnarray} \tag{ A.4 }$$

Defining ${p}_{X}^{{\prime} }={p}_{X}/\left(1-{p}_{\phi }\right)$ (X ∈ E, C, D), we obtain

$$\begin{eqnarray}&&\left\{\begin{array}{l}\tfrac{{{dp}}_{E}^{{\prime} }}{{dt}}=\left(1-{p}_{\phi }\right)\left[-{p}_{E}^{{\prime} }{p}_{C}^{{\prime} }\tanh \left(\tfrac{{Ex}\left({\bar{P}}_{C}\right)-{Ex}\left({\bar{P}}_{E}\right)}{2K}\right)-{p}_{E}^{{\prime} }{p}_{D}^{{\prime} }\tanh \left(\tfrac{{Ex}\left({\bar{P}}_{D}\right)-{Ex}\left({\bar{P}}_{E}\right)}{2K}\right)\right],\\ \tfrac{{{dp}}_{C}^{{\prime} }}{{dt}}=\left(1-{p}_{\phi }\right)\left[-{p}_{C}^{{\prime} }{p}_{E}^{{\prime} }\tanh \left(\tfrac{{Ex}\left({\bar{P}}_{E}\right)-{Ex}\left({\bar{P}}_{C}\right)}{2K}\right)-{p}_{C}^{{\prime} }{p}_{D}^{{\prime} }\tanh \left(\tfrac{{Ex}\left({\bar{P}}_{D}\right)-{Ex}\left({\bar{P}}_{C}\right)}{2K}\right)\right],\\ \tfrac{{{dp}}_{D}^{{\prime} }}{{dt}}=\left(1-{p}_{\phi }\right)\left[-{p}_{D}^{{\prime} }{p}_{E}^{{\prime} }\tanh \left(\tfrac{{Ex}\left({\bar{P}}_{E}\right)-{Ex}\left({\bar{P}}_{D}\right)}{2K}\right)-{p}_{D}^{{\prime} }{p}_{C}^{{\prime} }\tanh \left(\tfrac{{Ex}\left({\bar{P}}_{C}\right)-{Ex}\left({\bar{P}}_{D}\right)}{2K}\right)\right],\end{array}\right.\end{eqnarray} \tag{ A.5 }$$

which is the master equations governing the evolution of the whole system.

In this appendix, we would like to present an extended pair approximation method for describing pairwise interaction dynamics of prisoner's dilemma with expellers on regular graphs. Let p_E, p_C, p_D and p_ϕ denote the density of expellers, cooperators, defectors and vacant sites in a population, respectively. Let p_EE, p_EC, p_CE, p_CC, p_ED, p_DE, p_DD, p_Eϕ, p_ϕE, p_ϕϕ, p_CD, p_DC, p_Cϕ, p_ϕC, p_Dϕ and p_ϕD represent the density of EE, EC, CE, CC, ED, DE, DD, Eϕ, ϕE, ϕϕ, CD, DC, Cϕ, ϕC, Dϕ and ϕD pairs, respectively. Then ${q}_{X| Y}={p}_{{XY}}/{p}_{Y}$ specifies the conditional probability that the neighboring site of a site of state Y is in state X. Herein, X and Y stand for the state of a site, which is occupied by either an expeller or a cooperator or a defector, or just vacant. For pair approximation method [15, 108–112], only frequencies of state pairs p_XY are tracked. The probabilities of larger configurations are expressed and approximated by the frequencies of pair configurations. Based on the symmetry condition p_XY = p_YX, the compatibility condition ${p}_{X}={\sum }_{Y}{p}_{{XY}}$, and the closure conditions, the whole system can be described by the following eight variables in pair approximation: p_E, p_C, p_EE, p_EC, p_ED, p_CC, p_CD and p_DD.

During the pairwise interaction process, the configuration frequencies do not change except if a defector is ejected by an expeller after they play games with each other. Let us now consider the case that an ED pair is randomly selected for interaction from the population. The expulsion of the chosen expeller to the chosen defector from its current site to another vacant site in the population may lead to the variation of p_ED, p_CD and p_DD. The expelled defector has k_E expellers, k_C cooperators, k_D defectors and k_ϕ = k − 1 − k_E − k_C − k_D vacant sites among the k − 1 remaining neighbors on a regular graph with connectivity k. The frequency of such a configuration is

$$\begin{eqnarray*}&&\begin{array}{l}\tfrac{\left(k-1\right)!}{{k}_{\phi }!\left(k-1-{k}_{\phi }\right)!}{{q}_{\phi | {DE}}}^{{k}_{\phi }}{\left(1-{q}_{\phi | {DE}}\right)}^{k-1-{k}_{\phi }}\tfrac{\left(k-1-{k}_{\phi }\right)!}{{k}_{E}!\left(k-1-{k}_{\phi }-{k}_{E}\right)!}{\left(\tfrac{{q}_{E| {DE}}}{1-{q}_{\phi | {DE}}}\right)}^{{k}_{E}}{\left(1-\tfrac{{q}_{E| {DE}}}{1-{q}_{\phi | {DE}}}\right)}^{k-1-{k}_{\phi }-{k}_{E}}\tfrac{\left({k}_{C}+{k}_{D}\right)!}{{k}_{C}!{k}_{D}!}\\ \quad \times \,{\left(\tfrac{{q}_{C| {DE}}}{{q}_{C| {DE}}+{q}_{D| {DE}}}\right)}^{{k}_{C}}{\left(\tfrac{{q}_{D| {DE}}}{{q}_{C| {DE}}+{q}_{D| {DE}}}\right)}^{{k}_{D}}.\end{array}\end{eqnarray*}$$

The frequency of the configuration, in which a vacant site has ${k}_{E}^{\prime}$ expellers, ${k}_{C}^{\prime}$ cooperators, ${k}_{D}^{\prime}$ defectors and ${k}_{\phi }^{\prime} =k-{k}_{E}^{\prime} -{k}_{C}^{\prime} -{k}_{D}^{\prime}$ vacant sites, is

$$\begin{eqnarray*}&&\begin{array}{l}\tfrac{k!}{{k}_{\phi }^{\prime} !\left(k-{k}_{\phi }^{\prime} \right)!}{{q}_{\phi | \phi }}^{{k}_{\phi }^{\prime} }{\left(1-{q}_{\phi | \phi }\right)}^{k-{k}_{\phi }^{\prime} }\tfrac{\left(k-{k}_{\phi }^{\prime} \right)!}{{k}_{E}^{\prime} !\left(k-{k}_{\phi }^{\prime} -{k}_{E}^{\prime} \right)!}{\left(\tfrac{{q}_{E| \phi }}{1-{q}_{\phi | \phi }}\right)}^{{k}_{E}^{\prime} }{\left(1-\tfrac{{q}_{E| \phi }}{1-{q}_{\phi | \phi }}\right)}^{k-{k}_{\phi }^{\prime} -{k}_{E}^{\prime} }\tfrac{\left({k}_{C}^{\prime} +{k}_{D}^{\prime} \right)!}{{k}_{C}^{\prime} !{k}_{D}^{\prime} !}{\left(\tfrac{{q}_{C| \phi }}{{q}_{C| \phi }+{q}_{D| \phi }}\right)}^{{k}_{C}^{\prime} }\\ \quad \times \,{\left(\tfrac{{q}_{D| \phi }}{{q}_{C| \phi }+{q}_{D| \phi }}\right)}^{{k}_{D}^{\prime} }.\end{array}\end{eqnarray*}$$

Therefore, the number of ED pairs increases by ${k}_{E}^{\prime} -{k}_{E}-1$, and thus p_ED increases by $2\left({k}_{E}^{\prime} -{k}_{E}-1\right)/\left({kN}\right)$ with probability

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{Prob}\left({\rm{\Delta }}{p}_{{ED}}=\displaystyle \frac{2\left({k}_{E}^{\prime} -{k}_{E}-1\right)}{{kN}}\right) & = & {p}_{{ED}}\displaystyle \frac{k!}{{k}_{\phi }^{\prime} !\left(k-{k}_{\phi }^{\prime} \right)!}{{q}_{\phi | \phi }}^{{k}_{\phi }^{\prime} }{\left(1-{q}_{\phi | \phi }\right)}^{k-{k}_{\phi }^{\prime} }\displaystyle \frac{\left(k-{k}_{\phi }^{\prime} \right)!}{{k}_{E}^{\prime} !\left(k-{k}_{\phi }^{\prime} -{k}_{E}^{\prime} \right)!}\\ & & \times {\left(\displaystyle \frac{{q}_{E| \phi }}{1-{q}_{\phi | \phi }}\right)}^{{k}_{E}^{\prime} }{\left(1-\displaystyle \frac{{q}_{E| \phi }}{1-{q}_{\phi | \phi }}\right)}^{k-{k}_{\phi }^{\prime} -{k}_{E}^{\prime} }\displaystyle \frac{\left({k}_{C}^{\prime} +{k}_{D}^{\prime} \right)!}{{k}_{C}^{\prime} !{k}_{D}^{\prime} !}\\ & & \times {\left(\displaystyle \frac{{q}_{C| \phi }}{{q}_{C| \phi }+{q}_{D| \phi }}\right)}^{{k}_{C}^{\prime} }{\left(\displaystyle \frac{{q}_{D| \phi }}{{q}_{C| \phi }+{q}_{D| \phi }}\right)}^{{k}_{D}^{\prime} }\displaystyle \frac{\left(k-1\right)!}{{k}_{\phi }!\left(k-1-{k}_{\phi }\right)!}{{q}_{\phi | {DE}}}^{{k}_{\phi }}\\ & & \times {\left(1-{q}_{\phi | {DE}}\right)}^{k-1-{k}_{\phi }}\displaystyle \frac{\left(k-1-{k}_{\phi }\right)!}{{k}_{E}!\left(k-1-{k}_{\phi }-{k}_{E}\right)!}{\left(\displaystyle \frac{{q}_{E| {DE}}}{1-{q}_{\phi | {DE}}}\right)}^{{k}_{E}}\\ & & \times {\left(1-\displaystyle \frac{{q}_{E| {DE}}}{1-{q}_{\phi | {DE}}}\right)}^{k-1-{k}_{\phi }-{k}_{E}}\displaystyle \frac{\left({k}_{C}+{k}_{D}\right)!}{{k}_{C}!{k}_{D}!}{\left(\displaystyle \frac{{q}_{C| {DE}}}{{q}_{C| {DE}}+{q}_{D| {DE}}}\right)}^{{k}_{C}}\\ & & \times {\left(\displaystyle \frac{{q}_{D| {DE}}}{{q}_{C| {DE}}+{q}_{D| {DE}}}\right)}^{{k}_{D}}.\end{array}\end{eqnarray} \tag{ B.1 }$$

The number of CD pairs increases by ${k}_{C}^{\prime} -{k}_{C}$, and thus p_CD increases by $2\left({k}_{C}^{\prime} -{k}_{C}\right)/\left({kN}\right)$ with probability

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{Prob}\left({\rm{\Delta }}{p}_{{CD}}=\displaystyle \frac{2\left({k}_{C}^{\prime} -{k}_{C}\right)}{{kN}}\right) & = & {p}_{{ED}}\displaystyle \frac{k!}{{k}_{\phi }^{\prime} !\left(k-{k}_{\phi }^{\prime} \right)!}{{q}_{\phi | \phi }}^{{k}_{\phi }^{\prime} }{\left(1-{q}_{\phi | \phi }\right)}^{k-{k}_{\phi }^{\prime} }\displaystyle \frac{\left(k-{k}_{\phi }^{\prime} \right)!}{{k}_{E}^{\prime} !\left(k-{k}_{\phi }^{\prime} -{k}_{E}^{\prime} \right)!}\\ & & \times {\left(\displaystyle \frac{{q}_{E| \phi }}{1-{q}_{\phi | \phi }}\right)}^{{k}_{E}^{\prime} }{\left(1-\displaystyle \frac{{q}_{E| \phi }}{1-{q}_{\phi | \phi }}\right)}^{k-{k}_{\phi }^{\prime} -{k}_{E}^{\prime} }\displaystyle \frac{\left({k}_{C}^{\prime} +{k}_{D}^{\prime} \right)!}{{k}_{C}^{\prime} !{k}_{D}^{\prime} !}\\ & & \times {\left(\displaystyle \frac{{q}_{C| \phi }}{{q}_{C| \phi }+{q}_{D| \phi }}\right)}^{{k}_{C}^{\prime} }{\left(\displaystyle \frac{{q}_{D| \phi }}{{q}_{C| \phi }+{q}_{D| \phi }}\right)}^{{k}_{D}^{\prime} }\displaystyle \frac{\left(k-1\right)!}{{k}_{\phi }!\left(k-1-{k}_{\phi }\right)!}{{q}_{\phi | {DE}}}^{{k}_{\phi }}\\ & & \times {\left(1-{q}_{\phi | {DE}}\right)}^{k-1-{k}_{\phi }}\displaystyle \frac{\left(k-1-{k}_{\phi }\right)!}{{k}_{E}!\left(k-1-{k}_{\phi }-{k}_{E}\right)!}{\left(\displaystyle \frac{{q}_{E| {DE}}}{1-{q}_{\phi | {DE}}}\right)}^{{k}_{E}}\\ & & \times {\left(1-\displaystyle \frac{{q}_{E| {DE}}}{1-{q}_{\phi | {DE}}}\right)}^{k-1-{k}_{\phi }-{k}_{E}}\displaystyle \frac{\left({k}_{C}+{k}_{D}\right)!}{{k}_{C}!{k}_{D}!}{\left(\displaystyle \frac{{q}_{C| {DE}}}{{q}_{C| {DE}}+{q}_{D| {DE}}}\right)}^{{k}_{C}}\\ & & \times {\left(\displaystyle \frac{{q}_{D| {DE}}}{{q}_{C| {DE}}+{q}_{D| {DE}}}\right)}^{{k}_{D}}.\end{array}\end{eqnarray} \tag{ B.2 }$$

Similarly, the number of DD pairs increases by ${k}_{D}^{\prime} -{k}_{D}$, and thus p_DD increases by $2\left({k}_{D}^{\prime} -{k}_{D}\right)/\left({kN}\right)$ with probability

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{Prob}\left({\rm{\Delta }}{p}_{{DD}}=\displaystyle \frac{2\left({k}_{D}^{\prime} -{k}_{D}\right)}{{kN}}\right) & = & {p}_{{ED}}\displaystyle \frac{k!}{{k}_{\phi }^{\prime} !\left(k-{k}_{\phi }^{\prime} \right)!}{{q}_{\phi | \phi }}^{{k}_{\phi }^{\prime} }{\left(1-{q}_{\phi | \phi }\right)}^{k-{k}_{\phi }^{\prime} }\displaystyle \frac{\left(k-{k}_{\phi }^{\prime} \right)!}{{k}_{E}^{\prime} !\left(k-{k}_{\phi }^{\prime} -{k}_{E}^{\prime} \right)!}\\ & & \times {\left(\displaystyle \frac{{q}_{E| \phi }}{1-{q}_{\phi | \phi }}\right)}^{{k}_{E}^{\prime} }{\left(1-\displaystyle \frac{{q}_{E| \phi }}{1-{q}_{\phi | \phi }}\right)}^{k-{k}_{\phi }^{\prime} -{k}_{E}^{\prime} }\displaystyle \frac{\left({k}_{C}^{\prime} +{k}_{D}^{\prime} \right)!}{{k}_{C}^{\prime} !{k}_{D}^{\prime} !}\\ & & \times {\left(\displaystyle \frac{{q}_{C| \phi }}{{q}_{C| \phi }+{q}_{D| \phi }}\right)}^{{k}_{C}^{\prime} }{\left(\displaystyle \frac{{q}_{D| \phi }}{{q}_{C| \phi }+{q}_{D| \phi }}\right)}^{{k}_{D}^{\prime} }\displaystyle \frac{\left(k-1\right)!}{{k}_{\phi }!\left(k-1-{k}_{\phi }\right)!}{{q}_{\phi | {DE}}}^{{k}_{\phi }}\\ & & \times {\left(1-{q}_{\phi | {DE}}\right)}^{k-1-{k}_{\phi }}\displaystyle \frac{\left(k-1-{k}_{\phi }\right)!}{{k}_{E}!\left(k-1-{k}_{\phi }-{k}_{E}\right)!}{\left(\displaystyle \frac{{q}_{E| {DE}}}{1-{q}_{\phi | {DE}}}\right)}^{{k}_{E}}\\ & & \times {\left(1-\displaystyle \frac{{q}_{E| {DE}}}{1-{q}_{\phi | {DE}}}\right)}^{k-1-{k}_{\phi }-{k}_{E}}\displaystyle \frac{\left({k}_{C}+{k}_{D}\right)!}{{k}_{C}!{k}_{D}!}{\left(\displaystyle \frac{{q}_{C| {DE}}}{{q}_{C| {DE}}+{q}_{D| {DE}}}\right)}^{{k}_{C}}\\ & & \times {\left(\displaystyle \frac{{q}_{D| {DE}}}{{q}_{C| {DE}}+{q}_{D| {DE}}}\right)}^{{k}_{D}}.\end{array}\end{eqnarray} \tag{ B.3 }$$

In the limit of large site sizes $N\to \infty$, we obtain the following differential equation set:

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\dot{p}}_{E}=\mathop{\mathrm{lim}}\limits_{N\to \infty }\tfrac{{\rm{\Delta }}{p}_{E}}{1/N}=0,\\ {\dot{p}}_{C}=\mathop{\mathrm{lim}}\limits_{N\to \infty }\tfrac{{\rm{\Delta }}{p}_{C}}{1/N}=0,\\ {\dot{p}}_{{EE}}=\mathop{\mathrm{lim}}\limits_{N\to \infty }\tfrac{{\rm{\Delta }}{p}_{{EE}}}{1/N}=0,\\ {\dot{p}}_{{EC}}=\mathop{\mathrm{lim}}\limits_{N\to \infty }\tfrac{{\rm{\Delta }}{p}_{{EC}}}{1/N}=0,\\ {\dot{p}}_{{ED}}=\mathop{\mathrm{lim}}\limits_{N\to \infty }\tfrac{{\rm{\Delta }}{p}_{{ED}}}{1/N}=\tfrac{2\left({k}_{E}^{\prime} -{k}_{E}-1\right)}{k}\displaystyle \sum _{{k}_{E}=0}^{k-1}\displaystyle \sum _{{k}_{E}^{\prime} =0}^{k}\mathrm{Prob}\left({\rm{\Delta }}{p}_{{ED}}=\tfrac{2\left({k}_{E}^{\prime} -{k}_{E}-1\right)}{{kN}}\right),\\ {\dot{p}}_{{CC}}=\mathop{\mathrm{lim}}\limits_{N\to \infty }\tfrac{{\rm{\Delta }}{p}_{{CC}}}{1/N}=0,\\ {\dot{p}}_{{CD}}=\mathop{\mathrm{lim}}\limits_{N\to \infty }\tfrac{{\rm{\Delta }}{p}_{{CD}}}{1/N}=\tfrac{2\left({k}_{C}^{\prime} -{k}_{C}\right)}{k}\displaystyle \sum _{{k}_{C}}^{k-1}\displaystyle \sum _{{k}_{C}^{\prime} =0}^{k}\mathrm{Prob}\left({\rm{\Delta }}{p}_{{CD}}=\tfrac{2\left({k}_{C}^{\prime} -{k}_{C}\right)}{{kN}}\right),\\ {\dot{p}}_{{DD}}=\mathop{\mathrm{lim}}\limits_{N\to \infty }\tfrac{{\rm{\Delta }}{p}_{{DD}}}{1/N}=\tfrac{2\left({k}_{D}^{\prime} -{k}_{D}\right)}{k}\displaystyle \sum _{{k}_{D}=0}^{k-1}\displaystyle \sum _{{k}_{D}^{\prime} =0}^{k}\mathrm{Prob}\left({\rm{\Delta }}{p}_{{DD}}=\tfrac{2\left({k}_{D}^{\prime} -{k}_{D}\right)}{{kN}}\right).\end{array}\right.\end{eqnarray} \tag{ B.4 }$$

After some algebra, equation (B.4) can be simplified to

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\dot{p}}_{E}=0,\\ {\dot{p}}_{C}=0,\\ {\dot{p}}_{{EE}}=0,\\ {\dot{p}}_{{EC}}=0,\\ {\dot{p}}_{{ED}}=\tfrac{2{p}_{{ED}}}{k}\left[{{kq}}_{E| \phi }-\left(k-1\right){q}_{E| {DE}}-1\right],\\ {\dot{p}}_{{CC}}=0,\\ {\dot{p}}_{{CD}}=\tfrac{2{p}_{{ED}}}{k}\left[{{kq}}_{C| \phi }-\left(k-1\right){q}_{C| {DE}}\right],\\ {\dot{p}}_{{DD}}=\tfrac{2{p}_{{ED}}}{k}\left[{{kq}}_{D| \phi }-\left(k-1\right){q}_{D| {DE}}\right].\end{array}\right.\end{eqnarray} \tag{ B.5 }$$

The above equation set requires a 'moment closure' by approximating ${q}_{X| {YZ}}\approx {q}_{X| Y}$. This means that only first order pair correlations are taken into account and hence the name 'pair approximation'. Thus from equation (B.5), we have

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\dot{p}}_{E}=0,\\ {\dot{p}}_{C}=0,\\ {\dot{p}}_{{EE}}=0,\\ {\dot{p}}_{{EC}}=0,\\ {\dot{p}}_{{ED}}=\tfrac{2{p}_{{ED}}}{k}\left[{{kq}}_{E| \phi }-\left(k-1\right){q}_{E| D}-1\right],\\ {\dot{p}}_{{CC}}=0,\\ {\dot{p}}_{{CD}}=\tfrac{2{p}_{{ED}}}{k}\left[{{kq}}_{C| \phi }-\left(k-1\right){q}_{C| D}\right],\\ {\dot{p}}_{{DD}}=\tfrac{2{p}_{{ED}}}{k}\left[{{kq}}_{D| \phi }-\left(k-1\right){q}_{D| D}\right].\end{array}\right.\end{eqnarray} \tag{ B.6 }$$

Further rearranging equation (B.6), we obtain

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\dot{p}}_{E}=0,\\ {\dot{p}}_{C}=0,\\ {\dot{p}}_{{EE}}=0,\\ {\dot{p}}_{{EC}}=0,\\ {\dot{p}}_{{ED}}=\left[-\tfrac{2}{k}\left(\tfrac{k}{{p}_{\phi }}+\tfrac{k-1}{{p}_{D}}\right)\right]{{p}_{{ED}}}^{2}+\left[\tfrac{2}{k}\left(k\tfrac{{p}_{E}-{p}_{{EE}}-{p}_{{EC}}}{{p}_{\phi }}-1\right)\right]{p}_{{ED}},\\ {\dot{p}}_{{CC}}=0,\\ {\dot{p}}_{{CD}}=\left[-\left(\tfrac{k}{{p}_{\phi }}+\tfrac{k-1}{{p}_{D}}\right)\tfrac{2{p}_{{ED}}}{k}\right]{p}_{{CD}}+\left[\tfrac{2\left({p}_{C}-{p}_{{EC}}-{p}_{{CC}}\right){p}_{{ED}}}{{p}_{\phi }}\right],\\ {\dot{p}}_{{DD}}=\left[-\left(\tfrac{k}{{p}_{\phi }}+\tfrac{k-1}{{p}_{D}}\right)\tfrac{2{p}_{{ED}}}{k}\right]{p}_{{DD}}+\left[\tfrac{2\left({p}_{D}-{p}_{{ED}}-{p}_{{CD}}\right){p}_{{ED}}}{{p}_{\phi }}\right].\end{array}\right.\end{eqnarray} \tag{ B.7 }$$

From equation (B.7), we can see that the ordinary differential equation of ${\dot{p}}_{{ED}}$ has the form of Bernoulli's equation. Dividing throughout by ${{p}_{{ED}}}^{2}$, one gets

$$\begin{eqnarray}&&\displaystyle \frac{d{{p}_{{ED}}}^{-1}}{{dt}}+\left[\displaystyle \frac{2}{k}\left(k\displaystyle \frac{{p}_{E}-{p}_{{EE}}-{p}_{{EC}}}{{p}_{\phi }}-1\right)\right]{{p}_{{ED}}}^{-1}=\displaystyle \frac{2}{k}\left(\displaystyle \frac{k}{{p}_{\phi }}+\displaystyle \frac{k-1}{{p}_{D}}\right).\end{eqnarray} \tag{ B.8 }$$

Multiplying throughout by ${e}^{\int \left[\tfrac{2}{k}\left(k\tfrac{{p}_{E}-{p}_{{EE}}-{p}_{{EC}}}{{p}_{\phi }}-1\right)\right]{dt}}$, equation (B.8) is transformed into

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\left({{p}_{{ED}}}^{-1}{e}^{\int \left[\tfrac{2}{k}\left(k\displaystyle \frac{{p}_{E}-{p}_{{EE}}-{p}_{{EC}}}{{p}_{\phi }}-1\right)\right]{dt}}\right)=\displaystyle \frac{2}{k}\left(\displaystyle \frac{k}{{p}_{\phi }}+\displaystyle \frac{k-1}{{p}_{D}}\right){e}^{\int \left[\tfrac{2}{k}\left(k\displaystyle \frac{{p}_{E}-{p}_{{EE}}-{p}_{{EC}}}{{p}_{\phi }}-1\right)\right]{dt}},\end{eqnarray} \tag{ B.9 }$$

which leads to a general solution

$$\begin{eqnarray}&&{{p}_{{ED}}}^{-1}{e}^{\int \left[\tfrac{2}{k}\left(k\displaystyle \frac{{p}_{E}-{p}_{{EE}}-{p}_{{EC}}}{{p}_{\phi }}-1\right)\right]{dt}}=\int \displaystyle \frac{2}{k}\left(\displaystyle \frac{k}{{p}_{\phi }}+\displaystyle \frac{k-1}{{p}_{D}}\right){e}^{\int \left[\tfrac{2}{k}\left(k\displaystyle \frac{{p}_{E}-{p}_{{EE}}-{p}_{{EC}}}{{p}_{\phi }}-1\right)\right]{dt}}{dt}+C,\end{eqnarray} \tag{ B.10 }$$

where C is a constant. Multiplying ${e}^{\int -\left[\tfrac{2}{k}\left(k\tfrac{{p}_{E}-{p}_{{EE}}-{p}_{{EC}}}{{p}_{\phi }}-1\right)\right]{dt}}$ solves for ${{p}_{{ED}}}^{-1}$:

$$\begin{eqnarray}&&{{p}_{{ED}}}^{-1}={e}^{\int -\left[\tfrac{2}{k}\left(k\displaystyle \frac{{p}_{E}-{p}_{{EE}}-{p}_{{EC}}}{{p}_{\phi }}-1\right)\right]{dt}}\left[\int \displaystyle \frac{2}{k}\left(\displaystyle \frac{k}{{p}_{\phi }}+\displaystyle \frac{k-1}{{p}_{D}}\right){e}^{\int \left[\tfrac{2}{k}\left(k\displaystyle \frac{{p}_{E}-{p}_{{EE}}-{p}_{{EC}}}{{p}_{\phi }}-1\right)\right]{dt}}{dt}+C\right].\end{eqnarray} \tag{ B.11 }$$

Therefore, we obtain

$$\begin{eqnarray}&&{p}_{{ED}}=\displaystyle \frac{{e}^{\int \left[\tfrac{2}{k}\left(k\tfrac{{p}_{E}-{p}_{{EE}}-{p}_{{EC}}}{{p}_{\phi }}-1\right)\right]{dt}}}{\int \tfrac{2}{k}\left(\tfrac{k}{{p}_{\phi }}+\tfrac{k-1}{{p}_{D}}\right){e}^{\int \left[\tfrac{2}{k}\left(k\tfrac{{p}_{E}-{p}_{{EE}}-{p}_{{EC}}}{{p}_{\phi }}-1\right)\right]{dt}}{dt}+C}.\end{eqnarray} \tag{ B.12 }$$

Assuming that at the beginning of each pairwise interaction stage, expellers with a fraction ${p}_{E}\left(t^{\prime} =0\right)$, cooperators with a fraction ${p}_{C}\left(t^{\prime} =0\right)$ and defectors with a fraction ${p}_{D}\left(t^{\prime} =0\right)$ are distributed on a regular graph with the vacancy rate ${p}_{\phi }=1-\left({p}_{E}+{p}_{C}+{p}_{D}\right)$. Thus ${p}_{{ED}}\left(t^{\prime} \right)$ is given by:

$$\begin{eqnarray}&&{p}_{{ED}}\left(t^{\prime} \right)=\displaystyle \frac{{e}^{\tfrac{2}{k}\left(k\tfrac{{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)}{{p}_{\phi }}-1\right)t^{\prime} }}{\tfrac{{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }}{k\left[{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)\right]{p}_{D}\left(0\right)-{p}_{\phi }{p}_{D}\left(0\right)}\left({e}^{\tfrac{2}{k}\left(k\tfrac{{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)}{{p}_{\phi }}-1\right)t^{\prime} }-1\right)+\tfrac{1}{{p}_{{ED}}\left(0\right)}}.\end{eqnarray} \tag{ B.13 }$$

Similarly, ${p}_{{CD}}\left(t^{\prime} \right)$ and ${p}_{{DD}}\left(t^{\prime} \right)$ can be respectively expressed by

$$\begin{eqnarray}&&{p}_{{CD}}\left(t^{\prime} \right)=\displaystyle \frac{\tfrac{k\left[{p}_{C}\left(0\right)-{p}_{{EC}}\left(0\right)-{p}_{{CC}}\left(0\right)\right]}{k\left[{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)\right]-{p}_{\phi }}\left({e}^{\tfrac{2}{k}\left(k\tfrac{{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)}{{p}_{\phi }}-1\right)t^{\prime} }-1\right)+\tfrac{{p}_{{CD}}\left(0\right)}{{p}_{{ED}}\left(0\right)}}{\tfrac{{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }}{k\left[{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)\right]{p}_{D}\left(0\right)-{p}_{\phi }{p}_{D}\left(0\right)}\left({e}^{\tfrac{2}{k}\left(k\tfrac{{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)}{{p}_{\phi }}-1\right)t^{\prime} }-1\right)+\tfrac{1}{{p}_{{ED}}\left(0\right)}},\end{eqnarray} \tag{ B.14 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{p}_{{DD}}\left(t^{\prime} \right) & = & \displaystyle \frac{\tfrac{{k}^{2}{p}_{D}\left(0\right)\left\{{p}_{D}\left(0\right)+{p}_{\phi }-\left[{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)\right]-\left[{p}_{C}\left(0\right)-{p}_{{EC}}\left(0\right)-{p}_{{CC}}\left(0\right)\right]\right\}}{\left\{k\left[{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)\right]-{p}_{\phi }\right\}\left[{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }\right]}\left({e}^{\tfrac{2}{k}\left(k\tfrac{{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)}{{p}_{\phi }}-1\right)t^{\prime} }-1\right)+\tfrac{{p}_{{DD}}\left(0\right)}{{p}_{{ED}}\left(0\right)}}{\tfrac{{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }}{k\left[{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)\right]{p}_{D}\left(0\right)-{p}_{\phi }{p}_{D}\left(0\right)}\left({e}^{\tfrac{2}{k}\left(k\tfrac{{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)}{{p}_{\phi }}-1\right)t^{\prime} }-1\right)+\tfrac{1}{{p}_{{ED}}\left(0\right)}}\\ & & +\displaystyle \frac{\left\{\tfrac{k\left[{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)\right]{p}_{D}\left(0\right)-{p}_{\phi }{p}_{D}\left(0\right)}{{p}_{{ED}}\left(0\right)\left[{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }\right]}\left(1+\tfrac{k\left[{p}_{C}\left(0\right)-{p}_{{EC}}\left(0\right)-{p}_{{CC}}\left(0\right)\right]}{k\left[{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)\right]-{p}_{\phi }}\right)-1-\tfrac{{p}_{{CD}}\left(0\right)}{{p}_{{ED}}\left(0\right)}\right\}\tfrac{{{kp}}_{D}\left(0\right)}{{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }}}{\tfrac{{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }}{k\left[{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)\right]{p}_{D}\left(0\right)-{p}_{\phi }{p}_{D}\left(0\right)}\left({e}^{\tfrac{2}{k}\left(k\tfrac{{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)}{{p}_{\phi }}-1\right)t^{\prime} }-1\right)+\tfrac{1}{{p}_{{ED}}\left(0\right)}}\\ & & \times \mathrm{ln}\left\{\displaystyle \frac{{p}_{{ED}}\left(0\right)\left[{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }\right]}{k\left[{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)\right]{p}_{D}\left(0\right)-{p}_{\phi }{p}_{D}\left(0\right)}\left({e}^{\tfrac{2}{k}\left(k\tfrac{{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)}{{p}_{\phi }}-1\right)t^{\prime} }-1\right)+1\right\}.\end{array}\end{eqnarray} \tag{ B.15 }$$

For rest of equations in equation (B.7), the current densities of players and the current pair densities can be simply solved as ${p}_{X}\left(t^{\prime} \right)={p}_{X}\left(0\right)\left(X\in E,C\right)$ and ${p}_{{XY}}\left(t^{\prime} \right)={p}_{{XY}}\left(0\right)\left(X,Y\in E,C\right)$, respectively.

From equations (B.13)–(B.15) we can observe that there exists a critical fraction of vacant sites

$$\begin{eqnarray}&&{{p}_{\phi }}^{c}=k\left[{p}_{E}\left(0\right)-{p}_{{EE}}\left(0\right)-{p}_{{EC}}\left(0\right)\right],\end{eqnarray} \tag{ B.16 }$$

depending on which the equilibrium state of the population at the end of pairwise interaction stage can be summarized into the following three generic cases:

• (i)  
  If ${p}_{\phi }\lt {{p}_{\phi }}^{c}$, the stationary state is given by

  $$\begin{eqnarray}&&\left\{\begin{array}{l}{p}_{E}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{E}\left(0\right),\\ {p}_{C}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{C}\left(0\right),\\ {p}_{{EE}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{{EE}}\left(0\right),\\ {p}_{{EC}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{{EC}}\left(0\right),\\ {p}_{{ED}}\left(t^{\prime} ={t}_{e}^{\prime} \right)=\tfrac{{{kp}}_{D}\left(0\right)\left[{p}_{{ED}}\left(0\right)+{p}_{E\phi }\left(0\right)\right]-{p}_{\phi }{p}_{D}\left(0\right)}{{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }},\\ {p}_{{CC}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{{CC}}\left(0\right),\\ {p}_{{CD}}\left(t^{\prime} ={t}_{e}^{\prime} \right)=\tfrac{{{kp}}_{D}\left(0\right)\left[{p}_{{CD}}\left(0\right)+{p}_{C\phi }\left(0\right)\right]}{{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }},\\ {p}_{{DD}}\left(t^{\prime} ={t}_{e}^{\prime} \right)=\tfrac{{k}^{2}{p}_{D}{\left(0\right)}^{2}\left\{{p}_{D}\left(0\right)+{p}_{\phi }-\left[{p}_{{ED}}\left(0\right)+{p}_{E\phi }\left(0\right)\right]-\left[{p}_{{CD}}\left(0\right)+{p}_{C\phi }\left(0\right)\right]\right\}}{{\left[{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }\right]}^{2}};\end{array}\right.\end{eqnarray} \tag{ B.17 }$$
• (ii)  
  If ${p}_{\phi }={{p}_{\phi }}^{c}$, the stationary state is given by

  $$\begin{eqnarray}&&\left\{\begin{array}{l}{p}_{E}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{E}\left(0\right),\\ {p}_{C}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{C}\left(0\right),\\ {p}_{{EE}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{{EE}}\left(0\right),\\ {p}_{{EC}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{{EC}}\left(0\right),\\ {p}_{{ED}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{{ED}}\left(0\right),\\ {p}_{{CC}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{{CC}}\left(0\right),\\ {p}_{{CD}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{{CD}}\left(0\right),\\ {p}_{{DD}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{{DD}}\left(0\right);\end{array}\right.\end{eqnarray} \tag{ B.18 }$$
• (iii)  
  If ${p}_{\phi }\gt {{p}_{\phi }}^{c}$, the stationary state is given by

  $$\begin{eqnarray}&&\left\{\begin{array}{l}{p}_{E}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{E}\left(0\right),\\ {p}_{C}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{C}\left(0\right),\\ {p}_{{EE}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{{EE}}\left(0\right),\\ {p}_{{EC}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{{EC}}\left(0\right),\\ {p}_{{ED}}\left(t^{\prime} ={t}_{e}^{\prime} \right)=0,\\ {p}_{{CC}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={p}_{{CC}}\left(0\right),\\ {p}_{{CD}}\left(t^{\prime} ={t}_{e}^{\prime} \right)=\tfrac{{p}_{D}\left(0\right){p}_{{CD}}\left(0\right)\left\{k\left[{p}_{{ED}}\left(0\right)+{p}_{E\phi }\left(0\right)\right]-{p}_{\phi }\right\}-{{kp}}_{D}\left(0\right){p}_{{ED}}\left(0\right)\left[{p}_{{CD}}\left(0\right)+{p}_{C\phi }\left(0\right)\right]}{{p}_{D}\left(0\right)\left\{k\left[{p}_{{ED}}\left(0\right)+{p}_{E\phi }\left(0\right)\right]-{p}_{\phi }\right\}-{p}_{{ED}}\left(0\right)\left[{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }\right]},\\ {p}_{{DD}}\left(t^{\prime} ={t}_{e}^{\prime} \right)=\tfrac{{p}_{D}\left(0\right){p}_{{DD}}\left(0\right)\left\{k\left[{p}_{{ED}}\left(0\right)+{p}_{E\phi }\left(0\right)\right]-{p}_{\phi }\right\}}{{p}_{D}\left(0\right)\left\{k\left[{p}_{{ED}}\left(0\right)+{p}_{E\phi }\left(0\right)\right]-{p}_{\phi }\right\}-{p}_{{ED}}\left(0\right)\left[{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }\right]}\\ \quad -\,\tfrac{{k}^{2}{p}_{D}{\left(0\right)}^{2}{p}_{{ED}}\left(0\right)\left\{{p}_{D}\left(0\right)+{p}_{\phi }-\left[{p}_{{ED}}\left(0\right)+{p}_{E\phi }\left(0\right)\right]-\left[{p}_{{CD}}\left(0\right)+{p}_{C\phi }\left(0\right)\right]\right\}}{{p}_{D}\left(0\right)\left\{k\left[{p}_{{ED}}\left(0\right)+{p}_{E\phi }\left(0\right)\right]-{p}_{\phi }\right\}\left[{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }\right]-{p}_{{ED}}\left(0\right){\left[{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }\right]}^{2}}\\ \quad +\,\tfrac{{p}_{D}\left(0\right)\left\{k\left[{p}_{{ED}}\left(0\right)+{p}_{E\phi }\left(0\right)\right]-{p}_{\phi }\right\}-\left[{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }\right]\left[{p}_{{ED}}\left(0\right)+{p}_{{CD}}\left(0\right)\right]+{{kp}}_{D}\left(0\right)\left[{p}_{{CD}}\left(0\right)+{p}_{C\phi }\left(0\right)\right]}{{p}_{D}\left(0\right)\left\{k\left[{p}_{{ED}}\left(0\right)+{p}_{E\phi }\left(0\right)\right]-{p}_{\phi }\right\}{\left[{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }\right]}^{2}-{p}_{{ED}}\left(0\right){\left[{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }\right]}^{3}}\\ \quad \times \,{{kp}}_{D}{\left(0\right)}^{2}\left\{k\left[{p}_{{ED}}\left(0\right)+{p}_{E\phi }\left(0\right)\right]-{p}_{\phi }\right\}\mathrm{ln}\left(1-\tfrac{{p}_{{ED}}\left(0\right)\left[{{kp}}_{D}\left(0\right)+\left(k-1\right){p}_{\phi }\right]}{{p}_{D}\left(0\right)\left\{k\left[{p}_{{ED}}\left(0\right)+{p}_{E\phi }\left(0\right)\right]-{p}_{\phi }\right\}}\right),\end{array}\right.\end{eqnarray} \tag{ B.19 }$$

where ${t}_{e}^{\prime}$ represents the ending time of pairwise interaction stage. Obviously, as w → 0, the time interval ${t}_{e}^{\prime}$ of pairwise interaction stage becomes much longer than the time that is required for the whole system enters into the equilibrium state in the pairwise interaction stage.

To demonstrate the predictive power of the extended pair approximation method for describing pairwise interaction dynamics of spatial prisoner's dilemma with expellers, we now consider a typical class of initial conditions: at the start of pairwise interaction process, a fraction of p_ϕ random sites of a regular graph with degree k are set to be vacant, and the remaining sites are occupied by either expellers, cooperators or defectors with equal probability. In this case, the initial values of local densities can be estimated as ${q}_{X| Y}(0)\approx {p}_{X}(0)=(1-{p}_{\phi })/3$ (X ∈ E, C, D and Y ∈ E, C, D, ϕ) and ${q}_{\phi | Y}(0)\approx {p}_{\phi }$ (Y ∈ E, C, D, ϕ). Thus the solutions of ${p}_{{ED}}\left(t^{\prime} \right)$, ${p}_{{CD}}\left(t^{\prime} \right)$ and ${p}_{{DD}}\left(t^{\prime} \right)$ becomes

$$\begin{eqnarray}&&{p}_{{ED}}\left(t^{\prime} \right)=\displaystyle \frac{\left(1-{p}_{\phi }\right){e}^{2\left(\tfrac{1+{p}_{\phi }-2{{p}_{\phi }}^{2}}{9{p}_{\phi }}-\tfrac{1}{k}\right)t^{\prime} }}{\tfrac{9\left[k+\left(2k-3\right){p}_{\phi }\right]}{k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}}\left({e}^{2\left(\tfrac{1+{p}_{\phi }-2{{p}_{\phi }}^{2}}{9{p}_{\phi }}-\tfrac{1}{k}\right)t^{\prime} }-1\right)+\tfrac{9}{1-{p}_{\phi }}},\end{eqnarray} \tag{ B.20 }$$

$$\begin{eqnarray}&&{p}_{{CD}}\left(t^{\prime} \right)=\displaystyle \frac{\tfrac{k\left(1+{p}_{\phi }-2{{p}_{\phi }}^{2}\right)}{k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}}\left({e}^{2\left(\tfrac{1+{p}_{\phi }-2{{p}_{\phi }}^{2}}{9{p}_{\phi }}-\tfrac{1}{k}\right)t^{\prime} }-1\right)+1}{\tfrac{9\left[k+\left(2k-3\right){p}_{\phi }\right]}{\left(1-{p}_{\phi }\right)\left[k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}\right]}\left({e}^{2\left(\tfrac{1+{p}_{\phi }-2{{p}_{\phi }}^{2}}{9{p}_{\phi }}-\tfrac{1}{k}\right)t^{\prime} }-1\right)+\tfrac{9}{{\left(1-{p}_{\phi }\right)}^{2}}},\end{eqnarray} \tag{ B.21 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{p}_{{DD}}\left(t^{\prime} \right) & = & \displaystyle \frac{\tfrac{{k}^{2}\left(1-{p}_{\phi }\right){\left(1+2{p}_{\phi }\right)}^{2}}{\left[k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}\right]\left[k+\left(2k-3\right){p}_{\phi }\right]}\left({e}^{2\left(\tfrac{1+{p}_{\phi }-2{{p}_{\phi }}^{2}}{9{p}_{\phi }}-\tfrac{1}{k}\right)t^{\prime} }-1\right)+1}{\tfrac{9\left[k+\left(2k-3\right){p}_{\phi }\right]}{\left(1-{p}_{\phi }\right)\left[k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}\right]}\left({e}^{2\left(\tfrac{1+{p}_{\phi }-2{{p}_{\phi }}^{2}}{9{p}_{\phi }}-\tfrac{1}{k}\right)t^{\prime} }-1\right)+\tfrac{9}{{\left(1-{p}_{\phi }\right)}^{2}}}\\ & & -\displaystyle \frac{\tfrac{3{{kp}}_{\phi }\left(1+2{p}_{\phi }\right)}{{\left[k+\left(2k-3\right){p}_{\phi }\right]}^{2}}\mathrm{ln}\left\{\tfrac{\left(1-{p}_{\phi }\right)\left[k+\left(2k-3\right){p}_{\phi }\right]}{k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}}\left({e}^{2\left(\tfrac{1+{p}_{\phi }-2{{p}_{\phi }}^{2}}{9{p}_{\phi }}-\tfrac{1}{k}\right)t^{\prime} }-1\right)+1\right\}}{\tfrac{9\left[k+\left(2k-3\right){p}_{\phi }\right]}{\left(1-{p}_{\phi }\right)\left[k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}\right]}\left({e}^{2\left(\tfrac{1+{p}_{\phi }-2{{p}_{\phi }}^{2}}{9{p}_{\phi }}-\tfrac{1}{k}\right)t^{\prime} }-1\right)+\tfrac{9}{{\left(1-{p}_{\phi }\right)}^{2}}},\end{array}\end{eqnarray} \tag{ B.22 }$$

respectively. While the current densities of players and the current pair densities can be respectively solved as ${p}_{X}\left(t^{\prime} \right)=\left(1-{p}_{\phi }\right)/3\left(X\in E,C\right)$ and ${p}_{{XY}}\left(t^{\prime} \right)={\left(1-{p}_{\phi }\right)}^{2}/9\left(X,Y\in E,C\right)$.

According to equations (B.20)–(B.22) the critical ratio of vacant sites can be given as

$$\begin{eqnarray}&&{{p}_{\phi }}^{c}=\displaystyle \frac{k-9+3\sqrt{{k}^{2}-2k+9}}{4k}\end{eqnarray} \tag{ B.23 }$$

by solving $\left(1+{p}_{\phi }-2{{p}_{\phi }}^{2}\right)/\left(9{p}_{\phi }\right)=1/k$. Still, the stationary state of the system at the end of pairwise interaction process can be classified into the following three generic cases:

• (i)  
  If ${p}_{\phi }\lt {{p}_{\phi }}^{c}$, the equilibrium state becomes

  $$\begin{eqnarray}&&\left\{\begin{array}{l}{p}_{E}\left(t^{\prime} ={t}_{e}^{\prime} \right)=\left(1-{p}_{\phi }\right)/3,\\ {p}_{C}\left(t^{\prime} ={t}_{e}^{\prime} \right)=\left(1-{p}_{\phi }\right)/3,\\ {p}_{{EE}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={\left(1-{p}_{\phi }\right)}^{2}/9,\\ {p}_{{EC}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={\left(1-{p}_{\phi }\right)}^{2}/9,\\ {p}_{{ED}}\left(t^{\prime} ={t}_{e}^{\prime} \right)=\tfrac{\left(1-{p}_{\phi }\right)\left[k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}\right]}{9\left[k+\left(2k-3\right){p}_{\phi }\right]},\\ {p}_{{CC}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={\left(1-{p}_{\phi }\right)}^{2}/9,\\ {p}_{{CD}}\left(t^{\prime} ={t}_{e}^{\prime} \right)=\tfrac{k{\left(1-{p}_{\phi }\right)}^{2}\left(1+2{p}_{\phi }\right)}{9\left[k+\left(2k-3\right){p}_{\phi }\right]},\\ {p}_{{DD}}\left(t^{\prime} ={t}_{e}^{\prime} \right)=\tfrac{{k}^{2}{\left(1-{p}_{\phi }\right)}^{2}{\left(1+2{p}_{\phi }\right)}^{2}}{9{\left[k+\left(2k-3\right){p}_{\phi }\right]}^{2}};\end{array}\right.\end{eqnarray} \tag{ B.24 }$$
• (ii)  
  If ${p}_{\phi }={{p}_{\phi }}^{c}$, the equilibrium state becomes

  $$\begin{eqnarray}&&\left\{\begin{array}{l}{p}_{E}\left(t^{\prime} ={t}_{e}^{\prime} \right)=\left(1-{p}_{\phi }\right)/3,\\ {p}_{C}\left(t^{\prime} ={t}_{e}^{\prime} \right)=\left(1-{p}_{\phi }\right)/3,\\ {p}_{{EE}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={\left(1-{p}_{\phi }\right)}^{2}/9,\\ {p}_{{EC}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={\left(1-{p}_{\phi }\right)}^{2}/9,\\ {p}_{{ED}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={\left(1-{p}_{\phi }\right)}^{2}/9,\\ {p}_{{CC}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={\left(1-{p}_{\phi }\right)}^{2}/9,\\ {p}_{{CD}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={\left(1-{p}_{\phi }\right)}^{2}/9,\\ {p}_{{DD}}\left(t^{\prime} ={t}_{e}^{\prime} \right)={\left(1-{p}_{\phi }\right)}^{2}/9;\end{array}\right.\end{eqnarray} \tag{ B.25 }$$
• (iii)  
  If ${p}_{\phi }\gt {{p}_{\phi }}^{c}$, the equilibrium state becomes

  $$\begin{eqnarray}\left\{\begin{array}{rcl}{p}_{E}\left(t^{\prime} ={t}_{e}^{\prime} \right) & = & \left(1-{p}_{\phi }\right)/3,\\ {p}_{C}\left(t^{\prime} ={t}_{e}^{\prime} \right) & = & \left(1-{p}_{\phi }\right)/3,\\ {p}_{{EE}}\left(t^{\prime} ={t}_{e}^{\prime} \right) & = & {\left(1-{p}_{\phi }\right)}^{2}/9,\\ {p}_{{EC}}\left(t^{\prime} ={t}_{e}^{\prime} \right) & = & {\left(1-{p}_{\phi }\right)}^{2}/9,\\ {p}_{{ED}}\left(t^{\prime} ={t}_{e}^{\prime} \right) & = & 0,\\ {p}_{{CC}}\left(t^{\prime} ={t}_{e}^{\prime} \right) & = & {\left(1-{p}_{\phi }\right)}^{2}/9,\\ {p}_{{CD}}\left(t^{\prime} ={t}_{e}^{\prime} \right) & = & \tfrac{{\left(1-{p}_{\phi }\right)}^{2}}{3\left(2+{p}_{\phi }\right)},\\ {p}_{{DD}}\left(t^{\prime} ={t}_{e}^{\prime} \right) & = & \tfrac{{\left(1-{p}_{\phi }\right)}^{2}\left[k+\left(2k-3\right){p}_{\phi }\right]\left\{\left[k+\left(2k-3\right){p}_{\phi }\right]\left[k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}\right]-{k}^{2}\left(1-{p}_{\phi }\right){\left(1+2{p}_{\phi }\right)}^{2}\right\}}{9{\left[k+\left(2k-3\right){p}_{\phi }\right]}^{2}\left\{\left[k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}\right]-\left(1-{p}_{\phi }\right)\left[k+\left(2k-3\right){p}_{\phi }\right]\right\}}\\ & & -\tfrac{{{kp}}_{\phi }{\left(1-{p}_{\phi }\right)}^{2}\left(1+2{p}_{\phi }\right)\left[k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}\right]}{3{\left[k+\left(2k-3\right){p}_{\phi }\right]}^{2}\left\{\left[k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}\right]-\left(1-{p}_{\phi }\right)\left[k+\left(2k-3\right){p}_{\phi }\right]\right\}}\\ & & \times \mathrm{ln}\left\{1-\tfrac{\left(1-{p}_{\phi }\right)\left[k+\left(2k-3\right){p}_{\phi }\right]}{k+\left(k-9\right){p}_{\phi }-2k{{p}_{\phi }}^{2}}\right\}.\end{array}\right.\end{eqnarray} \tag{ B.26 }$$

Figures B1 and B2 show that simulation and analytical results almost perfectly coincide with each other both in steady state and during evolutionary process of pairwise interaction dynamics, and thus confirm the validity of our method. Besides, from equations (B.20) and (B.21) we can obtain

$$\begin{eqnarray}&&\left\{\begin{array}{l}\tfrac{{{dp}}_{{ED}}\left(t^{\prime} \right)}{{dt}^{\prime} }=0\ \mathrm{if}\ {p}_{\phi }={{p}_{\phi }}^{c},\\ \tfrac{{{dp}}_{{ED}}\left(t^{\prime} \right)}{{dt}^{\prime} }\lt 0\ \mathrm{if}\ {p}_{\phi }\ne {{p}_{\phi }}^{c},\end{array}\right.\end{eqnarray} \tag{ B.27 }$$

and

$$\begin{eqnarray}&&\left\{\begin{array}{l}\tfrac{{{dp}}_{{CD}}\left(t^{\prime} \right)}{{dt}^{\prime} }=0\ \mathrm{if}\ {p}_{\phi }={{p}_{\phi }}^{c},\\ \tfrac{{{dp}}_{{CD}}\left(t^{\prime} \right)}{{dt}^{\prime} }\gt 0\ \mathrm{if}\ {p}_{\phi }\ne {{p}_{\phi }}^{c}.\end{array}\right.\end{eqnarray} \tag{ B.28 }$$

Here, the overall impacts of pairwise interaction dynamics on strategy distribution are to decrease interaction frequency between expellers and defectors while increase interaction frequency between cooperators and defectors if only ${p}_{\phi }\ne {{p}_{\phi }}^{c}$ (figures B1 and B2), which reveals the origin of evolutionary advantages of expellers over cooperators as well as defectors. For ${p}_{\phi }={{p}_{\phi }}^{c}$, the current local densities become invariant with time step $t^{\prime}$ during the whole pairwise interaction process.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size