Competition for resources may reinforce the evolution of altruism in spatially structured populations

1.   Introduction

Though many models for the evolution of cooperation do incorporate mechanisms of competition, competition for space is usually only modeled implicitly. Such models tend to suggest that such competition (because it is also among relatives it effectively increases the cost of cooperation) would reduce the selection for cooperation ^[1,2]. However, approaches based on constructed fitness functions do not take into account all potential spatio-dynamical effects that act on populations ^[3]. Since consumer-resource interactions may have profound effects on the spatial patterns that result it is not so obvious what will be the net effect on selection. In the very least, it violates the implicit hypothesis of spatial homogeneity that underlies many kin selection models (which assume that genetic composition may vary spatially, but not densities or ecology in general).

In an earlier study using an explicitly spatial individual-based simulation method ^[4], we found that competition for resources can reinforce the evolution of altruism because it can segregate groups of related individuals, thus reinforcing kin selection and leading to more cooperation. As in Turing patterns ^[5], the spatial pattern depends strongly on the differences between the sizes of the kernels of competition and cooperation.

Here, we focus on the effects of competiton for resources on the evolution of cooperation using the framework of continuous-time probabilistic cellular automata ^[6,7]. This approach imposes certain constraints on the type of interactions that can be modeled. The fact that it is impossible to incorporate continuously varying local variables such as resource densities is offset by the fact that the resulting models are amenable to correlation dynamics (CD) approximations, which greatly extends the analytical toolbox that can be used to study them ^[8,9,10].

Most importantly, an analysis of the dynamics of spatial correlations allows to quickly assess the invasion capacity of arbitrary mutants in a resident system, which would require an excessive amount of simulations of the underlying cellular automaton. The advantage of having even an approximation of invasion fitness is that all the tools of adaptive dynamics theory at our disposal ^[11,12,13].

Here we will first construct a probabilistic cellular automaton (PCA) of individuals that live on a lattice, can invest in helping neighbours but also compete for resources necessary for reproduction. These resources are modeled as discrete units that appear randomly and diffuse through space. Then, a correlation dynamics (CD) approximation of this model is derived (including one or more resident types and a rare mutant) that is used to work out pairwise invasibility plots (PIPs) that allow to infer the ES level of investment in altruism. The difference with our previous study ^[4] is that here we assume that resource use characteristics do not evolve; we focus on the effect of resource competition per se.

The evolution of altruism is typically studied using variations of Hamilton's kin selection model ^[14]. A fitness function is set up for individuals that live in well-defined groups, which includes the putative effects of altruistic and competitive interactions on the fitness of a focal individual, taking into account the relatedness among the group members ^[17,18]. There exist various mathematical methods to work out the net effects as they percolate through the population ^[19,20] but most of these share the characteristics that (1) individuals live in discrete groups (metapopulations), (2) the underlying ecological interactions are not explicitly modeled and (3) fitness functions are assumed rather than derived. In particular the use of the Price equation ^[15] has allowed these methods to help analyse the potential consequences of the relative effects of altruism and competition ^[2]. However, in reality, populations are often organised in ill-defined and non-persistent groups that are not set up in advance but result from their ecological interactions, such as competition for space and resources. If the groups themselves are dynamic in terms of sizes and connections, methods based on explicit ecological mechanisms are needed to prise apart the selective effects of the various mechanisms ^[21].

Here we will concentrate on the effect of competition for resources on the evolution of altruism, assuming uptake rate to be a given parameter. By varying the uptake rate as well as the characteristics of the resource (such as turnover and diffusion rates) we can then study the effect of resource competition on the level of altruism and the relative effects of competition and cooperation.

The structure of this article is as follows. First we will present the basic ingredients of the probabilistic cellular automaton (PCA) that forms the starting point. Then we will work out the so-called Mean-Field approximation of this PCA, which is the set of ODEs that results when spatial pattern is not taken into account. This gives some insight in the ecological dynamics of the model and the effects of certain parameters.

We set up the model such that in the Mean-Field limit the consumer population has standard Lotka-Volterra dynamics. For simplicity, we assume that the resource has no population dynamics of its own; biologically, this resembles a consumer-detritus interaction where the resource items 'rain' onto the system; resource items that are present in the system then either decay or are consumed. Again for simplicity, we assume that the consumption of a single resource item is sufficient for reproduction. The consumers thus compete for space but also for resources; because the resources diffuse, local depletion will make itself felt to consumers further away. Because the consumers themselves do not move, they behave like plants or sessile organisms such as corals. The resulting model may be the simplest that is compatible with the concept of a competition kernel that is found in most population models defined in continuous space ^[22].

We continue by carrying out simulations of the PCA to gain insight in the spatial patterns that may arise. To study the evolutionary dynamics of the model, we derive the correlation dynamics (CD) equations of a multi-strain variant of the model. This multi-strain variant is then used to numerically analyse Pairwise Invasibility Plots (PIPs) to assess how ESS levels of altruism (if at all) depend on the ecological characteristics to the resource, which forms the novel aspect of this work. Finally, we discuss the results in relationship with previous work and suggest directions for future development.

2.   The Probabilistic Cellular Automaton

Space is represented by a (regular) network of sites (most simulations assume a square lattice where every site is connected to four or eight neighbours but other networks are possible too). A site may be empty (indicated by the symbol $\circ$), occupied by a resource item ($R$) or by one of a number ($n$) of phenotypes (indicated $\Phi_i$, where $i\in \{1, \ldots, n\}$). The events that modify the state of the lattice are the births and deaths of individuals (to keep things simple, individuals stay put and do not move), and the appearance, diffusion and decay of resource items. Note that there is no direct competition for resources among neighbours in the model; as in most biologically relevant settings, competition is indirect, because it causes a reduction in (local) resource density. The full set of events that define the model is given in Table 1 (for more information on how these events give rise to PCAs, see ^[10,23]).

The set of states a given site can be in is thus defined as

The density of a state is defined as the probability of finding a given site in this state and is indicated by $\left[{{x}}\right]$. In a similar fashion, the density of a pair of states $xy$ is defined as the probability of finding a given pair of neighbours in this combination of states and is denoted $\left[{{xy}}\right]$ (where $x, y \in \Omega$).

Events change the states of sites or pairs of sites on the lattice. The essential aspect of this PCA is that the mortality of an $i$ individual (if $i\in \Phi = \{ \Phi_1, \ldots, \Phi_n\}$) depends on its neighbours. Specifically, the death rate of the $i$ in an $ij$ pair depends on its neighbour $j$. Here, we will assume that

where $B(j)$ represents the benefit the $i$ individual receives from its $j$ neighbour (decreasing its mortality rate); the $i$ individual itself pays a cost $C(i)$ (increasing its mortality rate) to help all of its neighbours. A fully selfish phenotype is defined by $B(i) = C(i) = 0$ while any phenotype $\Phi_i$ with positive $B(i)$ (and associated $C(i)$) is to some degree altruistic. One of the main questions in evolutionary ecology is to understand which level of $B(i)$ (if positive at all) is favoured by natural selection.For notational simplicity, empty space and resource items are modeled as 'neighbours' that have no effect on mortality, that is $B(\circ) = C(\circ) = B(R) = C(R) = 0$. We have chosen this form over the alternative $d_i(j) = d_0(1+ C(i)- B(j))$ to avoid having to worry about negative mortality rates. When $C(i)$ and $B(j)$ are small relative to $d_0$, the two forms converge, which can be useful in the analysis, and to compare our results to other findings.

Note that although the notion of altruism was originally phrased in terms of cost and benefits associated with definite acts, such as whether to sacrifice oneself to save someone else, our analysis follows the modern convention of defining it in terms of more general fitness effects (see, e.g., ^[2]). Thus, an altruist is an individual that decreases its fitness in order to boost the fitness of one or more other individuals. In this setup, altruists always reduce their own fitness, irrespective of whether there are actual beneficiaries around or not.

It is not so difficult to derive the differential equations that govern the dynamics of the densities of single types. In correlation dynamics, however, one does not stop here; it is customary to consider (at least) changes in the densities of pairs of neighbours $\left[{{xy}}\right]$. However, before doing this we will consider what happens populations are (assumed to be) homogeneously distributed.

3.   Mean field dynamics

If space is ignored, a set of events like the one defined in Table 1 gives rise to a set of ordinary differential equations. The simplest way to obtain these is by assuming that the density of a pair of neighbours is the product of their densities, that is, by setting

and working out the bookkeeping of all the changes governing the $\operatorname{d\!}{}\left[{{x}}\right]/ \operatorname{d\!}{}{t}$ driven by the events in Table 1. More information about how this is done can be found in the cited literature (here, in particular, we have followed the formalism of ^[24]).

Under these conditions, resource dynamics is governed by

Since the consumption of one resource unit is directly linked to a reproduction event, the dynamics of phenotype $i$ is governed by

where the term [death rate of $i$] needs a bit more explanation. We assume that it depends on the individuals' social environments but to see the overall dynamics, let us assume for the moment that it is constant,

so that

which shows that in absence of social interaction competition for resources ($R$) is the sole mechanism that regulates population size (otherwise $\operatorname{d\!}{}{\left[{{ \Phi_i}}\right]}/ \operatorname{d\!}{}{}t$ can never be zero).

A differential equation for the density of empty sites, $\left[{{ \circ}}\right]$, could be derived in a similar fashion but this is not necessary as the density of empty sites follows from the conservation principle

More details can be found in Appendix A, but we report some elementary observations concerning the non-spatial model here. In absence of any consumption, the resource will settle at an equilibrium density of

This is an important compound parameter, and we will denote it $R^*$.

The presence of a single strain of consumers will cause it to settle at a stable equilibrium given by

provided there is enough resource, that is, if

The equilibrium may be significantly modified by social interactions, but these expressions give a useful first idea about the effects of the ecological parameters. In particular, it should be noted that, in spite of not affecting $R^*$, proportionally changing resource production and decay rates ($b$ and $\delta$) will affect consumer density ($\overline{{\left[{{ \Phi_1}}\right]}}$).

In this article we assume that every neighbour $y$ may have an effect on the focal individual's death rate. For simplicity, we model death rate as a pair event, such that the mortality rate of individual $x$ in an $xy$ pair not depends on its own phenotype but also on that of its neighbour $y$ (see Table 1). This implies that the overall mortality rate of $x$ phenotypes depend on the average densities of the phenotypes in the neighbourhood of the $x$'s. In the non-spatial version, these average neighbourhood densities are simply the global densities, so that

where $z$ is the number of neighbours ($z = 4$ in the default case) and $d_{ \Phi_i}(y)$ is given by Eq (2.2). This leads to

Note that this form of altruism differs slightly from the standard representation, where costs and benefits are either assumed to be measured in units of arbitrary 'fitness' or at least affect the same fitness-associated life history parameter (usually it is tacitly assumed that this is reproduction but it has been shown that such life history details may matter ^[25]). In the current setting, cooperative interactions specifically affect mortality rates only: the altruistic phenotyoe has a higher mortality rate but decreases that of its neighbours. This is not necessarily the most general case, but it is just as arbitrary as the other possible cases and it has the additional advantage of conceptually separating the effects of competition (affecting reproduction) from those of cooperative interactions (affecting survival).

4.   Altruists cannot invade a non-spatial population of nonaltruists

We can now consider the invasion rate of an altruistic strain in a non-spatial environment dominated by nonaltruists. For notational simplicity we can equate a strain's index with its investment in altruism (so that $\Phi = \{r, m\}$, $B(\Phi_m) = B(m) > 0$ and $C(\Phi_m) = C(m) > 0$, while $B(\Phi_r) = C(\Phi_r) = 0$). Since the $m$-type is much rarer than the $r$-type (implying $\left[{{m}}\right] \ll 1$), using Equations (3.3) and (3.11) we can deduce that the per-capita growth rate of the mutant ($m$) equals

which, given that the sum in brackets approximates unity if $\left[{{m}}\right]$ is very small (using Eq 3.6)), is approximated by

whereas the resident type ($r$), under the same conditions, has a per-capita growth rate of

which is always higher. Note that as the nonaltruist resident will be at equilibrium, the RHS of Eq (4.3) will be equal to zero (thus setting the resource density), which implies that the altruist's invasion rate (Eq (4.2)) is negative. Increasing the density of altruists would benefit the nonaltruists as much as it would the altruists, so the altruists will always do worse. However, it is well known that in spatial settings, the two types may segregate so that altruists may experience more altruists in their neighbourhood than do nonaltruists which may give them a sufficient advantage ^[10,26]. What is different here, is that not only altruistic support may vary locally, but also competition for resources.

If competition is only for space, a population will grow to fill up most of it, leading to a rather homogeneous overall distribution (with altruists reaching higher local densities than nonaltruists thanks to the net benefit they receive). Competition for resources may disrupt this pattern, as high-density patches resulting from much cooperation will necessarily also consume more resources, reducing population growth before space is full. In what follows we extend this observation by focusing specifically on how competition for resources affects the spatial dynamics of the different types.

5.   Spatial patterns

The spatial patterns generated by the cellular automaton are quite sensitive to certain parameters (as could be already expected from the non-spatial analysis) as the altruistic interactions that underlie the spatial dynamics give rise to a strong Allee effect ^[27].

Now we can add a second, mutant strain with a different level of altruism and show that these can, under certain circumstances, invade. As we have shown that altruists can never invade in non-spatial settings, such simulations show that spatial structure can, in principle, select for altruistic behaviour. However, finding out where evolution stops (i.e., what is the evolutionarily stable (ES) level of altruism) is not so obvious and would require us to repeat these simulations with different combinations of resident and mutant strains ^[13]. Then, to find out how the ES level depends on parameters would require repetition of this process for every possible combination of parameters that one wants to study, which amounts to an effort that is computationally intensive but still may yield litte insight. For instance, note when the altruist is rare, it is subject to strong demographic stochasticity ^[28,29], complicating fitness assessments. To carry out the adaptive dynamic analysis, therefore, we will derive a multi-strain version of the CD model. Even if explicit solutions do not exist, the resulting equations are much faster to solve numerically and more precise. More importantly, expressions for invasion fitness can be derived that give insight into the components of selection ^[3].

That spatial structure can favour altruism has of course been long known (see ^[3] for a review) but how this process is affected by resource competition is largely unexplored ^[4]. PCA simulations that we carried out show that often clusters of altruists are able to grow into and displace less altruist strains. A typical snapshot with such an invading cluster is shown in Figure 1. The temporal dynamics will be shown below, but if the simulation were continued long enough the altruists would eventually completely dominate the system. In what follows we will assess in particular the effect of resource dynamics on mutant invasion and the resulting ESS level of altruism using a Correlation Dynamics (CD) approximation, focusing on the effects of resource dynamics (turnover and diffusion). Then we will return to a PCA simulation to study the resulting spatial patterns.

6.   Correlation dynamics

The Correlation Dynamics method (often called the 'Pair Approximation Method' but this terminology is misleading ^[3,30]) produces a set of differential equations for spatial correlations (typically between neighbouring sites). These correlations allow the assessment the effects of spatial structure to a first approximation without having to resort to full-scale spatial simulations ^[8,10,23]. The correlation dynamics approach generates a rather large number of differential equations but they can be used to apply the basic methodology of adaptive dynamics.

We used a Mathematica (Wolfram Inc.) package that one of us developed (MvB) to derive and manipulate the Correlation Dynamics (CD) equations from the PCA specifications. These equations are described in Appendix B (the package we used, which is based on the methods outlined in ^[10,23,24], is available upon request).

The CD model is too complicated to be explicitly solved, but it can easily be numerically analysed. Since the CD model only tracks short-distance correlations, it will diverge from PCA dynamics if these involve significant larger-distance correlations, which are generated by many large-scale patterns ^[23]. However, as a first approximation, it can describe aspects of the dynamics that the MF approximation definitely cannot. Most significantly, altruists can never invade in the non-spatial version whereas the CD correctly predicts the mutant will invade (Figures 1 and 2). Of course, the PCA model is stochastic and will show demographic stochasticity whereas the CD approximation is deterministic. The idea is that increasing the system (i.e., lattice) size will decrease the stochastic fluctuations, leaving only the systematic error of the CD approximation.

As can be seen in Figure 2, compared with explicit simulations the per-capita rate of increase of the mutant derived from the CD model (after a transition period to allow the resident system to equilibrate) gives a good indication of invasion fitness,

for a sufficiently large stretch of time $T$ of a time series $\left[{{m}}\right]_t$ from a numerical simulation of the CD model. (In mean field situations this produces equations (4.2) and (4.3), but in spatially structured populations more complicated expressions arise ^[3].) Here we will use definition (6.1) to numerically calculate a Pairwise Invasibility Plot (PIP), a standard tool from Adaptive Dynamics ^[11,13]. From the PIP that corresponds to the basic set of parameter values used previously (Figure 3), and assuming, in addition, that the costs and benefits are related through a public goods game scheme $B(\Phi_i) = 10\; C(\Phi_i)$ ^[16] we can infer that the ES level of altruism will be around $C(\Phi^*) \approx 0.06$. This implies a significiant investment in altruistic behaviour: isolated individuals have a $6\%$ higher mortality rate than nonaltruistic individuals. Altruistic behaviour is thus only favoured when an altruist is helped by sufficiently many other altruists in its neighbourhood.

7.   Invasion fitness

Useful insight is obtained by analysing the invasion conditions of arbitrary mutants strains. A mutant phenotype $m$ will invade if its per capita fitness

is positive, given that the system is dominated by the resident type $r$. Here $\left[{{x|m}}\right]$ stands for the local density of $x$, that is, for the probability of finding an $x$ at a given neighbouring site of a mutant (which may be significantly different from the global density $\left[{{x}}\right]$). Note that per-capita fitness measures like Eq (7.1) are not identical to 'true' invasion fitness (Eq (6.1)) but they do have the same sign structure ^[25], implying that invasion fitness is positive whenever per-capita fitness is positive. Note also that in order to work out the sign of Eq (7.1) we need to know the local densities $\left[{{ R|m}}\right]$ and $\left[{{m|m}}\right]$, and these are difficult to calculate exactly except in special cases ^[8,10,21]. Here, we will use numerical approximations obtained from the full CD model.

Figure 4 shows the advantage of using the CD approach when studying invasion. Because a mutant is initially rare, it is subject to substantial demographic stochasticity and as a consequence its invasion fitness fluctuates widely (panel A); the invasion fitness of the resident (which has a much bigger population) fluctuates much less (panel B). The only way to reduce these fluctuations is to increase the system size, increasing simulation time. The CD model sidesteps the problems associated with demographic stochasticity, as it assumes that system size is infinite.

The mutant will be able to invade when its per-capita fitness $F_r(m)$ (given by Eq (7.1)) is positive. In a system with only selfish residents ($C(r) = B(r) = 0$), altruistic mutants can invade non-altruist populations if

or, using that $\left[{{ \circ|m}}\right] + \left[{{ R|m}}\right] + \left[{{r|m}}\right] = 1 - \left[{{m|m}}\right]$,

where the right hand side is actually the version of Hamilton's Rule that incorporates the demography of the modeled population ^[10]. This can be seen most easily when it is assumed that costs and benefits are small (both $B(m) \ll 1$ and $C \ll 1$) as then the condition converges to

which is the version of Hamilton's Rule that results if costs and benefits linearly affect mortality (rather than reproduction, as is usually assumed).

The probability that a neighbouring site is occupied by another mutant, $\left[{{m|m}}\right]$, thus takes on the role of the coefficient of relatedness (cf. ^[10]). If the mutant has no effect on its (local) environment the standard version of Hamilton's Rule emerges when Condition (7.4) is analysed. Then the mutant should pick costs $C(m)$ and benefits $B(m)$ such that

is maximised for a fixed $\left[{{m|m}}\right]$. However, though the coefficient of relatedness is often assumed to be a fixed and given constant, in spatially extended systems it is something the mutant itself can influence through its use of space, that is, $\left[{{m|m}}\right]$ depends on the characteristics of the mutant's clusters when it invades ^[3,10]. Furthermore, in the model we are considering here, the mutant may also have an effect on local resource dynamics, which affects the left hand side of Condition (7.3).

Any trait $T_m$ of the mutant that affects its vital rates would affect the local densities $\left[{{ R|m}}\right]$ and $\left[{{m|m}}\right]$, which adds terms to the optimisation condition ^[3]. In standard evolutionary approaches, local environmental variables (including relatedness parameters) are assumed to be constant (which typically demands the use of qualifiers such as 'small effect' when talking about mutations). The real challenge, however, is to gain insight into how mutants influence their own environment, which becomes all the more urgent when one wants to discuss the relative effects of competition and cooperation. Even with extensive simulation of variations of the PCA model it may be difficult to synthesise their effects into a comprehensive picture.

Optimisation charicature

If it is assumed that the mutant does not influence its own environment, maximisation of invasion fitness amounts to minimisation of the mortality rate. Evaluating $\operatorname{d\!}{} F_r(m) / \operatorname{d\!}{} m = 0$ then leads to the condition

where the prime stands for differentiation with respect to $m$. Setting $B(m) = \alpha C(m)$ (so that the individuals essentially play a public goods game with multiplication factor $\alpha$ ^[16]) this leads to a quadratic equation in $C(m)$ that can be solved explicitly (see Appendix C). The solution is rather messy, and depends on $\alpha$ and $\left[{{m|m}}\right]$ in a complicated way (Figure 5). However, there is a pattern. First, the investment is always zero for low relatedness ($\left[{{m|m}}\right] < 1/\alpha$), and increases monotonically. Investment depends also on the multiplication factor $\alpha$; for large values of relatedness the optimum investment decreases with $\alpha$ (because then a smaller investent already produces a big enough benefit). However, the relatedness $\left[{{m|m}}\right]$ typically is lower than $1/z$ ^[10] so that we do not expect high levels of altruism to evolve. This agrees well with the ES level of altruism that was predicted by the PIP shown in Figure 3.

Assuming the local environment is fixed gives useful insight into the basic balance of costs and benefits, but does not yet take the fact into account that strains may modify the local spatial pattern and thus affect local densities $\left[{{x|m}}\right]$ they experience. Local resource density obviously has a direct effect as the resource uptake rate (which equals birth rate in our model) of the mutant equals $u \left[{{ R|m}}\right]$) but as our earlier work suggested, by segregating the clusters, resource competition may have the additional effect of decreasing $\left[{{r|m}}\right]$ relative to $\left[{{m|m}}\right]$.

In principle, one would like to use the correlation dynamics model to calculate the local densities but in practice this is not feasible as only but the simplest CD models allow explicit calculation of local densities ^[8,23]. Our previous studies suggested that clusters even when they are spatially segregated may still compete for resources if the kernel of competition is sufficiently large. In the present framework this implies that resource dynamics may play a role: in the PCA the direct kernel of competition is the same as the kernel of cooperation (the $z = 4$ Von Neumann neighbourhood) but competition for resources may be felt farther away when resources diffuse, resulting in a competition kernel that is effectively larger.

Analysis of the non-spatial ('mean field') version of the model suggests that the equilibrium resource density in absence of consumption, $R^*$, is only dependent on the relative values of resource productivity and decay (Eq (3.7)). However, the density of the consumers (Eq (3.9)) does also depend on the absolute resource productivity so that consumer density changes even if the resource productivity/decay ratio stays constant. And indeed, PIPs for different values of $b$ and $\delta$ but with the same ratio may lead to strongly different level of ES altruism (Figure 6A). Slowing down the overall turnover rate $b + \delta$ intensifies competition, as food items diffuse farther during their existence, effectively increasing the competition kernel. Thus a harsher environment results, but the reduced consumer densities that result allow kin selection to compensate to a certain extent, favouring increased levels of altruism. Varying the resource diffusion rate $D$ has a similar effect, but less strong (Figure 6B).

Based on these insights one might expect that increasing the uptake rate will have a similar effect (as these are also likely to decrease resource levels) but this is actually not the case. The reason is that increasing uptake rate also increases the rate of reproduction. It allows the consumers to make more efficient use of the scarce resources and this more than compensates for the reduced resource density. The effect is that the global consumer density will increase, clusters are more likely to mix and individuals are less likely to be surrounded by of relatives. This seems to hold for any ecological parameter that boosts population growth (results not shown). Altruism is only maintained by natural selection under relatively harsh ecological conditions, preventing a population to grow to high densities and occupy all space.

8.   Characteristic cluster size

Though a cursory glance at the spatial distributions suggests they seem rather random they are in fact highly structured. Consider the distribution shown in Figure 7A, which resulted from a simulation with six resident phenotypes that are selectively neutral so that the colour pattern in the Figure gives an indication of the relatedness structure in the population. The phenotypes are definitely not well-mixed, and seem to form characteristic clusters (an animation of this simulation is available as Electronic Supplementary Information). If we pick one of the strains and count its density in plots of varying sizes (here, we use windows of form of 'taxicab' (or 'Manhattan') circles with different taxicab radiuses ^[33]), we observe that there is a great variance among smaller plots (some contain many individuals, others few, see Figure 7B) but as the plots grow to a certain size the variance suddenly decreases (Figure 7C), suggesting that a characteristic lenght scale of phenotype clusters of, in this case, about $4$ or $5$ sites (see ^[31,32] for more information about the so-called 'coherence scale' of assessing characteristic spatial structures).

Another way of gaining insight is measuring the densities of food items and of the different phenotypes increasing distances away from a given type (Figures 8A–C). In the simulation with $6$ phenotypically identical strains, one can see that in the close neighbourhood of an individual, most space is occupied by individuals of the same strain and that the density of food items is much lower. Only a number of steps away the densities start to converge to the global averages. This reinforces the insight that in spatially explicit systems, what individuals encounter may be very different from the global average. Not only do they interact mainly with relatives, the density of resources and free space may be different.

Different strains are spatially segregated, as often reported, but they still do compete for space and resources. To assess the intensity of this competition, however, is more complicated than just analysing a spatial pattern, as it is a highly dynamical process that involves demographic processes and the generation, diffusion, consumption and decay of the resource. To our knowledge no mathematical method is currently available to analyse characteristic spatio-temporal patterns.

9.   Discussion

The results we present here reinforce our earlier observation ^[4] that competition for resources may have unexpected effects on the evolution of altruism, at least from a standard kin selection perspective. Traditionally, kin selection models are based on evaluating the effects of various interactions on the fitness of an average member on the population (the 'focal individual') ^[19,20]. Competition, in this framework, is assumed to have an overall negative effect on the fitness of such focal individuals. An individual does not only benefit from its altruistic relatives, it also competes with them, and because of their proximity even more intensely than with unrelated individuals. Competition is then expected to lower the investment in altruistic behaviour ^[1,2].

That the relative scales of cooperation and competition are important has been shown before ^[18,40]. If their spatial scales are equal, the benefits of cooperation and the costs of competition tend to cancel out, but if the scale of competition is larger than that of cooperation, a group of related individuals may nevertheless reap the benefits of increased cooperation, because they may make their presence felt farther out (and can thus 'export' their characteristics ^[34,35]).

In our earlier work ^[4], we observed that competition may reinforce altruism under certain conditions, because it can increase the spatial segregation of unrelated individuals, thus preventing altruistic clusters from being exploited by non-related cheater individuals. The study presented here lends more support to this idea, which is important because the ecological setup is rather different from that of our earlier study (the most important of which are that in our previous work resource uptake evolved whereas here it is a fixed trait).

What is clear from our results is that competition for resources may affect the evolution of altruism, but not always in immediately intuitive ways. For instance, if the resource uptake rate goes up (which would imply that clusters of relatives may make their presence felt further away), the ESS investment in altruism goes down (Figure 6A). The reason for this effect is that the increased uptake rate leads to higher overall densities, which counteracts the effect of segregation. Higher uptake rates benefit all strains equally; different evolutionary outcomes might result if strains can differentially invest in uptake rates, as indicated by our previous study ^[4]. This is something we leave for a future study.

We have demonstrated how one can apply the theory of characteristic lenght scales ^[31,32] to assess the characteristic cluster size of a rare mutant strain. However, as this involves deciding when certain variances drop below a given threshold, it is not so obvious how to give a more precise mathematical definition. At the moment it is not possible to analyse how phenotypic traits affect characteristic cluster size except by simulation. Obtaining analytical insight into this link is worthy of a study in itself.

What we did not observe, like we did in our previous study ^[4], is the emergence of a clear quasi-stationary Turing-like pattern. There are reasons for expecting such patterns as the distribution is generated by positive and negative feedback processes operating on different spatial scales ^[36,37]. The reason for this is probably that the emergent structures remain relatively small and are therefore strongly subject to demographic stochasticity. Nevertheless, the Allee effect that results from the altruistic interactions and which punishes isolated individuals makes the clusters drift around much less than expected by chance. The size of the clusters and the degree of mixing strongly depends on ecological parameters. The less 'harsh' the environment is, the higher overall densities and the more mixing of strains. The evolutionary consequences are demonstrated in the sequence of results displayed in Figure 6A, that is, reduced altruism.

We have made various simplifying assumptions in order to keep our model as tractable as possible. Real systems are inevitably more complicated. For instance, we assumed that the resource is biologically inert, but in reality resources are often other species (or produced by them), and these have spatial population dynamics of their own. How these dynamics will affect the evolutionary outcome is an open question, but the present study shows, if anything, that the properties of such resources will have an effect.

In our previous study ^[4], we also observed a strong interaction between the evolution of resource use and the evolution of altruism: under certain conditions, altruism did not evolve if resource use does not evolve in tandem. Here, with resource uptake rate being fixed and given, a certain level of altruism always evolves. Nevertheless, our results suggest that the level of altruism does depend on uptake rate (results not shown) so the outcome may change if resource uptake evolves too. We leave this for a future study.

To study evolution of cooperative behaviour it has become customary to construct a priori variations of Hamilton's Rule (often in a framework based on Price's Equation ^[38,39], but there exist other formal systems as well ^[19]). However, these do not capture all spatio-dynamical effects, in particular the ecological mechanisms that actually create the group structures and relatedness patterns in spatially extended populations ^[3]. What is needed is a more comprehensive insight into how cooperative behaviour affects the ecological mechanisms that create the spatial patterns and the spatial scales they imply. In our study, competition for resources reinforces kin selection because the ecological dynamics happens to segregate the clusters of relatives that compete, but we still do not have a good toolbox to study these mechanisms. It is important to develop such tools, as our results clearly demonstrate that the underlying ecology should be considered as well.

Acknowledgments

MvB should acknowledge support of program Investissements d'Avenir launched by the French Government and implemented by ANR with the references ANR-10-LABX-54 MEMOLIFE and ANR-11-IDEX-0001-02 PSL Research 432 University. MvB also expresses his gratitude to the Institute des Hautes Etudes Scientifiques (IHES, Bures-sur-Yvette, France) for its hospitality.

Conflict of interest

The authors declare no conflicts of interest in this paper.

A.   Non-spatial resident dynamics

The basic ecological model (non-spatial and non-evolving) is governed by two differential equations,

if the death rate is kept constant.

Working out the equilibrium condition $\operatorname{d\!}{}\left[{{R}}\right]/ \operatorname{d\!}{} t = \operatorname{d\!}{}\left[{{ \Phi}}\right]/ \operatorname{d\!}{} t = 0$ yields

This equilibrium is always stable.

Adding the altruistic interactions introduces density dependence in the mortality rate, which creates an Allee effect.

The consequences of this Allee effect are that potentially there are now three equilibria, two stable ones separated by an unstable one. This bistability causes the interesting hysteresis effects reported in some studies ^[27].

B.   Correlation Dynamics equations

The idea behind the Correlation Dynamics approach is that one can write down differential equations for the states of pairs of neighbouring sites, just as one can do for single sites, which gives the dynamics of the mean densities. For the sake of argument consider a system with only a single population which for conciseness we will denote $x$ here. Then there are four types of pairs, $\circ \circ$, $\circ x$, $x \circ$ and $xx$. The densities of $\circ x$ and $x \circ$ pairs are necessarily the same, so we need only one equation for those.

Take the events that can destroy an $\circ x$ pair: the $x$ can reproduce into the empty site (thus producing an $xx$ pair) --- but the same outcome also results if another $x$ among the $z-1$ other neighbours reproduces into the empty site (neighbours of other types can reproduce into the empty site too). If the density of $\circ x$ pairs is denoted $\left[{{ \circ x}}\right]$ (that is the probability of finding a random pair in the $\circ x$ state), then this density will therefore decrease with an overall rate

where $\left[{{x}}|{{ \circ x}}\right]$ stands for the probability that a given other neighbour of the $\circ$ site is also occupied by an $x$.

In this simple version $\circ x$ pairs can be created in two ways: either the left $x$ dies in an $xx$ pair or a third neighbour reproduces into the right empty site of a $\circ \circ$ pair (we do not incorporate movement or dispersal in this model). Thus the components of the differential equations resulting from processes that create $\circ x$ pairs are

and

giving

The other equations are constructed in a similar fashion (in practice we use a Mathematica (TM) package for doing that).

At this point, the equations depend on pair densities ($\left[{{ \circ \circ}}\right]$, $\left[{{ \circ x}}\right]$, $\left[{{x \circ}}\right]$ and $\left[{{xx}}\right]$) but also on conditional densities like $\left[{{x}}|{{ \circ x}}\right]$. The latter indicates the conditional probability that a given other neighbour of the empty site in an $\circ x$ pair is occupied. Formally, this can be expressed as

which shows this quantity also depends on triple densities (i.e., the frequencies of chains of three neighbours in particular states). To avoid having to derive differential equations for the triple densities (which would not have solved the problem, as these depend on higher order structures again), most studies 'close' the system by assuming that the conditional probability can be approximated by ignoring the state of the far end,

which results in an expression in terms of the (known) frequencies of singletons and pairs, and which is what we did in this article (the singletons can be calculated easily from the pairs).

Even when it is taken into account that symmetric pairs have the same dynamics, the number of differential equations can get fairly big. If the number of states is $N$ number if state variables becomes $(N+1)N/2$, so the model with mutant and resident, which has $z = 4$ states, necessitates the construction of $10$ differential equations). For more information, see the references cited.

C.   Optimum altruism given fixed environments

The optimum level of altruism for a strain $\Phi_m$ invading a fully selfish world can be easily calculated if $\left[{{R|m}}\right]$ and $\left[{{m|m}}\right]$ are known and fixed, let us say, $\left[{{R|m}}\right] = \bar R$ and $\left[{{m|m}}\right] = \bar \rho$ ($\bar \rho$ thus standing for the coefficient of relatedness). Assuming altruism is a scalar trait, given by $m$, that affects $C(m)$ and $B(m)$ (but not $\left[{{R|m}}\right]$, nor $\left[{{m|m}}\right]$), we can take the derivative the mutant's invasion fitness (given by Eq (7.1)) with respect to $m$, to obtain

where the primes denote derivatives by $m$.

The optimum $m^\circ$, if it exists, should then satisfy

After some algebra this results in the condition

Making the public games assumption $C(m) = m$ and $B(m) = \alpha m$ finally results in a quadratic expression in $m$,

which can be solved explicitly, and whose positive solution is

which has been used to draw Figure 5.

Note that on a larger timescale, local densities will depend on $m$, but we do not have expressions for these. The solution (14.5) should hold if costs and benefits are not too big.