Spatial adaptive dynamics

In the following, we assume that the total population size is constant, and that spatially structured populations are represented by lattices in which each site is occupied by one individual. Each individual interacts with a limited number of local neighbours, and we assume this number, $k$, to be the same for all individuals. We first consider a case where there are two types of players in the structured population: a mutant type with trait value $y$, and a resident type with trait value $x$ (where $x$ and $y$ denote investment strategies in a continuous game). If the mutant has $j$ other mutants among its $k$ neighbours, the mutant payoff is $$\begin{array}{r}{\pi }_m(j)=[(k-j)P(y,x)+jP(y,y)]/k.\end{array}$$ Similarly, the payoff of a resident with $j$ mutant neighbours is given by $$\begin{array}{r}{\pi }_r(j)=[jP(x,y)+(k-j)P(x,x)]/k.\end{array}$$ The payoffs, ${\pi }_m(j),{\pi }_r(j)$, of mutants and residents from interactions with their $k$ neighbours determines the birth rates as $b_m(j)=\exp (w{\pi }_m(j))$ and $b_r(j)=\exp (w{\pi }_r(j))$, where $w>0$ denotes the strength of selection. The birth rate is proportional to the probability of taking over an empty site for which a given mutant or resident individual competes. For $w\ll 1$ selection is weak and differences in payoffs result in minor differences in birthrates, and hence in small differences in probabilities of winning competition for an empty site. With strong selection, $w\gg 1$, payoff differences are amplified in the corresponding birthrates. This exponential payoff-to-birthrate map has several convenient features:

1. ensures positive birthrates,
2. admits easy conversion to probabilities for reproduction,
3. selection can be arbitrarily strong,
4. for weak selection the more traditional form of birthrates $b_i(j)\approx 1+w{\pi }_i(j)$ with $i=m,r$, is recovered.

Note that in the limit $y\to x$ differences in birthrates vanish. However, this does not imply weak selection. Instead, selection strength is determined by the magnitude of the invasion fitness gradient at $x$, which is proportional to $w$.

In well-mixed populations the current state is simply given by the frequency of mutants and residents, respectively. In contrast, in structured populations the state space is immense because it involves all possible configurations. Pair approximation offers a convenient framework to account for corrections arising from spatial arrangements. Instead of simply tracking the frequencies of mutants and residents, pair approximation considers the frequencies of neighbouring strategy pairs. We denote the frequencies of mutant-mutant, mutant-resident, resident-mutant and resident-resident pairs by $p_{mm},p_{mr},p_{rm}$, and $p_{rr}$, respectively, with $p_{mm}+p_{mr}+p_{rm}+p_{rr}=1$. The frequency of mutants is given by $p_m=p_{mm}+p_{mr}$ and that of residents by $p_r=p_{rr}+p_{rm}$. For consistency, the frequency of mutants can also be expressed as $p_m=p_{mm}+p_{rm}$ by considering links that point towards mutants rather than away from them. Consequently $p_{rm}=p_{mr}$ must hold and the dynamics is reduced to two dynamical variables.

The most informative quantities are the global mutant frequency, $p_m=p_{mm}+p_{mr}$, and the local mutant density $q_{m|m}=p_{mm}/p_m$, i.e. the conditional probability that a neighbour of a mutant is also a mutant. Note that for rare mutants, $p_m\ll 1$, their local densities need not be small. The derivation of the corresponding dynamical equations depends on the details of the microscopic updating.

Death-birth process

For death-birth updating, first the birthrate of all individuals is calculated by taking the payoffs from interactions with their $k$ neighbours into account (we assume regular spatial structures in which each individual has the same number of neighbours $k$). An individual is then uniformly at random selected to die, and its $k$ neighbours compete to repopulate the vacant site with their offspring. The probability of success is proportional to the birthrate. We note that this setup implies that competition occurs at the local scale, and that life expectancy is the same for all individuals. This is in contrast to birth-death updating where competition occurs globally and individuals with successful neighbours have a lower life expectancy.

Note that payoffs, and hence birth rates, are based on interactions with all neighbours, including the neighbour that may subsequently be chosen to die (uniformly at random) and its vacant site subject to recolonization by the offspring of one of its neighbours. To determine the dynamics of $p_m$ and $q_{m|m}$, we first note that configurations only change when a resident is replaced by a mutant, or when a mutant is replaced by a resident.

The frequency of mutants and mutant-mutant pairs increases whenever a resident dies and a mutant neighbor repopulates the vacated site. For a resident with $j$ mutant neighbours this happens with probability: $$\begin{array}{}\text{(1)}& T_{\text{db}}^{+}(j)=& (1-p_m)\cdot (\genfrac{}{}{0ex}{}{k}{j})q_{m|r}^j(1-q_{m|r}{)}^{k-j}\frac{j{b}_m((k-1)q_{m|m})}{j{b}_m((k-1)q_{m|m})+(k-j)b_r((k-1)q_{m|r})}\end{array}$$ with $q_{m|r}=(1-q_{m|m})p_m/(1-p_m)$, and where $b_r(v),b_m(v)$ denote the birthrates of residents and mutants with an average number of $v$ mutants among their neighbours. The first term in Eq. $\text{(1)}$ indicates the probability that a resident dies, the second denotes the probability that it had $j$ mutant neighbours and the last term is the probability that one of them reproduces. Similarly, the frequency of mutants and mutant-mutant pairs decreases if a mutant dies and one of its resident neighbors reproduces. For a mutant with $j$ mutant neighbours this happens with probability: $$\begin{array}{}\text{(2)}& & T_{\text{db}}^{-}(j)=p_m(\genfrac{}{}{0ex}{}{k}{j})q_{m|m}^j(1-q_{m|m}{)}^{k-j}\frac{(k-j)b_r(1+(k-1)q_{m|r})}{j{b}_m(1+(k-1)q_{m|m})+(k-j)b_r(1+(k-1)q_{m|r})}.\end{array}$$ All other transitions do not alter the composition of the population, and hence the $T_{\text{db}}^{\pm }(j)$ for $j=0,...,k$ define the rate of change of the frequency of mutants: $$\begin{array}{}\text{(3)}& {\dot{p}}_m& =\sum _{j=0}^k(T_{\text{db}}^{+}(j)-T_{\text{db}}^{-}(j)).\end{array}$$ In order to derive the rates of change of $q_{m|m}$ it helps to start with changes in $p_{mm}$. First, consider a resident with $j$ mutant neighbours that has been successfully replaced by a mutant. This happens with probability $T_{\text{db}}^{+}(j)$ and increases the number of $mm$-pairs by $j$ or, equivalently, their frequency $p_{mm}$ by $2j/(Nk)$, where $Nk/2$ denotes the total number of undirected links in a regular graph of size $N$ and degree $k$. Similarly, with probability $T_{\text{db}}^{-}(j)$ a mutant with $j$ mutant neighbours is replaced by a resident, which reduces the frequency $p_{mm}$ by $2j/(Nk)$. Thus the rate of change of $p_{mm}$ becomes: $$\begin{array}{}\text{(4)}& {\dot{p}}_{mm}& =\sum _{j=0}^k\frac{2j}{k}(T_{\text{db}}^{+}(j)-T_{\text{db}}^{-}(j)),\end{array}$$ where the term $1/N$ has been omitted because the rates of change in both $p_m$ and $p_{mm}$ are proportional to the inverse population size and hence this factor can be absorbed through a constant rescaling of time. Finally, using ${\dot{q}}_{m|m}=({\dot{p}}_{mm}-q_{m|m}{\dot{p}}_m)/p_m$ results in: $$\begin{array}{}\text{(5)}& {\dot{q}}_{m|m}& =\frac{1}{p_m}\sum _{j=0}^k(T_{\text{db}}^{+}(j)-T_{\text{db}}^{-}(j))(\frac{2j}{k}-q_{m|m}).\end{array}$$ This yields the dynamical equations for the global, Eq. $\text{(3)}$, and local, Eq. $\text{(5)}$, frequencies of mutants.

Spatial invasion fitness

The invasion fitness of a rare mutant trait $y$ in a resident $x$ is defined by $f(x,y)=\frac{{\dot{p}}_m}{p_m}$ in the limit $p_m\to 0$. We first note that the leading order of the global dynamics, Eq. $\text{(3)}$, is $O(p_m)$, while the local dynamics, Eq. $\text{(5)}$, scales with $O(1)$ and hence happens much faster when mutants are rare, $p_m\ll 1$. This results in a convenient separation of time scales, so that we can use the equilibrium $q_{m|m}^{\ast }$ of Eq. $\text{(5)}$ in Eq. $\text{(3)}$ to calculate the invasion fitness. In the limit $p_m\to 0$, the sum over $T_{\text{db}}^{+}(j)$ in Eq. $\text{(3)}$ and divided by $p_m$ reduces to: $$\begin{array}{}\text{(6)}& \underset{p_m\to 0}{lim}\frac{1}{p_m}\sum _{j=0}^k{T}_{\text{db}}^{+}(j)& =\frac{k(1-q_{m|m}^{\ast })b_m((k-1)q_{m|m}^{\ast })}{b_m((k-1)q_{m|m}^{\ast })+(k-1)b_r(0)}\end{array}$$ using $q_{m|r}=p_m(1-q_{m|m}^{\ast })/(1-p_m)$. In contrast, the sum over $T_{\text{db}}^{-}(j)$ simplifies only marginally by cancelling the common factor $p_m$ and $q_{m|r}\to 0$. Thus the invasion fitness becomes: $$\begin{array}{}\text{(7)}& f(x,y)& =\frac{k(1-q_{m|m}^{\ast })b_m((k-1)q_{m|m}^{\ast })}{b_m((k-1)q_{m|m}^{\ast })+(k-1)b_r(0)}-\sum _{j=0}^k(\genfrac{}{}{0ex}{}{k}{j}){q_{m|m}^{\ast }}^j(1-q_{m|m}^{\ast }{)}^{k-j}\frac{(k-j)b_r(1)}{j{b}_m(1+(k-1)q_{m|m}^{\ast })+(k-j)b_r(1)}.\end{array}$$ In order to calculate $q_{m|m}^{\ast }$, we first note that in the limit of rare mutants, $p_m\to 0$, and for mutant traits $y$ close to the resident trait $x$, Eq. $\text{(5)}$ somewhat simplifies to $$\begin{array}{rl}{\dot{q}}_{m|m}& =\frac{k(1-q_{m|m})b_m((k-1)q_{m|m})}{b_m((k-1)q_{m|m})+(k-1)b_r(0)}(\frac{2}{k}-q_{m|m})-\\ \text{(8)}& & \sum _{j=0}^k(\genfrac{}{}{0ex}{}{k}{j})q_{m|m}^j(1-q_{m|m}{)}^{k-j}\frac{(k-j)b_r(1)}{j{b}_m((k-1)q_{m|m}+1)+(k-j)b_r(1)}(\frac{2j}{k}-q_{m|m}).\end{array}$$ Second, a Taylor expansion of the right-hand-side of Eq. $\text{(8)}$ in $y$ around $x$ yields, up to first order: $$\begin{array}{rl}{\dot{q}}_{m|m}=& \frac{1}{k^3}(1-q_{m|m})(2k^2(1-(k-1)q_{m|m})+w(y-x)(k-1)(1-q_{m|m})(2-(k-4)q_{m|m})\times \\ \text{(9)}& & ((k-1)q_{m|m}{\partial }_{z}P(x,z)+k{\partial }_{z}P(z,x)){|}_{z=x})+O((y-x{)}^2).\end{array}$$ Thus, in this approximation, the equilibrium $q_{m|m}^{\ast }$ is given by the roots of a third order polynomial (plus the trivial, uninteresting root $q_{m|m}=1$). To circumvent further analytical challenges, the zeroth order approximation of $q_{m|m}^{\ast }$ in $y$ near $x$ can be obtained by solving for the roots of Eq. $\text{(9)}$ for $y=x$, which yields $q_{m|m}^{\ast }=q_{m|m}^{\ast (0)}+O(y-x)$ with $q_{m|m}^{\ast (0)}=1/(k-1)$.

Next, the first order approximation of $q_{m|m}^{\ast }$ is obtained by implicit differentiation of Eq. $\text{(9)}$ with respect to $y$ (keeping in mind that $q_{m|m}$ is a function of $y$) and evaluation at $y=x$. Setting the expression to zero yields an equation for the zeroth and first order coefficients, $q_{m|m}^{\ast (0)}$ and $q_{m|m}^{\ast (1)}$, of the Taylor expansion at the equilibrium $q_{m|m}^{\ast }$ for $y$ near $x$: $$\begin{array}{}\text{(10)}& -\frac{1}{k^3}(2k^2(k-2(k-1)q_{m|m}^{\ast (0)})q_{m|m}^{\ast (1)}-w(1-q_{m|m}^{\ast (0)}{)}^2(k-1)(2-(k-4)q_{m|m}^{\ast (0)})((k-1)q_{m|m}^{\ast (0)}{\partial }_{z}P(x,z)+k{\partial }_{z}P(z,x)){|}_{z=x})& =0.\end{array}$$ Solving for the first order coefficient $q_{m|m}^{\ast (1)}$ using $q_{m|m}^{\ast (0)}=1/(k-1)$ then yields: $$\begin{array}{}\text{(11)}& q_{m|m}^{\ast (1)}=& w\frac{k^2-4}{2(k-1{)}^2{k}^2}(k{\partial }_{z}P(x,z)+{\partial }_{z}P(z,x)){|}_{z=x}.\end{array}$$ Assembling all the pieces finally results in the first order expansion in $y$ of the local pair density equilibrium, $q_{m|m}^{\ast }$, around $x$. The upshot is that this allows to simplify the invasion fitness to a function of the mutant and resident traits only. The invasion fitness of mutants, $f(x,y)$, defined as their per capita growth rate, ${\dot{p}}_m/p_m$, in the limit $p_m\to 0$, then becomes: $$\begin{array}{}\text{(12)}& f(x,y)=& \frac{k(1-q_{m|m}^{\ast })b_m((k-1)q_{m|m}^{\ast })}{b_m((k-1)q_{m|m}^{\ast })+(k-1)b_r(0)}-\sum _{j=0}^k(\genfrac{}{}{0ex}{}{k}{j})(q_{m|m}^{\ast }{)}^j(1-q_{m|m}^{\ast }{)}^{k-j}\frac{(k-j)b_r(1)}{j{b}_m(1+(k-1)q_{m|m}^{\ast })+(k-j)b_r(1)},\end{array}$$ where $b_m(v)$ and $b_r(v)$ denote the birth rates of mutants and residents, respectively, with an average number of $v$ mutants in their neighbourhood.

Even though the solution to ${\dot{q}}_{m|m}=0$ is analytically inaccessible, in general, the equilibrium $q_{m|m}^{\ast }$ can be approximated using a Taylor expansion if $|y-x|\ll 1$: $$\begin{array}{}\text{(13)}& q_{m|m}^{\ast }=& \frac{1}{k-1}+w(y-x)\frac{k^2-4}{2(k-1{)}^2{k}^2}({\partial }_{z}P(x,z)+k{\partial }_{z}P(z,x)){|}_{z=x}+O((y-x{)}^2).\end{array}$$ It follows that in the limit $y\to x$, mutants with at least one resident neighbour have, on average, one mutant neighbour among their $k-1$ other neighbours. Note, mutants with no resident neighbours are uninteresting because they are unable to initiate a change in the population configuration. Interestingly, this limit of the local pair configuration is fairly robust with respect to changes in the updating process (c.f. Eq. $\text{(26)}$ for birth-death updating). Moreover, in this limit a rare mutation with positive invasion fitness is guaranteed to eventually take over.

Selection gradient, convergent and evolutionary stability

Using Eqs. $\text{(12)}$ and $\text{(13)}$ the selection gradient, $D_{\text{db}}(x)=\frac{\partial f(x,y)}{\partial y}{|}_{y=x}$, as well as its Jabobian, $C{S}_{\text{db}}(x^{\ast })=\frac{d{D}_{\text{db}}(x)}{dx}{|}_{x=x^{\ast }}$, and the Hessian of fitness, $E{S}_{\text{db}}(x^{\ast })=\frac{{\partial }^{2}f{(}^{\ast },y)}{\partial y^2}{|}_{y=x^{\ast }}$, at a singular point $x^{\ast }$ can be calculated as: $$\begin{array}{}\text{(14)}& D_{\text{db}}(x)=& w\frac{k-2}{k(k-1)}(k{\partial }_{y}P(y,x)+{\partial }_{y}P(x,y)){|}_{y=x}\text{(15)}& C{S}_{\text{db}}(x^{\ast })=& w\frac{k-2}{k(k-1)}(k{\partial }_y^{2}P(y,x^{\ast })+{\partial }_y^{2}P(x^{\ast },y)+(k+1){\partial }_{y,z}P(y,z)){|}_{z=y=x^{\ast }}\text{(16)}& E{S}_{\text{db}}(x^{\ast })=& C{S}_{\text{db}}(x^{\ast })-w\frac{(k-2{)}^2(k+1)}{k^2(k-1)}{\partial }_{y,z}P(y,z){|}_{z=y=x^{\ast }}.\end{array}$$

Birth-death process

For structured populations with birth-death updating, a parent is first selected from the entire population with a probability proportional to its birthrate, and then its offspring replaces one of the parent's $k$ neighbours, selected uniformly at random. Thus, with birth-death updating, competition occurs at the scale of the entire population, rather than just locally. neighbors of individuals with high birthrates tend to be short-lived, whereas those with neighbors having low birthrates tend to live longer. This results in non-uniform life expectancies, with high turn-over in high payoff regions and low turn-over in low payoff regions. For strong selection the low birthrates can result in almost frozen regions. In general, because of global competition, outcomes of birth-death processes are more similar to outcomes in well-mixed populations than are outcomes of the corresponding death-birth processes (always assuming that selection acts on birth rates).

The setup for pair approximation with birth-death updating is basically the same as in the death-birth updating: the frequency of mutants and mutant pairs increase whenever a mutant reproduces and its offspring replaces a resident neighbour. For a mutant with $j$ mutant neighbours this happens with probability: $$\begin{array}{}\text{(17)}& T_{\text{bd}}^{+}(j)& =(\genfrac{}{}{0ex}{}{k}{j})q_{m|m}^j(1-q_{m|m}{)}^{k-j}\cdot \frac{p_m{b}_m(j)}{\mathcal{B}}\cdot \frac{k-j}{k},\end{array}$$ where $\mathcal{B}=p_m{b}_m(k{q}_{m|m})+(1-p_m)b_r(k{q}_{m|r})$ represents the average birth rate in the population (recall that $q_{m|r}=p_m(1-q_{m|m})/(1-p_m)$), and $b_r(i),b_m(i)$ denote the birth rates of residents and mutants with an (expected) number of $i$ mutant neighbours. The first term in Eq. $\text{(17)}$ represents the probability that the mutant parent has $j$ mutant neighbours, the second term indicates the probability that this mutant is selected for reproduction and the last term the probability that one of its $k-j$ resident neighbours is replaced. Similarly, the frequencies of mutants and mutant pairs decrease if a resident is selected for reproduction and replaces a mutant. For a resident with $j$ mutant neighbours this happens with probability: $$\begin{array}{}\text{(18)}& T_{\text{bd}}^{-}(j)& =(\genfrac{}{}{0ex}{}{k}{j})q_{m|r}^j(1-q_{m|r}{)}^{k-j}\cdot \frac{(1-p_m)b_r(j)}{\mathcal{B}}\cdot \frac{j}{k}.\end{array}$$ All other transitions do not alter the composition of the population. Together $T_{\text{bd}}^{\pm }(j)$ define the rate of change of the frequency of mutants: $$\begin{array}{}\text{(19)}& {\dot{p}}_m& =\sum _{j=0}^k(T_{\text{bd}}^{+}(j)-T_{\text{bd}}^{-}(j)).\end{array}$$ The derivation of the rate of change of $q_{m|m}$ is a bit trickier for birth-death updating. Let us again start by focussing on changes in $p_{mm}$. First, consider a mutant that has successfully replaced a resident neighbor. This resident neighbour has one mutant neighbour (the reproducing individual) as well as an expected $q_{m|r}(k-1)$ further mutants among its remaining $k-1$ neighbours. Hence $p_{mm}$ increases at a rate proportional to $(1+q_{m|r}(k-1))/(Nk/2)$, where $Nk/2$ indicates the normalization, i.e. the total number of (undirected) links in a regular graph of size $N$ and degree $k$. Similarly, if a resident has successfully replaced a mutant neighbour, $p_{mm}$ decreases at a rate $q_{m|m}(k-1)/(Nk/2)$, because the mutant neighbour has itself $q_{m|m}(k-1)$ mutant neighbours (note, one is a resident with certainty, i.e. the reproducing individual). Thus, $$\begin{array}{}\text{(20)}& {\dot{p}}_{mm}& =\sum _{j=0}^k\frac{2}{k}(T_{\text{bd}}^{+}(j)(1+q_{m|r}(k-1))-T_{\text{bd}}^{-}(j)q_{m|m}(k-1)).\end{array}$$ Again, the term $1/N$ in Eqs. $\text{(19)}$ and $\text{(20)}$ has been absorbed through a constant rescaling of time. Similarly, both equations share the common factor $1/\mathcal{B}$, which can be absorbed through a non-linear rescaling of time (because $\mathcal{B}>0$). Neither scaling changes the direction of selection or the location of equilibrium points. Note that the summations in Eqs. $\text{(19)}$ and $\text{(20)}$ can be carried out. Finally, using ${\dot{q}}_{m|m}=({\dot{p}}_{mm}-q_{m|m}{\dot{p}}_m)/p_m$ and $q_{m|r}=(1-q_{m|m})p_m/(1-p_m)$, we obtain: $$\begin{array}{}\text{(21)}& {\dot{p}}_m& =p_m(1-q_{m|m})b_m(q_{m|m}(k-1))-(1-p_m)q_{m|r}{b}_r(1+q_{m|r}(k-1))\text{(22)}& {\dot{q}}_{m|m}& =\frac{1}{k}((1-q_{m|m})(2+2(k-1)q_{m|r}-k{q}_{m|m})b_m((k-1)q_{m|m})+(1-\frac{1}{p_m})q_{m|m}{q}_{m|r}(k-2)b_r(1+(k-1)q_{m|r})),\end{array}$$

Spatial invasion fitness

The invasion fitness of mutants, $f(y,x)$, is the per capita growth rate, ${\dot{p}}_m/p_m$, in the limit $p_m\to 0$: $$\begin{array}{}\text{(23)}& f(x,y)& =(1-q_{m|m}^{\ast })\exp [-wP(x,x)]\times \text{(24)}& & (\exp [\frac{w}{k}((k-(k-1)q_{m|m}^{\ast })P(y,x)+(k-1)q_{m|m}^{\ast }P(y,y))]-\exp [\frac{w}{k}((k-1)P(x,x)+P(x,y))]).\end{array}$$ Here $q_{m|m}^{\ast }$ is the solution of Eq. $\text{(22)}$ in the limit $p_m\to 0$: for $p_m\ll 1$ the time scales of the local (fast) and global (slow) dynamics again separate (see Eqs. $\text{(21)}$, $\text{(22)}$), such that the local densities of mutants can be assumed to be at equilibrium $q_{m|m}^{\ast }$. In the limit $p_m\to 0$ the local dynamics, Eq. $\text{(22)}$, reduces to $$\begin{array}{rl}{\dot{q}}_{m|m}& =-\frac{1}{k}(1-q_{m|m})\exp [-\frac{w}{k}(P(x,y)-P(x,x))]((k-2)q_{m|m}+(k{q}_{m|m}-2)\times \\ \text{(25)}& & \exp [\frac{w}{k}((k-(k-1)q_{m|m})P(y,x)+(k-1)q_{m|m}P(y,y)-P(x,y)-(k-1)P(x,x))]),\end{array}$$ but analytical solutions remain inaccessible. However, since adaptive dynamics is based on the assumption that differences between residents and mutants are small, $|x-y|\ll 1$, we consider a Taylor expansion in $y$ around $x$ of the right-hand-side of Eq. $\text{(25)}$ to obtain the first order approximation of $q_{m|m}^{\ast }$: $$\begin{array}{}\text{(26)}& q_{m|m}^{\ast }=& \frac{1}{k-1}+w(y-x)\frac{k-2}{2(k-1{)}^2}{\partial }_{z}P(z,x){|}_{z=x}+O((y-x{)}^2).\end{array}$$ Note that the zeroth order approximation is the same as for death-birth updating (c.f. Eq. $\text{(13)}$).

Selection gradient, convergent and evolutionary stability

The selection gradient for birth-death updating, $D_{\text{bd}}(x)$, is now obtained by inserting Eq. $\text{(26)}$ into Eq. $\text{(23)}$ and evaluating at $y=x$: $$\begin{array}{}\text{(27)}& D_{\text{bd}}(x)=& w\frac{k-2}{k-1}{\partial }_{y}P(y,x){|}_{y=x}.\end{array}$$ In contrast to death-birth updating spatial structure does not affect the sign of the selection gradient as compared to well-mixed populations (see adaptive dynamics in a nutshell) and hence spatial structure does not affect the existence or location of singular points $x^{\ast }$. The Jacobian of the selection gradient and the Hessian of the fitness at $x^{\ast }$ are: $$\begin{array}{}\text{(28)}& C{S}_{\text{bd}}(x^{\ast })& =w\frac{k-2}{k-1}({\partial }_{y,z}P(y,z)+{\partial }_y^{2}P(y,x^{\ast })){|}_{z=y=x^{\ast }}\text{(29)}& E{S}_{\text{bd}}(x^{\ast })& =w\frac{k-2}{k(k-1)}(2{\partial }_{y,z}P(y,z)+k{\partial }_y^{2}P(y,x^{\ast })){|}_{z=y=x^{\ast }}=\frac{2}{k}C{S}_{\text{bd}}(x^{\ast })+w\frac{(k-2{)}^2}{k(k-1)}{\partial }_y^{2}P(y,x^{\ast }){|}_{z=y=x^{\ast }}\end{array}$$ The conditions for convergence stability are again the same as in well-mixed populations, whereas the criteria for evolutionary stability are generally different.

Evolutionary divergence of coexisting traits

For mutant traits $y$ capable of invading a resident population $x$, i.e. if the invasion fitness $f(y,x)>0$, ecological (frequency) dynamics drives the spread of the mutant trait, possibly replacing the resident. However, if the "resident" is equally capable of invading the "mutant", $f(x,y)>0$, then ecological dynamics drives the two traits to a co-existence equilibrium of two resident traits $x$ and $y$. The direction of selection for further evolution in such populations with two (or more) co-existing traits is again determined by the selection gradients. Each of the coexisting resident traits, say $x_1,\dots ,x_m$, has a corresponding selection gradient $D_i$, $i=1,\dots ,m$, which is derived from the invasion fitness $f(y,x_1,\dots ,x_m)$ of a rare mutant $y$ occurring in the resident $x_i$: $$\begin{array}{}\text{(30)}& D_i(x_1,...,x_m)=\frac{\partial f(y,x_1,...,x_m)}{\partial y}{|}_{y=x_i}\end{array}$$ The invasion fitness of a rare mutant in trait $i$ generally depends on all residents traits, and hence so do the selection gradients $D_i$. Given two resident traits emerging from evolutionary diversification in a monomorphic ancestral population, the most basic question is whether the selection gradients in the two residents tend to increase the phenotypic distance between the coevolving traits, i.e., whether the selection gradients lead to evolutionary divergence.

In structured populations with two residents $x_1$ and $x_2$ an analytical derivation of the invasion fitness of a mutant $y$, $f(y,x_1,x_2)$, is based on the local configuration of the neighbourhoods of the two residents as well as of the mutant, which together determine their respective payoffs and birthrates. Because of the time scale separation between the slow global frequency dynamics and the fast equilibration of local configurations, the neighbourhood of the two residents, $x_1$ and $x_2$, is approximated by the equilibria $q_{x_1|x_1}^{\ast }$ and $q_{x_2|x_2}^{\ast }$, see Eq. $\text{(13)}$. These local configurations determine the frequencies of $x_1$ and $x_2$ at their ecological equilibrium because the condition $p_{x_1{x}_1}+p_{x_1{x}_2}+p_{x_2{x}_1}+p_{x_2{x}_2}=1$ for the pair frequencies expressed in terms of $p_{x_1}$ yields $p_{x_1}{q}_{x_1|x_1}+2p_{x_1}(1-q_{x_1|x_1})+(1-p_{x_1})q_{x_2|x_2}=1$ or $$\begin{array}{}\text{(31)}& p_{x_1}& =\frac{1-q_{x_2|x_2}}{2-q_{x_1|x_1}-q_{x_2|x_2}}.\end{array}$$

If a $x_i$ resident dies, the neighbourhood of a mutant neighbour $y_j=x_j+\mathrm{\Delta }$, where $\mathrm{\Delta }\ll 1$ denotes a small mutational change in the parental trait $x_j$, includes

1. the demised $x_i$,
2. one $y_j$ neighbour, as well as
3. $q_{x_j|x_j}(k-2)$ neighbours with trait $x_j$ and
4. $(1-q_{x_j|x_j})(k-2)$ neighbours with the `other' trait $x_{\overline{j}}$, i.e. if $j=1$ then $\overline{j}=2$ and vice versa.

Note that (2) follows because rare mutants $y_i$ also undergo rapid local equilibration, which again yields $q_{y_i|y_i}^{\ast }=1/(k-1)$ up to leading order, while (3) and (4) assume that the local correlations are otherwise the same as for their $x_j$ parent. The resulting birthrate of the $y_j$ neighbour is given by $$\begin{array}{}\text{(32)}& b_{x_i}(y_j)& =\exp (\frac{w}{k}(P(y_j,x_i)+P(y_j,y_j)+(k-2)(q_{x_j|x_j}P(y_j,x_j)+(1-q_{x_j|x_j})P(y_j,x_{\overline{j}})))).\end{array}$$ If the birthrate of $y_j$ exceeds that of $x_j$ it eventually takes over as a new resident replacing $x_j$. Thus, divergence focusses on fitness differences between $y_j$ and $x_j$. On average, the birthrate of $y_i$ is thus given by $b_{x_1,x_2}(y_i)=p_{x_1}{b}_{x_1}(y_i)+(1-p_{x_1})b_{x_2}(y_i)$ and the divergence becomes $$\begin{array}{}\text{(33)}& (D_1(x_1,x_2),D_2(x_1,x_2))& =\underset{\mathrm{\Delta }\to 0}{lim}(\frac{b_{x_1,x_2}(x_1+\mathrm{\Delta })-b_{x_1,x_2}(x_1)}{\mathrm{\Delta }},\frac{b_{x_1,x_2}(x_2+\mathrm{\Delta })-b_{x_1,x_2}(x_2)}{\mathrm{\Delta }}).\end{array}$$ The resulting gradient field is depicted in the mutual invasibility plots for $\mathrm{\Delta }={10}^{-3}$ in Fig.5 in the main text.