Altruism in populations at the extinction transition

TL;DR

This paper models the evolution of cooperation in spatial populations, revealing how cooperation can persist at extinction thresholds when benefits outweigh costs, with phase transitions characterized by spatial correlations.

Contribution

It introduces a birth-death process model showing how cooperation survives at extinction transitions, emphasizing the role of benefit-cost ratios and spatial correlations.

Findings

Full cooperation occurs at the extinction transition when benefit exceeds cost.

Spatial correlations are sufficient for sustained cooperation.

Phase diagrams are similar on lattices and in pair approximation.

Abstract

We study the evolution of cooperation as a birth-death process in spatially extended populations. The benefit from the altruistic behavior of a cooperator is implemented by decreasing the death rate of its direct neighbors. The cost of cooperation is the increase of a cooperator's death rate proportional to the number of its neighbors. When cooperation has higher cost than benefit, cooperators disappear. Then the dynamics reduces to the contact process with defectors as the single particle type. Increasing the benefit-cost ratio above 1, the extinction transition of the contact process splits into a set of nonequilibrium transitions between four regimes when increasing the baseline death rate $p$ as a control parameter: (i) defection only, (ii) coexistence, (iii) cooperation only, (iv) extinction. We investigate the transitions between these regimes. As the main result, we find that full cooperation is established at the extinction transition as long as benefit is strictly larger than cost. Qualitatively identical phase diagrams are obtained for populations on square lattices and in pair approximation. Spatial correlations with nearest neighbors only are thus sufficient for sustained cooperation.

Equations

$$S=\{c_i,d_i{\}}_{i=1}^N\equiv \{x_i{\}}_{i=1}^N,x_i\in \{c_i,d_i\},$$

$$e_i=1\stackrel{{\pi }_b(c_i,S)}{\to }{c}_i=1$$

$$e_i=1\stackrel{{\pi }_b(d_i,S)}{\to }{d}_i=1$$

$${\pi }_b(c_i,S)=e_i\sum _{j\in N_i}{c}_j/k_j\equiv e_i{\tilde{c}}_i,$$

$${\pi }_b(d_i,S)=e_i\sum _{j\in N_i}{d}_j/k_j\equiv e_i{\tilde{d}}_i,$$

$$c_i=1\stackrel{{\pi }_d(c_i,S)}{\to }{e}_i=1,$$

$$d_i=1\stackrel{{\pi }_d(d_i,S)}{\to }{e}_i=1,$$

$${\pi }_d(c_i,S)=$$

$$=$$

$${\pi }_d(d_i,S)=$$

$$=$$

$$\frac{\partial (⟨c⟩+⟨d⟩)}{\partial p}>0.$$

$$p(x)=\{\begin{array}{ll}\frac{2x-1}{L}& \text{if }x\le L/2\\ \frac{2(L-x)+1}{L}& \text{otherwise.}\end{array}$$

$$\begin{array}{rl}{\partial }_{t}P(S,t)=\sum _{i=1}^N& \sum _{x_i\in \{c_i,d_i\}}\{(E_{x_i}^{-}-1)\left[{\pi }_b(x_i,S)P(S,t)\right]\\ & +(E_{x_i}^{+}-1)\left[{\pi }_d(x_i,S)P(S,t)\right]\},\end{array}$$

$$\frac{d}{dt}⟨c_i⟩=$$

$$=$$

$$\frac{d}{dt}⟨d_i⟩=$$

$$=$$

$$\frac{d}{dt}⟨c_i{c}_j⟩=$$

$$=$$

$$\frac{d}{dt}⟨c_i{d}_j⟩=$$

$$=$$

$$\frac{d}{dt}⟨d_i{d}_j⟩=$$

$$=$$

$$⟨dd⟩=(1-p)⟨d⟩,$$

$$p\left[1-k(1-p)\right]⟨d⟩+(k-1)⟨ded⟩=0.$$

$$p\le 1-\frac{1}{k}.$$

$$⟨cc⟩=\frac{1-p}{1-ϵp}⟨c⟩,$$

$$\begin{array}{rl}& \frac{p}{1-ϵp}\left[p(1-ϵ)-(k-1)(1-p)\right]⟨c⟩\\ & +(k-1)\left(⟨cec⟩+pϵ⟨ccc⟩\right)=0,\end{array}$$

$$p\le 1-\frac{1-ϵ}{k-ϵ}\ge 1-\frac{1}{k}.$$

$$⟨x_i{x}_j⟩\simeq ⟨x_i⟩⟨x_j⟩,i\ne j$$

$$⟨x_i⟩\simeq ⟨x_j⟩\equiv ⟨x⟩,\text{for all }i.$$

$$\frac{d}{dt}⟨c⟩=$$

$$\frac{d}{dt}⟨d⟩=$$

$$⟨c⟩=0\mathrm{and}⟨d⟩=1-p,$$

$$⟨c⟩=\frac{1-p}{1-ϵp}\mathrm{and}⟨d⟩=0.$$

$$\frac{d}{dt}{⟨c⟩}_1\simeq -p(1-p)(1-ϵ){⟨c⟩}_1,$$

Full text

I Introduction

Altruism or cooperativity tr71 describe behavior that is more in favor of others than of the actor herself. Alarm calls are an example of altruistic behavior: Increasing the risk of becoming prey itself first, one individual of a group warns the others of a predator approaching ta87 . At first glance, the observation of altruism sustained over generations appears incompatible with Darwin’s theory of natural selection, featuring the survival of the fittest darwin1859 ; spencer1864 . If non-altruists acting only to their own benefit have an advantage over altruists in terms of reproductive success, altruistic traits eventually disappear.

The question of sustained altruism and cooperativity has been addressed in the framework of evolutionary game theory, in particular by work on the Prisoner’s Dilemma and Public Goods Games da80 ; axha81 . In these and other games, the time evolution of the system is assumed to take place as a sequence of two elementary steps: (i) the combined behavioral choices of the participants lead to an assignment of a payoff to each player which (ii) determines the choice of their future strategies or roles. In the simplest case, with two possible strategies, cooperation and defection, the dilemma arises as follows. Regardless of the other agent’s move, an agent’s best (highest payoff) move is always defection. On the other hand, the sum of all players payoffs is maximal when all cooperate. Therefore, natural selection always favors defection da80 , despite cooperation is the best global strategy.

The aforementioned social dilemma is frequently analyzed by means of the replicator equation ohno06 ; rocusa09 ; brpa18 describing the time evolution of the fraction of players holding one of the two strategies. If the fitness of an individual equals its payoff, the resulting replicator equation for the Prisoner’s Dilemma has only two steady-state solutions, the only-defector and the only-cooperator solutions, the former being the only stable one. Nevertheless, the prevalence of cooperation is still possible within the context of evolutionary games, provided appropriate reciprocity mechanisms are included in the dynamics tr71 ; axha81 ; leke06 ; za14 ; wakojuta15 : Direct reciprocity, indirect reciprocity, kin selection, group selection, and network structure. If compared to the well-mixed situation, the new mechanisms include update rules that favor the interactions among cooperators.

The network structure mechanism was one of the first reciprocity mechanisms studied in the literature. It refers to the restriction of agents interactions among neighbors. In a two-dimensional regular network, the survival of altruists was explained in terms of their ability of preventing the exploitation of defectors through the formation of clusters noma92 ; namaiw97 ; tamaiw97 ; iwnale98 . Further progress in the field considered births and deaths: The second step of the dynamics, the one that allows a change of the strategy, is now interpreted as a death of a player followed by a birth. The new ecologic perspective allowed to assess the importance of new relevant issues, such as the fluctuation of the population density miwi00 ; lilihiiwboni13 ; Huang:2015 ; mcfrhawano18 , the movement of agents bara98 ; yiliwalulimali14 ; busz05 ; vasiar07 ; lero10 , the spatial distribution of neighbors and their number ifkido04 ; sena09 , among others. Recent works also consider networks of interactions ohno06 ; asgola08 ; fuwanoha09 ; notaan10 ; pegoszflmo13 ; wakojuta15 ; kiyoki16 , focus on the critical properties of the system hasz05 ; szvusz05 ; vuszsz06 , include other novel dynamic rules raamnaoh10 ; pesz10 ; pisapa12 ; liliclgu15 ; pachadno15 ; iyki16 ; szpe17 , analyze the formation of patterns noma92 ; szfa07 ; lanoha08 ; funoha10 ; waha11 ; yabasa18 , and evaluate the effect on the population growing as external pressure rises Sella:2000 . The latter aspect has been widely analyzed in the context of competing species gamere13 ; dodi13 , but has not received much attention in relation with the prevalence of altruism.

Although general considerations about the prevalence of altruism in the context of the Public Goods Games can be inferred from the numerous studies on the topic ohhalino06 ; ko11 , the behavior of cooperation turns out to be very dependent on the specific dynamics considered ha06 ; rocusa09a . This is the case when trying to evaluate the importance of the spatial heterogeneity and the formation of clusters of cooperators: Many studies namaiw97 ; iwnale98 ; lefedi03 ; thel03 ; lilihiiwboni13 explain the coexistence of cooperation and defection using the so called pair approximation, an approach that goes one step beyond mean field by tracking the dynamics of pairs of neighbors. However, pair approximation still assumes spatial homogeneity of the system. Hence, there is no need for the formation of clusters of cooperators for explaining their long-term survival.

Recent works on the evolution of cooperation suggest the need of giving up on certain common statements of evolutionary game theory gaferutacusamo12 ; grrosemitr12 ; pejorawabosz17 ; sa18 . Particularly, some experiments on the dynamics of human cooperation show that people choose their strategy regardless the payoff of others grgrmisetrcumosa14 . Similar conclusions are given in the context of living beings doissi17 . See also recent experimental and numerical works on related topics waszpe12 ; kudi13 ; nadrfo16 ; leahha17 .

Here we study the evolution of cooperation in the framework of interacting particle systems. We model birth and death in a spatially extended population as a contact process and ask the following: What is the phase diagram of the contact process with an additional — cooperative— type of particle that supports survival at neighboring sites? Our approach provides a natural framework to assess the effects of different mechanisms on the behavior of the system and on the survival of cooperativity, such as the dynamics of interactions, the fluctuation in the population size, the presence or absence of cooperation clusters, and the spatial variation of parameters, among others.

The organization of the work is as follows. In Sec. II we introduce the agent-based model of a population of cooperators and defectors living on a generic network. For later sections the main focus is on the square lattice, where the system has only three relevant parameters: the total number of sites $N$, the death parameter $p$, and the cost-of-altruism parameter $\epsilon$. Section III includes stochastic simulations. We obtain the phase diagram in the parameter space $(p,\epsilon)$ showing the steady-state configurations of the system. The effect of $p$ being spatially dependent is also addressed. In Sec. IV the system is described theoretically. Three complementary formulations, using main-field or pair-approximation approaches, are given. They aim at describing the system under different physical conditions. Finally, a discussion and outlook of the main results are included in Sec. V.

II Definition of the model

The model describes the evolution of a population on an arbitrary network with $N$ nodes. The set of neighbors of a node $i$ is denoted by $N_{i}$. The network is symmetric (undirected), so that $j\in N_{i}$ implies $i\in N_{j}$; also $i\notin N_{i}$ (no self-loops). Each agent in the population is either a cooperator $C$ or a defector $D$, with $c_{i}$ and $d_{i}$ being their respective numbers at node $i$. A site or node of the network holds at most one agent ($C$ or $D$) but it may also be empty ($E$), hence $0\leq c_{i}+d_{i}\leq 1$ and $c_{i}d_{i}=0$. Thus the state of the system $S$ is given by

[TABLE]

where $X$ is either a cooperator or a defector, and $x_{i}$ its number at site $i$. Moreover, the number $e_{i}=1-x_{i}$ gives $1$ if site $i$ is empty and [math] if $x_{i}=1$. From condition $c_{i}d_{i}=0$ we also have $e_{i}x_{i}=0$.

A state transition is either the birth or the death of one agent at a site $i$. At the birth of a cooperator we set $c_{i}=1$ at a previously empty site $i$,

[TABLE]

Likewise for the birth of a defector, $d_{i}=1$ is set at an empty site $i$,

[TABLE]

These transitions occur at a rate proportional to the fraction of neighboring sites occupied by the agent type to be born, as

[TABLE]

where $k_{i}=|N_{i}|$ is the degree (number of those neighbors) of node $i$, and $\tilde{x}_{i}\equiv\sum_{j\in N_{i}}x_{j}/k_{j}$. The death of an agent is a state transition setting $c_{i}=0$ or $d_{i}=0$ at a previously occupied site $i$,

[TABLE]

with respective rates

[TABLE]

where now $\bar{x}_{i}\equiv k_{i}^{-1}\sum_{j\in N_{i}}x_{j}$. Agents die at a baseline rate $p$. This rate is reduced, however, by the fraction of adjacent sites occupied by a cooperator. The death rate of a cooperator, on the other hand, has an additional positive term proportional (with factor $1-\epsilon$) to the fraction of adjacent agents. This way, the parameter $\epsilon$ accounts for the cost of the altruistic act, the limit of $\epsilon=0$ corresponding to maximum cost where the altruist definitely loses its life for saving that of its neighbor. The other limit is costless altruism at $\epsilon=1$.

In the absence of cooperators, or in the absence of defectors with $\epsilon=0$, the model reduces to the contact process ha74 ; madi05 equivalent to the SIS (susceptible-infected-susceptible) model of epidemics he89 ; albrdrwu08 . The equivalence is obtained by mapping each empty site to a susceptible individual and each site with a defector to an infected individual.

III Simulations

Let us first illustrate and numerically analyze the dynamics on periodic square lattices. As defined above, the model features non-ergodicity. Eventually both types of agents go extinct in a finite size system. In the simulations in this section, a slightly modified version of the model is employed: We set to zero the death rate of an agent currently being the only one of its type (C or D). This allows us to take long-term measurements of concentrations and distributions without having to restart the dynamics. Given the rates, simulations are performed with a standard Gillespie algorithm gi76 ; gi07 .

III.1 Square lattice with homogeneous parameters

Figure 1 shows the parameter dependence of the stationary mean concentrations of agents. At $\epsilon=0$, cooperators are absent in the whole range of $p$, while the concentration of defectors is positive for $p<p_{c}\approx 0.62$ and vanishes for $p>p_{c}$. Now fixing $0<\epsilon<1$ and increasing $p$ from [math] to $1$, the concentration of defectors $\left\langle d\right\rangle$ still decreases with $p$. Before $\left\langle d\right\rangle$ reaches zero, however, the concentration of cooperators $\left\langle c\right\rangle$ becomes positive. Simulations on square lattices of smaller size ($N=20^{2}$, $N=30^{2}$) and checks with $N=100^{2}$ yield results almost identical to those of Fig. 1.

In the coexistence regime of cooperators and defectors (green area in Figure 1), the growth of cooperation outweighs the decline of defection. Here the total concentration of agents grows with $p$,

[TABLE]

Figure 2(a) explicitly shows this non-monotonicity by plotting $\left\langle c\right\rangle+\left\langle d\right\rangle$ versus $p$ for different choices of $\epsilon$.

Let us now take a closer look at the transitions between the regimes observed in Figure 1. To this end, we record the distributions in the number of agents (each type separately) and consider their changes under parameter variation. Figure 3 shows this analysis for five transitions (a)-(e), also marked in the bottom panel of Figure 1. Transitions in Figs. 3(a) and 3(e) are extinctions of one type of agent in the absence of the other type. However, the transitions are distinguishable by the approximate exponents of the algebraic decay of distributions, giving $1/4$ for the extinction of defectors versus $3/7$ for cooperators. This indicates that, even in the absence of defectors and close to the extinction transition (e), the dynamics of cooperators is essentially different from the contact process.

Differences in the distributions of the order parameter (Figure 1), however, do not contradict transitions (a)-(e) belonging to the same universality class. Transitions (a), (b) and (e) fulfill the premises of the directed percolation conjecture, cf. section 3.3.6 in hi00 . Transitions (c) and (d) do not fulfill the assumption of a unique absorbing state because only one type of agent goes extinct at the transition. In preliminary numerical exploration (results not shown here), we have found the scaling of the order parameter (concentration of agents) compatible with the value $0.580(4)$ for exponent $\beta$ in directed percolation in two dimensions wazhligade13 . We conjecture that all transitions (a)-(e) belong to the universality class of directed percolation.

III.2 Spatially dependent parameter $p$

Let us study a variation of the model with a spatial dependence of the parameter $p$, a way of mimicking ecological conditions ch00 ; okle13 . For an agent at lattice site $(x,y)$, $x,y\in\{1,\dots,L\}$, the death rate is based on the parameter value

[TABLE]

For $L$ even, the minimum value $1/L$ is assumed by $p(x)$ at $x=1$ and $x=L$; its maximum value $1-1/L$ is obtained at $x=L/2$ and $x=L/2+1$. The parameter $\epsilon$ remains spatially homogeneous, here $\epsilon=0.75$.

Figure 4(a) shows the concentration of agents as a function of lattice coordinate $x$, i.e. averaged over lattice coordinate $y$ and time. We see that the effect of parameter $p$ is local. The $p$-dependence of $\left\langle c\right\rangle$ and $\left\langle d\right\rangle$ observed under spatially homogeneous $p$ in Section III.1 qualitatively matches that of the scenario with spatially dependent $p$.

IV Analytic approximations

In this section, we derive three complementary theoretical descriptions of our model, defined in Sec. II. The first two ones are based on a mean-field approximation, while the third one uses the pair approximation. As will be shown, the different approaches have different ranges of applicability and explain the prevalence/extinction and even the coexistence of altruism and defection under different physical and biological conditions. In the case of the pair approximation, a very similar phase diagram to the numerical one shown in Fig. 1 is obtained.

Our starting point is the master equation for the probability $P(S,t)$ of finding the system in state $S$ at time $t$. By means of a probabilistic balance in the continuum time limit ka92 , and using the rates given by Eqs. (4)-(9), the master equation reads as

[TABLE]

where the operators $E^{\pm}_{x_{i}}$ act on a generic function $f(x_{1},\dots,x_{i},\dots,x_{N})$ as $E^{\pm}_{x_{i}}f(x_{1},\dots,x_{i},\dots,x_{N})=f(x_{1},\dots,x_{i}\pm 1,\dots,x_{N})$, with $x_{k}\in\{c_{k},d_{k}\}$, $k=1,\dots,N$.

By taking moments of the master equation (12) we can derive equations for the mean numbers of cooperators and defectors in site $i$, $\left\langle c_{i}\right\rangle$ and $\left\langle d_{i}\right\rangle$. After using the relation $e_{i}=1-c_{i}-d_{i}$ and some manipulations, we obtain

[TABLE]

for $i=1,\dots,N$. Since the first moments are coupled to the second ones through correlations between neighbors, it is also convenient to derive equations for the two node correlations for neighboring sites, i.e. $\left\langle x_{i}x_{j}\right\rangle$ with $j\in N_{i}$:

[TABLE]

where $\tilde{x}_{i}$ and $\bar{x}_{i}$ are defined just after Eqs. (5) and (9), respectively. The two remaining moments, $\left\langle c_{i}e_{j}\right\rangle$ and $\left\langle d_{i}e_{j}\right\rangle$ can be obtained from the previous ones by means of the identity $1=e_{i}+c_{i}+d_{i}$, as $\left\langle c_{i}e_{j}\right\rangle=\left\langle c_{i}\right\rangle-\left\langle c_{i}c_{j}\right\rangle-\left\langle c_{i}d_{j}\right\rangle$ and $\left\langle d_{i}e_{j}\right\rangle=\left\langle d_{i}\right\rangle-\left\langle d_{i}d_{j}\right\rangle-\left\langle d_{i}c_{j}\right\rangle$.

Although the system of Eqs. (IV)-(IV) are exact and valid for any structure of neighbors (network), it is not closed, due to the presence of three nodes correlations. Therefore, in order to have a closed set of equations, three approximations are explored. The first two ones make use of the mean-field approximation, where two node correlations are ignored, and the third one uses pair approximation. Furthermore, we restrict ourselves to regular networks where $k_{i}=k$ for all $i$, so as to simplify the description (now $\tilde{x}_{i}=\bar{x}_{i}=k^{-1}\sum_{j\in N_{i}}x_{j}$).

IV.1 Exact relations

Before proceeding with the approximations, some exact relations will be derived. They apply for homogeneous steady-state configurations.

Consider first the case of only defectors. Since $\left\langle c\right\rangle=0$, we also have $\left\langle cc\right\rangle=\left\langle cd\right\rangle=\left\langle ce\right\rangle=0$. Using Eq. (IV), together with $\left\langle de\right\rangle=\left\langle d\right\rangle-\left\langle dd\right\rangle$, we have

[TABLE]

and, with Eq. (IV) and the identity $\left\langle dde\right\rangle+\left\langle ddd\right\rangle=\left\langle dd\right\rangle$,

[TABLE]

which implies, in order to have positive solutions, $1-k(1-p)\leq 0$, that is

[TABLE]

This is an overestimation of the extinction probability of defectors, for all $\epsilon\in[0,1]$. For $\epsilon=0$, where the model is the SIS model, and the square lattice ($k=4$), the previous estimation is $0.75$ while the one from the simulations is around $0.62$ saol02 ; vofama09 , see also Fig. 1.

For the only-cooperator case, it is $\left\langle d\right\rangle=0$ and $\left\langle dd\right\rangle=\left\langle cd\right\rangle=\left\langle de\right\rangle=0$. Using Eq. (IV) together with $\left\langle ce\right\rangle=\left\langle c\right\rangle-\left\langle cc\right\rangle$, we get

[TABLE]

and, with Eq. (IV) and the identity $\left\langle cce\right\rangle+\left\langle ccc\right\rangle=\left\langle cc\right\rangle$,

[TABLE]

which now implies $p(1-\epsilon)-(k-1)(1-p)\leq 0$ or

[TABLE]

Again, this is an overestimation of the critical probability extinction when there are only cooperators in the system. The critical value here is bigger or equal to the one of Eq. (20), as expected due to the altruistic benefit. Equation (23) also provides an estimation of the dependence of the critical probability on $\epsilon$. In particular, it tends to $1$ for $\epsilon\to 1$, in agreement with the numerical simulations of Fig. 1.

Equations (20) and (23), and also the other relations, are the same for $\epsilon=0$ provided we interchange the types of particles, because the model with only defectors and only cooperators coincide in this limit. This can be seen from the rates defining the dynamics in Eqs. (4)-(9): The rates for defectors in the absence of cooperators are the same as the rates for cooperators in the absence of defectors at $\epsilon=0$.

IV.2 Global mean-field approximation

For the global mean-field case, equivalent to the dynamics on a complete graph in the limit of infinite system size, correlations among nodes are absent. In general, assuming the mean-field approximation implies the following two approximations:

[TABLE]

This is also a good approximation when there is no correlation expected between the agents, for instance when there is one kind of agent and the distribution of empty sites is homogeneously distributed. Then, the concentrations $\left\langle c\right\rangle$ and $\left\langle d\right\rangle$ of cooperators and defectors evolve, according to Eqs. (IV) and (IV), as

[TABLE]

The system (26) and (27) can be used now to analyze the homogeneous steady-state solutions. Requiring stationarity, $\frac{d}{dt}\left\langle c\right\rangle=\frac{d}{dt}\left\langle d\right\rangle=0$, we find the trivial solution $\left\langle c\right\rangle=\left\langle d\right\rangle=0$ (all sites empty) and, two other, nontrivial ones, namely

[TABLE]

The trivial solution is clearly unstable, since the coefficient $1-p$ of the less degree terms in Eqs. (26) and (27) is positive for $p<1$. However, it is an absorbing state, and their presence becomes important for small system sizes, as already mentioned in Sec. III.

In order to assess the stability of the solution with only defectors, consider the perturbation of Eq. (28): $\left\langle c\right\rangle=0+\left\langle c\right\rangle_{1}$ and $\left\langle d\right\rangle=1-p+\left\langle d\right\rangle_{1}$ with $\left\langle c\right\rangle_{1}\sim\left\langle d\right\rangle_{1}$. Then, up to linear order in the perturbations, we have

[TABLE]

The first equation, and hence the second one, have $\left\langle c\right\rangle_{1}=\left\langle d\right\rangle_{1}=0$ as the steady solution, revealing the stable character of (28). Proceeding similarly with the only-cooperators solution, Eq. (29), we obtain the system

[TABLE]

which now reveals the unstable character of the solution, since the solution of Eq. (33) increases exponentially with time. According to this analysis, in well-mixed populations, cooperators go extinct.

IV.3 Local mean-field approximation

We can go one step beyond the global mean-field approximation by considering situations where the concentrations of cooperators and defectors change from site to site. In particular, we suppose situations where the site dependence can be encoded through a vector $\mathbf{r}$, which is nothing but the vector of space position in a regular graph. This way, we deduce in the sequel a macroscopic description that removes one of the approximation of the global mean field, namely that of Eq. (25), but still neglects correlations, Eq. (24). The procedure is similar to the one used in Ref. khlohe17 .

By looking at the dynamics on a length scale $L$ much larger than the typical distance between sites $l$, the relevant quantities become local concentrations:

[TABLE]

In a regular graph in $\mathbb{R}^{d}$, for example, $\kappa(\mathbf{r})$ and $\delta(\mathbf{r})$ give the number of cooperators and defectors inside a region of volume $l^{d}$ centered at position $\mathbf{r}$. The new quantities are assumed to be smooth functions of $\mathbf{r}$, a property that allows us to relate any density of site $j\in N_{i}$ and position $\mathbf{l}$, say $\chi(\mathbf{r}+\mathbf{l})=\kappa(\mathbf{r}+\mathbf{l})$ or $\chi(\mathbf{r}+\mathbf{l})=\delta(\mathbf{r}+\mathbf{l})$ at position $\mathbf{r}$, with that of site $i$, $\chi(\mathbf{r})$, as

[TABLE]

Hence, we have

[TABLE]

where we have assumed $\sum_{k\in N_{i}}\mathbf{l}_{k}\simeq 0$, which is an exact expression for a regular square lattice and quiet a good approximation for isotropic configurations. Moreover,

[TABLE]

which is valid, again, under isotropic configurations of sites.

With approximations (24), (37), and (38), the exact system (IV) and (IV) becomes the following reaction-diffusion system

[TABLE]

As expected, we recover the mean-field description for homogeneous solutions, hence we still have the solutions given in Eqs. (28) and (29). However, an important benefit of the present description, if compared to that of the global mean-field approximation, is the possibility of studying the latter solutions under local perturbations, in contrast to homogeneous and global ones done in the previous subsection.

Consider the homogeneous solution of Eq. (28), $\kappa_{0}=0$ and $\delta_{0}=1-p$. Following the standard linear stability analysis, we seek solutions of the form $\kappa=\kappa_{0}+\kappa_{1}$ and $\delta=\delta_{0}+\delta_{1}$, with $\kappa_{1}\sim\delta_{1}\ll\delta_{0}$. After linearizing and seeking solutions of the form $\chi_{1}=\tilde{\chi}_{1}e^{i\boldsymbol{\xi}\cdot\mathbf{r}}$, system (IV.3) and (IV.3) becomes

[TABLE]

The steady state solution for any wavelength $\mathbf{\xi}$ is the trivial one, meaning that the solution of only defectors is linearly stable: any initial and small spatial perturbation in the number of defectors (and also cooperators) tends to zero as time increases.

Proceeding similarly with the solution of Eq. (29), we get

[TABLE]

In this case, the stability of the system depends on the value of $\xi$. Setting $\xi=2\pi/L$, the smallest allowed value for the given boundary conditions, the solution (29) is stable for $p<p_{c}^{*}$ with

[TABLE]

where we have used the approximation $L/l\simeq N^{1/d}$. This means that, under this approximation, the only-cooperators solution is stable for systems small enough. For $N\to\infty$ it is $p_{c}^{*}\to 1$, and the solution is always unstable, and we recover the result of mean field.

Although the local mean-field approximation could in principle be seen as very crude, it shows the importance of taking into account the system size while describing altruism, as already pointed out in Ref. sa18 . In this case, the inclusion of spatial dependence, while still neglecting correlations, stabilizes the only-cooperators solution for $p<\tilde{p}_{c}$. Moreover, the results suggest the existence of other solutions, spatially non-homogeneous ones, and the possibility of discontinuous (first-order) transitions among them. This is because the only-defectors solution keeps always linearly stable, with no other stable solution close to it.

IV.4 Pair approximation

The previous mean-field approaches are expected to fail when the concentration of defectors and cooperators are locally correlated. Since births occur among neighboring sites, correlations are expected to be important, in general. Hence, we reconsider system (IV)-(IV), and try to express the three nodes moments as a function of the one and two nodes mean values. Although different approaches are possible (see for instance khlohe17 ), we explore here the so-called pair approximation. Pair approximation has been extensively applied to a variety of stochastic processes defined on a network, aiming at describing different situations as diverse as spin dynamics olmesa93 ; gl13 , opinion dynamics vaeg08 ; scbe09 ; vacasa10 ; pecasato18 ; pecasato18a , epidemics ledu96 ; eake02 ; tasigrhoki12 , and population dynamics namaiw97 ; iwnale98 ; lilihiiwboni13 ; lefedi03 ; thel03 . In each of the cases, the pair approximation assumes that the probability of a given node quantity $x_{i}$ conditioned to the values of a neighboring site $x_{j}$ and to a next-neighboring site $x_{k}$ is independent of the latter maogsasa92 : $\text{Prob}(x_{i}|x_{j}x_{k})\simeq\text{Prob}(x_{i}|x_{j})$. In other words, the state of a neighbor of a given node is considered to be independent of the state of another neighbor. In our model, where $x_{i}\in\{c,d,e\}$ takes the values [math] or $1$, the mean values $\left\langle x_{i}x_{j}\right\rangle$ and $\left\langle x_{i}x_{j}x_{k}\right\rangle$ are essentially the respective probabilities of the given quantities, hence, under the pair approximation $\left\langle x_{i}x_{j}x_{k}\right\rangle=\text{Prob}(x_{i}|x_{j}x_{k})\text{Prob}(x_{j}x_{k})\simeq P(x_{i}|x_{j})\text{Prob}(x_{j}x_{k})$, we have

[TABLE]

Note that the order of appearance of the variables inside the brackets is important: $x_{i}$ refers to a node which is a neighbor of $x_{j}$ and $x_{j}$ is a neighbor of $x_{k}$. Observe that the pair approximation keeps the correlations regardless of the occupancy of the middle node, namely $\sum_{x_{j}\in\{c,d,e\}}\left\langle x_{i}x_{j}x_{k}\right\rangle=\left\langle x_{i}1x_{k}\right\rangle\simeq\sum_{x_{j}\in\{c,d,e\}}\frac{\left\langle x_{i}x_{j}\right\rangle\left\langle x_{j}x_{k}\right\rangle}{\left\langle x_{j}\right\rangle}\neq\left\langle x_{i}\right\rangle\left\langle x_{j}\right\rangle$, in general.

For simplicity, we consider homogeneous situations for which system (IV)-(IV), within the pair approximation of Eq. (46), becomes

[TABLE]

where $\left\langle xy\right\rangle$ is for any two adjacent nodes with particles $x$ and $y$. Hence, $\left\langle xy\right\rangle=\left\langle yx\right\rangle$.

IV.4.1 Steady-state solutions

The system (47)-(51) has several steady-state solutions. The most obvious one it the trivial solution, without particles, $\left\langle c\right\rangle=\left\langle d\right\rangle=0$. This is the absorbing state we have already mentioned.

By setting all time derivatives of Eqs. (47)-(51) to zero and $c=0$, we obtain the steady-state solution for defection only:

[TABLE]

valid for

[TABLE]

Observe that the previous inequality is the same as the one in Eq. (20), derived using exact relations. However, in this case, $p=1-1/k$ is the exact critical value for the extinction of defectors in the absence of cooperators, within the pair approximation.

The only-cooperators solution is obtained from Eqs. (47)-(51) as a steady-state solution with $d=0$. Now,

[TABLE]

for

[TABLE]

The equality of the last relation holds for $\epsilon\to 0$. For the other limiting value of $\epsilon$, i.e. $\epsilon\to 1$, the upper allowed value of $p$ is $1$, as we also obtained exactly.

Other steady-state solutions describing coexistence, but close to the previous ones, can also be found as follows. First, we notice that the system (47)-(51), under the steady-state condition, can be reduced to a nonlinear system of only two equations with $\left\langle c\right\rangle$ and $\left\langle d\right\rangle$ as unknown quantities. Second, we seek solutions close to the one-type ones, i.e. $\left\langle c\right\rangle\simeq\frac{1-k(1-p)-\epsilon p^{2}}{(1-k+p)(1-\epsilon p)}$ and $\left\langle d\right\rangle\simeq 0$ for the only-cooperators case and $\left\langle c\right\rangle\simeq 0$ and $\left\langle d\right\rangle\simeq\frac{1-k(1-p)}{1-k+p}$ for the only-defectors. For the former case, the resulting equations are linear and nontrivial solutions appear below the following line:

[TABLE]

with

[TABLE]

For the only-defectors case, the resulting set of equation is nonlinear, but one can still find a condition for a nontrivial solution to exist. Now, the nontrivial solutions are above $\epsilon_{d}(p)$, which has the following approximate expression

[TABLE]

with

[TABLE]

Since $\epsilon_{c}(p)\geq\epsilon_{d}(p)$, the coexistence solutions are in the region in between the two lines, as shown in Fig. 5 for the square lattice ($k=4$).

Finally, there may be other solutions describing coexistence not necessarily close to the only-cooperator nor only-defectors ones. This can be shown explicitly for $\epsilon=1$, for which we can obtain explicit expressions. After some algebra, we get

[TABLE]

valid for

[TABLE]

Moreover, it can be seen that this solution is linearly unstable, with only one unstable mode. However, the characteristic time of the unstable mode is much slower than the others, meaning that the system can stay close to the solution for a long time.

IV.4.2 Stability of the steady-state solutions

The stability of the only-defectors and only-cooperators solutions have been studied by means of a modified linear stability analysis of system (47)-(51), following several steps. First, using the identities $\left\langle ce\right\rangle=\left\langle c\right\rangle-\left\langle cc\right\rangle-\left\langle cd\right\rangle$ and $\left\langle de\right\rangle=\left\langle d\right\rangle-\left\langle cd\right\rangle-\left\langle dd\right\rangle$, all mean values are expressed in terms of $\left\langle c\right\rangle$, $\left\langle d\right\rangle$, $\left\langle cc\right\rangle$, $\left\langle cd\right\rangle$, and $\left\langle dd\right\rangle$. Second, the homogeneous solution is linearly perturbed as

[TABLE]

with $\mathbf{u}=(\left\langle c\right\rangle,\left\langle d\right\rangle,\left\langle cc\right\rangle,\left\langle cd\right\rangle,\left\langle dd\right\rangle)$ the vector of the homogeneous solutions, $\mathbf{u}_{0}$ is the vector of the unperturbed solutions, $\mathbf{u}_{1}$ is the perturbation vector, and $\gamma$ a perturbative parameter. Third, the proposed solution is replaced in (47)-(51) and the resulting system is expanded up to linear order in $\gamma$. Contrary to the usual linear perturbation schemes, we obtain a nonlinear closed system of equations for the unknown perturbation quantities $u_{1,i},$ for $i=1,\dots,5$. For both, the only-cooperators and only-defectors solutions, the equation for the perturbation can be written as

[TABLE]

with $M$ being a matrix and $\beta$ a linear function of $\left\langle cc\right\rangle/\left\langle c\right\rangle$, $\left\langle cd\right\rangle/\left\langle c\right\rangle$, $\left\langle cd\right\rangle/\left\langle d\right\rangle$, and $\left\langle dd\right\rangle/\left\langle d\right\rangle$, whose explicit form depends on the solution considered. In any case, $\beta$ is a bounded function, since $0\leq\left\langle xy\right\rangle/\left\langle x\right\rangle\leq 1$ for $x,y\in\{c,d\}$. Finally, the asymptotic behavior of $\mathbf{u}_{1}(t)$ for $t\to\infty$, hence the stable or unstable character of $\mathbf{u}_{0}$, can be determined from the spectra of $M(\beta)$ for any $\beta$, using the following lemma.

Lemma: If all eigenvalues of $M(\beta)$ have negative real parts for all values of $\beta$, then $\mathbf{u}_{0}$ is linearly stable.

Proof: Given a time $t>0$ and an integer $n>0$, we define $t_{i}=\frac{t}{M}i$ for $i=0,\dots,n$. Thanks to the Mean Value Theorem, it is $\mathbf{u}_{1}(t_{i})=\left(I+M_{1}\frac{t}{n}\right)\mathbf{u}(t_{i-1})$ for $i\geq 1$, where $M_{i}$ is the value of $M$ for a time in $(t_{i-1},t_{i})$ and use has been made of Eq. (69) to evaluate the time derivative. By iteration, $\mathbf{u}_{1}(t_{i})=\left[\prod_{k=1}^{i}\left(I+M_{k}\frac{t}{n}\right)\right]\mathbf{u}_{0}$. Denoting by $\lVert\cdot\rVert$ any vector norm, we have $\lVert\mathbf{u}_{1}(t)\rVert=\lVert\left[\prod_{k=1}^{n}(I+M_{k}\frac{t}{n})\right]\mathbf{u}_{0}\rVert\leq\lVert(I+\tilde{M}\frac{t}{n})^{n}\mathbf{u}_{0}\rVert$ where $\tilde{M}$ is such that $\lVert(I+\tilde{M}\frac{t}{n})\mathbf{u}_{0}\rVert=\max_{k}\lVert\left(I+M_{k}\frac{t}{n}\right)\mathbf{u}_{0}\rVert$. Taking $n\to\infty$, $\lVert\mathbf{u}_{1}(t)\rVert\leq\lVert e^{\tilde{M}t}\mathbf{u}_{0}\rVert$ which tends to zero as $t\to\infty$, as all eigenvalues of $\tilde{M}$ have negative real parts. $\square$

Using the lemma, we see that the only-cooperators solution is stable above line $\epsilon_{c}(p)$ given by Eq. (56), and the only-defectors solution is stable below line $\epsilon_{d}(p)$ given approximately by Eq. (60). This implies that the instability of the one-type solutions is due to the presence of coexistence solutions which become stable. Numerical evaluation of the time evolution of system (47)-(51) confirms the theoretical analysis.

V Discussion

As the root of the present work, we introduce a basic stochastic model of a spatially extended population of altrustic and non-altruistic agents, called cooperators and defectors. The population evolves by a birth-death process. In line with the considerations by Huang and colleages Huang:2015 , an agents’ interaction with another agent influences the death rates only. Agents’ interactions are altruistic acts. They lower the death rate of the recipient while increasing the death rate of the donor, relative to a baseline death rate $p$ for agents in isolation, $p$ being one of two parameters of the model. The benefit-cost ratio of the altruistic is encoded in the second parameter $\epsilon$.

Results are obained as (1) stochastic simulations of finite systems and (2) stationary solutions and their stability in approximate descriptions by rate equations. The pair approximation, neglecting all spatial correlations except those of nearest neighbors, yields our main result: For any benefit-cost ratio above $1$, the stable stationary solutions in dependence of baseline death rate $p$ display (i) a regime of co-existence of cooperators and defectors and (ii) a regime of a population of cooperators only. In the $(p,\epsilon)$ parameter plane, these regimes and the related transitions appear as a continuation of the known extinction transition for a spatial population without cooperative interaction (also known as contact process, asynchronous SIS model). The latter case corresponds to benefit-cost ratio of exactly 1 ($\epsilon=0$).

The phase diagram from pair approximation is fully qualitatively consistent with that from stochastic simulation with finite square lattices. Simulations of sufficiently large instances of $k$-regular random graphs yield an equivalent phase diagram (results not shown here); this holds also for preliminary simulation results on other graphs, including scale-free baal99 and small-world networks wast98 . Thus we speculate that the observed type of $(p,\epsilon)$ phase diagram is generic, holding for most types of connected sparse graphs. For dense graphs, however, we expect mean-field behavior without stable cooperation seen in Sec. IV.2.

Consider a spatially extended population subject to a decline in livability, which in reality may be a reduction of food resources or an increase of predators. In our model, this scenario is represented by increasing $p$ and comes with the following prediction. Initially without cooperators, the concentration of agents decreases until reaching a transition point with the onset of co-existence. In this regime, the concentration of defectors further decreases; this decrease is overcompensated by the increase in cooperators. Thus in the co-existence regime, there is net population growth under increasing $p$ Sella:2000 . Further increase of $p$ first leads to a regime with a population containing cooperators only and then an extinction phase where zero population size is the only stable solution.

From earlier studies, both theoretical Sella:2000 and experimental, increasing baseline death rate has been known to enhance cooperation. Perturbing populations of yeast cells by dilution shocks, Sánchez and Gore find populations with larger fractions of cooperative cells (providing digestive enzyme to the population) more likely to survive sago13 . Datta and co-authors observe cooperation promoted when a population expands the space it occupies, cooperators forming a wave of invaders dakocvdugo13 . According to the rule by Ohtsuki and colleagues ohhalino06 , cooperation supersedes defection when the benefit-cost ratio is larger than the agent’s number of neighbors $z$. While their theory assumes a population of constant size and each agent with a constant number $z$ of neighbors, we here see cooperation enhanced when the number of neighbors (occupied adjacent sites) is reduced dynamically due to a shrinking population density. When the population most “needs” it, i.e. at low density close to extinction, cooperation appears as a stable stationary solution for any benefit-cost ratio above 1. Future work may check if a rule relating benefit-cost ratio and neighborhood size characterizes the appearance of stable cooperation also in the present model with varying population size. For experimentally testing the present model’s predictions, the unperturbed steady state of a population would have to be observed.

Giving up the spatial homogeneity of the baseline death rate $p$, we have investigated the scenario of a gradient between low $p$ (high livability) and high $p$ (low livability). The regimes encountered previously by tuning $p$ for the whole system are now found simultaneously at their corresponding spatial position. In particular, high concentration of cooperators is found next to the region uninhabited due to large death rate $p$. There is a region of co-existence where total population concentration increases with $p$ also spatially. Cooperation arises when and where needed to avoid extinction.

Acknowledgments

Partial financial support has been received from the Agencia Estatal de Investigacion (AEI, Spain) and Fondo Europeo de Desarrollo Regional under Project PACSS RTI2018-093732-B-C21 (AEI/FEDER,UE) and the Spanish State Research Agency, through the Maria de Maeztu Program for units of Excellence in R&D (MDM-2017-0711). KK acknowledges funding from MINECO through the Ramón y Cajal program and through project SPASIMM, FIS2016-80067-P (AEI/FEDER, EU).