Adaptive topographies and equilibrium selection in an evolutionary game.

It has long been known in the field of population genetics that adaptive topographies, in which population equilibria maximise mean population fitness for a trait regardless of its genetic bases, do not exist. Whether one chooses to model selection acting on a single locus or multiple loci does matter. In evolutionary game theory, analysis of a simple and general game involving distinct roles for the two players has shown that whether strategies are modelled using a single 'locus' or one 'locus' for each role, the stable population equilibria are unchanged and correspond to the fitness-maximising evolutionary stable strategies of the game. This is curious given the aforementioned population genetical results on the importance of the genetic bases of traits. Here we present a dynamical systems analysis of the game with roles detailing how, while the stable equilibria in this game are unchanged by the number of 'loci' modelled, equilibrium selection may differ under the two modelling approaches.

Introduction

Biological fitness typically depends on complicated phenotypes, which depend on multiple genetic loci. This raises an interesting modelling dilemma. At one extreme, one may model selection acting on phenotypes as if they were under simple genetic control at a single haploid locus; this is the ‘phenotypic gambit’ [1] widely used in evolutionary modelling, and referred to as evolutionary game theory when applied to model social interactions [2]. If multiple loci do underly a phenotype then, to make accurate evolutionary predictions, such models should capture inter-locus fitness interactions. Yet, they can be of much greater complexity, having to account for a number of phenotypes that may be exponential in the number of loci involved. At the other extreme, a very simple model may be formulated that considers selection acting independently on frequencies of different alleles at different loci. Such a model would be more tractable, but neglects important quantities such as linkage disequilibrium between loci. Hence, it may give incorrect predictions. An intermediate solution is also possible, through the adoption of multilocus population genetics [3-5].

In this paper, we examine the consequences of these two extreme approaches to modelling a simple general and classical problem: interactions in a social game where the players are assigned distinct roles [6]. Such interactions occur in many contexts, such as those where one individual possesses a territory and the other does not [6], interactions between adult reproductives and helpers [7], or between parents and offspring [8]. Even where payoffs are the same from both individuals’ perspectives, ‘uncorrelated asymmetries’ can lead to different behaviours being stable in the distinct roles, and these have previously been analysed in terms of evolutionary stable strategies at the strategy level [2, 6]. Recently, a new analysis of a social game with roles played between relatives has taken the independent gene-level view, and has shown that this gives the same attracting equilibria as the strategy-level view [9]; thus these equilibria correspond with the fitness-maximising evolutionary stable strategies of the game, regardless of whether they arise from a model using one locus or two. This is intriguing on several fronts. First, modelling selection at the strategy-level is akin to modelling selection acting on a larger number of genes competing for a single locus. Results from population genetics show that ‘adaptive topographies’ that take no account of the underlying genetic-basis of fitness do not exist; moving from modelling a trait using a single locus, to modelling that trait using multiple loci, can lead to population equilibria that do not correspond with population fitness maxima [10]. Second, the dimensions of the phase spaces of the two dynamical systems describing these different modelling levels are different, which means that one should not expect their behaviour to be the same. We show in this paper that a projection of the higher-dimensional system onto the phase space of the other still does not lead to a topologically equivalent system. We present a dynamical systems analysis of both systems in order to elucidate their differences. In particular, we show that they do not have equal numbers of equilibria, but for both models there are always at most two stable coexisting equilibria, and the same stable equilibria exist in both models under the same parametrisations. Despite their identical stable equilibria, equilibrium selection behaviour in the two can differ; seemingly equivalent initial conditions in the two systems can lead them to converge to different stable equilibria.

Analysis

Donation games with roles played between relatives

We consider the donation game with potentially non-additive payoffs as presented in Table 1. Interactions are structured such that there is an ‘uncorrelated asymmetry’ [2]; that is, players occupy distinct behavioural roles, and have different strategies according to the role they occupy. Interactions are further structured such that they occur between genetic relatives [9]. This form of the game provides insight into a particular problem of biological interest, namely selection of non-additive social behaviours between relatives; however, the payoff matrix is equivalent to the original fully general payoff matrix [2], as shown in [9], and hence could capture other biologically-interesting interactions. If we wish to model changes in frequencies, rather than changes in value of a trait that the population is monomorphic for, as in adaptive dynamics [11], then there are two different ways of modelling such a game as a dynamical system. On the one hand, the dynamical system can describe the evolution of the frequencies of all possible strategies, which we shall label genotypes. The set of all possible genotypes for the donation game, denoted , contains four elements, namely, . The genotype dynamics is modelled by the rates of change of the frequencies f • of individuals with these four genotypes, with . The equations are of the form Here, w • is the inclusive fitness of an individual with a given genotype and is the population mean fitness defined as The inclusive fitnesses of individuals with the different genotypes are given in the equivalent neighbour-modulated fitness form [12] by Here, 0 ≤ r ≤ 1 is the average degree of relatedness between interacting individuals within the population, giving the probability that interactants have identical genotypes over-and-above that given by the population frequencies of individuals with these genotypes [13]. In the formulation above, we used the fact that the sum of the frequencies is 1, that is, f DD = 1 − (f CC + f CD + f DC).

Table 1

Payoffs for the non-additive donation game.

	C	D
C	b − c + d	−c
D	b	0

As discussed in [9], the simplest alternative model is based on selection acting independently on the frequencies of sub-strategies for each behavioural role. Here, one describes the rate of change in the frequency of each allele for each role; we shall sometimes refer to this as the gene dynamics view. There are two roles (players) in the donation game and two different alleles, namely, cooperative C and defective D. Since the frequencies in both roles again add up to 1, we only consider the frequencies ϕ C1 and ϕ C2 of cooperative alleles occuring in each of the two roles; the frequencies of defective alleles, and corresponding equations, follow from the equalities ϕ D1 = 1 − ϕ C1 and ϕ D2 = 1 − ϕ C2. The replicator dynamics [14] is then defined as where i, j ∈ {1, 2} and i ≠ j. The inclusive fitnesses of cooperative and defective alleles for each case are now given by Here, ω C1, ω D1, ω C2 and ω D2 are defined in terms of the two distintinct behavioural contexts, which are mutually exclusive [9]. For example, ω C1 is the inclusive fitness of a cooperation allele specifying behaviour for role 1; the first two terms give the expected direct pay-off to the cooperative allele 1 arising from the behaviour it encodes (cooperation), depending on whether allele 2 is of type C or D, while the last term gives the expected pay-off for allele 2 in the social partner, weighted by the relatedness of that partner to the focal individual who occupies role 1.

Formulations (1) and (6) both describe the evolution of four different frequencies, but the dynamical systems are not the same. In particular, note that the state of system (1) is determined by three of the four frequencies, since f CC + f CD + f DC + f DD = 1; this means that the phase space of system (1) has dimension three. The state of system (6), however, is already determined by two of the four frequencies, since ϕ C1 + ϕ D1 = 1 and ϕ C2 + ϕ D2 = 1, and the dimension of its phase space is only two. Due to this difference in dimensions, the two systems cannot be topologically equivalent [15, 16] and one should expect that the behaviour of the two systems is not the same. One may be tempted to believe that the higher-dimensional system (1) implies the behaviour of system (6), because ϕ C1 and ϕ C2 should evolve in the same way as f CC + f CD and f CC + f DC, respectively. However, it is not hard to show that also in this sense the two systems are not topologically equivalent. While the proof is straightforward, it is rather tedious and not very insightful. Therefore, we provide this proof in Supporting Information (S1 Appendix).

Despite this lack of topological equivalence between systems (1) and (6), it might be assumed that the two systems have the same number of stable equilibria and predictions of the long-term behaviour made using either model give the same results. In this paper, we explain in detail that systems (1) and (6) can, in fact, give conflicting predictions for the long-term behaviour. We discuss the equilibria and stability properties for system (1) and compare predictions from system (1) with predictions from system (6) by defining ϕ C1 = f CC + f CD and ϕ C2 = f CC + f DC. In the following section, we first consider system (6).

Analysis of equilibrium states for the gene dynamics model

A detailed analysis of the equilibria for system (6) in their most general form has already been provided in [9]. We present here the analysis as is standard in dynamical systems theory [15, 17] and determine stability properties using the Jacobian matrix; this same approach will also be used for the analysis of system (1). The two-dimensional system (6) can be rewritten explicitly in terms of the two variables ϕ C1 and ϕ C2 and the parameters b, c, d and r as Recall that the dynamics for ϕ D1 and ϕ D2 can readily be deduced from the relationships ϕ D1 = 1 − ϕ C1 and ϕ D2 = 1 − ϕ C2. We focus here on the cases d > 0 and d < 0, where we assume b, c > 0 and 0 < r < 1.

Equilibria are found as solutions that simultaneously satisfy and . The equality holds if Similarly, the equation is satisfied if ϕ C2 = 0, ϕ C2 = 1 or ϕ C1 = e*. Hence, system (6) always has the four equilibria (ϕ C1, ϕ C2) = (0, 0), (1, 0), (0, 1) and (1, 1), and there exists a fifth equilibrium provided 0 < e* < 1, that is, provided b + d > 0. Here, we used the assumptions b > 0 and 0 < r < 1 made earlier in this section; note that b + d > 0 amounts to assuming that any negative non-additive payoffs d, arising when two donators interact, are not large enough to offset the benefits b arising from donation. The conditions on the right-hand side of the double-implication above correspond to conditions A.4 and A.5 derived in [9].

The stability of these equilibria is determined by the eigenvalues of the Jacobian matrix, obtained from linearizing system (6) about the respective equilibria. Let us define The Jacobian matrix at an equilibrium (ϕ C1, ϕ C2) is then defined as Hence, the Jacobian matrices for (ϕ C1, ϕ C2) = (0, 0) and (ϕ C1, ϕ C2) = (1, 1) are diagonal matrices with double eigenvalues e(0) = rb − c and −e(1) = c − d − r(b + d), respectively. Therefore, the origin is stable if and only if e(0) < 0 ⇔ r < c/b. If we again assume b + d > 0, then (1, 1) is stable if and only if −e(1) < 0 ⇔ r > (c − d)/(b + d). We conclude that (0, 0) and (1, 1) are both stable precisely when (e*, e*) exists. This equilibrium (ϕ C1, ϕ C2) = (e*, e*), has the anti-diagonal Jacobian matrix with eigenvalues ±e*[1 − e*](1 + r)d. Hence, (e*, e*) is always a saddle equilibrium. Finally, the Jacobian matrices for (ϕ C1, ϕ C2) = (1, 0) and (ϕ C1, ϕ C2) = (0, 1) are diagonal matrices with both the same eigenvalues, namely, −e(0) = c − rb and e(1) = r(b + d) − (c − d). Therefore, (1, 0) and (0, 1) are sources in the parameter regime where (0, 0) and (1, 1) are both stable. Otherwise, they are typically saddles, because stability of (1, 0) and (0, 1) requires (c − d)/(b + d) > c/b and this can only hold if d < 0; see also [9].

To illustrate the global behaviour of the gene dynamics model (1), let us consider an example of a situation where the equilibrium (e*, e*) exists; as parameters, we choose b = 2, c = 0.5, d = 0.25 > 0, and r = 0.185. Fig. 1 shows the phase portrait for this case in the (ϕ , ϕ )-plane. The (gray) arrows indicate the direction of the flow and clearly show a situation of bistability, with both (0, 0) and (1, 1) (blue dots) attracting nearby points. Note that (1, 0) and (0, 1) (red dots) are both sources, because nearby points all move away from these two equilibria. The basins of the two attracting equilibria are separated by two trajectories of points that flow from the respective two sources to the saddle equilibrium (e*, e*) ≈ (0.4388, 0.4388); all other points near (e*, e*) flow away from (e*, e*). These two special trajectories form the stable manifold of (e*, e*) that acts as a separatrix for the two attractors at (0, 0) and (1, 1).

Fig 1

Phase portrait illustrating bistability of the equilibria (0, 0) and (1, 1) for the gene system (6) with parameters b = 2, c = 0.5, d = 0.25 > 0, and r = 0.185.

The (solid blue) curve through the saddle equilibrium (e*, e*) ≈ (0.4388, 0.4388) is the stable manifold of (e*, e*) that separates the two basins of attraction for (0, 0) and (1, 1).

Behaviour of the gene dynamics model as r varies from 0 to 1

We are primarily interested in how the stability of the equilibrium states change as the degree r of relatedness varies between 0 and 1. We again refer to the results in [9] for comparison. The two cases d > 0 and d < 0 are different and we first consider the case d > 0. If d > 0 then (c − d)/(b + d) < c/b. The existence and stability properties of the equilibria are illustrated in Fig. 2(a), where we plot the equilibria in the (ϕ C1, ϕ C2)-plane versus r. For r < (c − d)/(b + d), the origin is a global attractor (solid blue line), (1, 1) is a repellor (dotted red line), (1, 0) and (0, 1) are saddles (dashed green lines), and (e*, e*) does not exist. When r = (c − d)/(b + d) a bifurcation occurs and the equilibrium (e*, e*) emerges from the (1, 1)-branch; this bifurcation is a transcritical bifurcation, but it is degenerate due to the symmetries of the model and not only the stability of (1, 1), but also of (1, 0) and (0, 1) changes. For (c − d)/(b + d) < r < c/b, both the origin and (1, 1) are attractors (solid blue lines), (1, 0) and (0, 1) are repellors (dotted red lines), and (e*, e*) is a saddle (dashed green line); as illustrated in Fig. 1, the stable manifold of (e*, e*) separates the basins of the two attractors in the system. At r = c/b, another transcritical bifurcation occurs, which is similarly degenerate, where (e*, e*) merges with the origin and again (1, 0) and (0, 1) change stability as well. For r > c/b, the origin is a repellor (dotted red line), (1, 1) is now the global attractor (solid blue line), (1, 0) and (0, 1) are again saddles (dashed green lines), and (e*, e*) no longer exists.

Fig 2

Bifurcation diagrams with 0 < r < 1 of the gene model with d > 0 (a) and d < 0 (b); here, we asssume c > b > 0 are such that 0 < (c − d)/(b + d) < 1 (which is automatically satisfied if d > 0, but not if d < 0).

The stability of the equilibria is indicated by solid (blue), dashed (green) and dotted (red) lines for attractors, saddles and repellors, respectively.

The situation for d < 0 is quite different for the parameter regime where (e*, e*) exists, because it gives rise to bistability of the off-diagonal equilibria (1, 0) and (0, 1). The corresponding bifurcation diagram is shown in Fig. 2(b). The equilibrium (e*, e*) can only exist if b and c are such that 0 < c/b < (c − d)/(b + d) < 1. As before, the origin is the global attractor as long as r < c/b; the equilibrium (1, 1) is a repellor, (1, 0) and (0, 1) are saddles, and (e*, e*) does not exist. There are again two (degenerate) transcritical bifurcations, one at r = c/b and one at r = (c − d)/(b + d), where the equilibrium (e*, e*) merges in opposite order with (0, 0) and (1, 1), respectively. This means that both the origin and (1, 1) are repellors when c/b < r < (c − d)/(b + d), and the equilibria (1, 0) and (0, 1) are the attractors. In this regime, (e*, e*) is again a saddle and its stable manifold separates the basins of (1, 0) and (0, 1); one branch of the stable manifold of (e*, e*) connects to (0, 0) and the other to (1, 1). As for the case d > 0, if r > (c − d)/(b + d) then the origin is a repellor, (1, 1) is the global attractor and (1, 0) and (0, 1) are saddles; the equilibrium (e*, e*) no longer exists.

Let us mention here that the equilibrium (e*, e*) only occurs if d ≠ 0. If d = 0 then (1, 0) and (0, 1) are always saddles, and the origin is the global attractor for r < c/b, with (1, 1) a repellor, while it is a repellor for r > c/b, when (1, 1) is the global attractor. The bifurcation at r = c/b is highly degenerate in this case.

Remark 1 The local stability analysis, along with the location of the stable manifolds of the saddles, completely determines the global behaviour of the system. In particular, the only attractors for system (6) are equilibria, which follows from the Poincaré-Bendixson theorem for planar systems: Any other attractor must be a periodic orbit. Since the lines {ϕ C1 = 0}, {ϕ C1 = 1}, {ϕ C2 = 0}, and {ϕ C2 = 1} are all invariant for system (6), the only possible rotation must occur around an equilibrium in the interior of the square [0, 1] × [0, 1]; however, (e*, e*) is always a saddle, which prevents the existence of a surrounding periodic orbit. See also [16, 18].

Analysis of equilibrium states for the genotype model

As we already mentioned, it might be assumed that stable equilibria of the gene dynamics model (6) should correspond to stable equilibria of the genotype model (1) and, more importantly, in the case of bistability, both systems should have the same predictions for corresponding initial conditions. Therefore, we now compare our findings for the gene dynamics model with a similar equilibrium analysis for the genotype model (1). Recall that the genotype dynamics is modelled as with . In its most general form, this system is fully determined by the dynamics of the frequencies f CC, f CD and f DC, with f DD given by the relationship f CC + f CD + f DC + f DD = 1. Its equilibria satisfy for all combinations . Note that we cannot have all f • = 0, or in other words, the origin (0, 0, 0, 0) is not a solution, because we require f CC + f CD + f DC + f DD = 1. We can show that there also exist no equilibria with all f • ≠ 0 and we have the following Lemma.

Lemma 1 If (f CC, f CD, f DC, f DD) is an equilibrium of (1) with 0 < r < 1 and d ≠ 0 then that is, at least one of its coordinates is zero.

The proof of Lemma 1 is given at the end of the Analysis section.

Equation (11) provides a systematic way to derive all possible equilibria of (1). Furthermore, we can use the number of nonzero coordinates as a guide to ensure we find all of them. This leads to the following complete classification of equilibria of (1).

Theorem 2 Consider the genotype dynamics modelled as system (1) with d ≠ 0 and 0 < r < 1. There are at most eight equilibria, of which at most seven can coexist for a small range of r depending on the choice of the parameters b, c > 0. Based on their numbers of nonzero coordinates, we distinguish three classes:

There are four equilibria with a single nonzero coordinate. These are (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), and (0, 0, 0, 1), which exist for all 0 < r < 1.

There are two equilibria with two nonzero coordinates, namely, , which exists for all 0 < r < 1, and The equilibrium E 23 only exists if c − b < d (for either d > 0 or d < 0), and its bounds of existence are defined by and If d > 0 then E 23 exists for ; this range becomes if d < 0.

There are also two equilibria with three nonzero coordinates, but these must have f CD = f DC; here, either f CC = 0 or f DD = 0. The first equilibrium is If we assume 2b − d > 0, then E 1 can only exist if b > c and its bounds of existence become and As before, E 1 exists only for if d > 0 and only for if d < 0. For the case 2b − d < 0, we must have d > 0 and E 1 can exist only if , with bounds if b > c and if b < c. The only other possible equilibrium is Existence of E 4 requires and the bounds on r become and Again, if d > 0 then E 4 exists only for and the range becomes if d < 0.

The proof of Theorem 2 is given at the end of the Analysis section.

Contradicting predictions from the gene dynamics and genotype models

As an illustration of the global behaviour of the genotype model (1), let us consider the same parameter values used for the gene dynamics model (6) in Fig. 1, namely, b = 2, c = 0.5, d = 0.25 > 0, and r = 0.185. The gene dynamics model (6) has five equilibria for this choice of parameters, which is the largest possible number of equilibria for this two-dimensional model. For the genotype model (1) we find six co-existing equilibria. This model is three dimensional and the (f CC, f CD, f DC)-coordinates of these equilibria are (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1), along with and E 23 ≈ (0.4110, 0, 0). A phase portrait is shown in Fig. 3. Note that the phase space is confined to the tetrahedron bounded by the three coordinate planes f CC = 0, f CD = 0, f DC = 0, and the plane f CC + f CD + f DC = 1. We find that two of the equilibria are stable, namely, (0, 0, 0) and (1, 0, 0). The equilibria (0, 1, 0) and (0, 0, 1) are sources, and and E 23 are saddles. The equilibrium E 23 is the only equlibrium with a two-dimensional stable manifold and it is this surface that separates the basins of the two attracting equilibria in phase space. We computed the stable manifold of E 23 by continuation of a two-point boundary value problem with the package Auto [19, 20]; the formulation of this computational method is described in [21, 22]. Our computation shows that the stable manifold of E 23 is an almost planar surface. It intersects the tetrahedron that defines the phase space of the gene dynamics model (6) in two curves along the sides f = 0 and f = 0, and the closure of this two-dimensional stable manifold includes the straight line f + f = 1 on the side f = 0.

Fig 3

Phase portrait in (f CC, f CD), f DC)-space illustrating bistability of the equilibria (0, 0, 0) and (1, 0, 0) for the genotype model (1) with parameters b = 2, c = 0.5, d = 0.25 > 0, and r = 0.185.

The (blue) surface emanating from the saddle equilibrium E 23 ≈ (0.4110, 0, 0) is the stable manifold of E 23 that separates the two basins of attraction for (0, 0, 0) and (1, 0, 0); compare also Fig. 1.

Despite the fact that there are more equilibria than for the gene dynamics model (6), the phase portrait of the genotype model (1) in Fig. 3 seems rather similar: comparing Fig. 1, there are two attracting equilibria separated by the stable manifold of a saddle equilibrium; since the equilibrium of the genotype dynamics model (1) is contained in the closure of the separatrix, its role in the global dynamics is determined by the stable manifold of E 23. Furthermore, just as for (e*, e*) in model (6), the equilibrium E 23 lies roughly in the middle between the two attracting equilibria. In order to compare the dynamics of these two systems (1) and (6) in more detail, we define the variables f C1: = f CC + f CD and f C2: = f CC + f DC as given by system (1). The two systems could be considered equivalent if any trajectory for system (1) would give rise to a projection onto (f C1, f C2)-coordinates that maps one-to-one to a trajectory for system (6). Note that there is a one-to-one correspondence between the equilibria (0, 0, 0), (1, 0, 0), (0, 1, 0) and (0, 0, 1) of (1) in class (i) and the equilibria (0, 0), (1, 0), (0, 1) and (1, 1) of (6). However, none of the equilibria of (1) map to the equilibrium (e*, e*) of (6). This can have dramatic consequences for the behaviour of the two systems. In particular, it should be possible to choose an initial condition in (f CC, f CD, f DC)-space that lies in the basin of (1, 0, 0), that is, to the right of the stable manifold of E 23, such that its projection onto the (f C1, f C2)-plane lies in the basin of (0, 0). An example to this effect is given in Fig. 4. Here, we consider again the parameters b = 2, c = 0.5, d = 0.25 > 0, and r = 0.185, and choose the initial condition (f CC, f CD, f DC) = (0.25, 0.25, 0.1). Under the flow of (1), this point converges to (1, 0, 0), as indicated by the (brown) curve in Fig. 4(a). However, the projection onto the (f C1, f C2)-plane of this trajectory starts from the initial condition (0.5, 0.35), which lies in the basin of (0, 0) with respect to the flow of (6); the projected trajectory (brown curve) and the trajectory as dictated by (6) (grey curve) are shown in Fig. 4(b).

Fig 4

The initial conditions for a trajectory of the genotype model (1) that converges to the attracting equilibrium (1, 0, 0) in (f CC, f CD, f DC)-space can project onto the two-dimensional phase plane f C1 = f CC + f CD and f C2 = f CC + f DC such that the corresponding trajectory for the gene dynamics model (6) converges to the equilibrium (0, 0).

Panel (a) shows the trajectory for (1) in (f CC, f CD, f DC)-space (brown curve) and panel (b) shows the corresponding projection overlayed on the phase portrait for the gene dynamics model (6); the expected trajectory as defined by (6) is shown in grey.

The example discussed above is a numerical illustration of the following important conjecture.

Proposition 3 Suppose the parameters b, c > 0, d ≠ 0 and 0 < r < 1 are chosen such that the gene dynamics model (6) exhibits bistability between attracting equilibria A 1 and A 2. Consider the basin of attraction of A 1, denoted , and an initial condition , that is, the trajectory through (ϕ C1, ϕ C2) converges to A 1. Let (f CC, f CD, f DC) be an initial condition of the genotype system (1) with f CC + f CD = ϕ C1 and f CC + f DC = ϕ C2; here, we use the same values for b, c, d and r as for (6). It is possible to choose (f CC, f CD, f DC) such that , but the trajectory associated with the flow of (1) does not converge to an attractor that corresponds to A 1.

Proposition 3 is motivated by the fact that there is no candidate equilibrium of (1) that corresponds to the equilibrium (e*, e*) of (6). This equilibrium is important, because its stable manifold acts as a separatrix between the basins of A 1 and A 2. If we assume that the genotype system (1) also exhibits bistability for the chosen parameter values, then there must exist an equilibrium of (1) and corresponding stable manifold that acts as a separatrix in a similar way. The projection of this stable manifold onto the (ϕ C1, ϕ C2)-plane will not be the same as the stable manifold of (e*, e*) and the mismatch causes the possible differences in dynamics of the two systems.

Stability properties of genotype equilibria in class (i) as r varies from 0 to 1

The example illustrated in Fig. 4 does not constitute a proof of Proposition 3, but clearly indicates its validity for a particular choice of parameters. Here, we give a detailed analysis for a class of parameters, where we only consider the case d > 0 and b > c; the case d < 0 can be obtained in a similar fashion.

For d > 0 and b > c, system (6) has coexisting equilibria (ϕ C1, ϕ C2) = (0, 0) and (1, 1) that are both stable in the regime where we used the notation and from Theorem 2. For the equilibrium (ϕ C1, ϕ C2) = (1, 1) is a source instead of a sink, and for the equilibrium (ϕ C1, ϕ C2) = (0, 0) is a source instead of a sink. Let us now investigate the stability properties of the corresponding equilibria (0, 0, 0) and (1, 0, 0) of (1).

The stability of equilibria of (1) is determined by the eigenvalues of the Jacobian matrix Note that f CC + f CD + f DC + f DD = 1 induces a dependency of f DD on all three coordinates. In particular, this means that the partial derivatives must be calculated using the formulation for in terms of f CC, f CD and f DC only. The evaluation at an equilibrium simplifies a lot due to the fact that for any f • ≠ 0.

For the equilibrium (0, 0, 0) almost all terms drop out and we get Hence, the eigenvalues of Jac(0, 0, 0) are on the diagonal and using (2)–(5) with (f CC, f CD, f DC) = (0, 0, 0), we find An equilibrium is stable if and only if all its eigenvalues have negative real part. Since we assume d > 0 and b > c, we find that the stability interval for (0, 0, 0) is Note that E 23 merges with (0, 0, 0) when ; this is a transcritical bifurcation that renders two of the three eigenvalues of (0, 0, 0) unstable. The saddle (0, 0, 0) becomes a source at a second transcritical bifurcation when E 1 merges with it at . We conclude that (0, 0, 0) of (1) destabilises at an r-value below the r-value at which (0, 0) of (6) destabilises.

Let us now consider the equilibrium (1, 0, 0). The Jacobian matrix becomes Analogous to the case for (0, 0, 0), we have . Due to the upper triangular structure, the eigenvalues of Jac(1, 0, 0) are also on the diagonal, with the first one determined by Here, we used the fact that (f CC, f CD, f DC) = (1, 0, 0). Hence, using (2)–(5), we find that the eigenvalues of (1, 0, 0) are given by Therefore, (1, 0, 0) is stable if and only if provided c − d > 0. For r > 0 small enough, (1, 0, 0) is a source; it becomes a saddle with one stable eigenvalue when r increases past , which causes the emergence of E 4 in a transcritical bifurcation. As r increases further, (1, 0, 0) becomes stable in a second transcritical bifurcation; this time, two eigenvalues change sign simultaneously (due to the symmetry f CD = f DC and the bifurcation gives rise to the equilibrium E 23. We conclude that (1, 0, 0) of (1) stabilises at , which lies above the r-value at which (1, 1) of (6) stabilises.

We can utilise this mismatch in stability intervals to illustrate Proposition 3 for a range of r-values with d > 0 and b > c. Consider and let (0, 0) of (6) be the attractor A 1 of Proposition 3. Then almost any initial condition of (6) satisfies the conditions of Proposition 3: almost all initial conditions (f CC, f CD, f DC) of (1) with f CC + f CD = ϕ C1 and f CC + f DC = ϕ C2 will not converge to (0, 0, 0), because (0, 0, 0) is not stable. (The only exceptions are initial conditions that lie on the one-dimensional stable manifold of the saddle (0, 0, 0).) Similarly, for , the equilibrium (1, 1) of (6) is stable, but (1, 0, 0) of (1) is not and Proposition 3 applies.

Remark 2 A complete stability analysis of the equilibria of system (1) is not included in this paper. As for system (6), we believe that the only attractors of system (1) are equilibria, but it is far from straightforward to provide a proof. In contrast to the discussion in [18], system (1) has a three-dimensional phase space and the Poincaré–Bendixson theorem does not apply. Hence, it is possible that an attracting periodic orbit, or even a chaotic attractor exists for system (1). However, the planes defined by {f CC = 0}, {f CD = 0}, {f DC = 0}, and {f CC + f CD + f DC = 0} are all invariant for system (1), as is the ‘diagonal’ plane {f CD = f DC}; each pair of planes intersects in lines that are also invariant, so that we can use the Poincaré–Bendixson theorem to claim that any such periodic orbit must lie in the interior of one of two tetrahedra, namely, the tetrahedron with vertices (f CC, f CD, f DC) = (0, 0, 0), (1, 0, 0), (0, 1, 0), and or the tetrahedron with vertices (f CC, f CD, f DC) = (0, 0, 0), (1, 0, 0), , and (0, 0, 1). System (1) is very similar to the class of equivariant vector fields discussed in [23], for which all trajectories are attracted by the invariant planes, but the global attractor could be a heteroclinic cycle between three saddles. It may be possible to use the approach in [23] for system (1), but we leave the complete analysis as a challenge for future research.

Analysis of equilibrium states for the genotype model (1)

We complete the analysis of the equilibria for the genotype model (1) in their most general form by providing the proofs of Lemma 1 and Theorem 2.

Proof of Lemma 1

Recall that equilibria of (1) satisfy f • = 0 or . Hence, if all f • ≠ 0, we must have for all • ∈ {CC, CD, DC, DD}. This means that so we must have equality of all inclusive fitnesses. Equations (2) and (3) give Here, we used the assumption d ≠ 0. Due to symmetry, even without requiring f CD = f DC, we also have Similarly, (2) and (5) give Using (3) and (5) leads to Similarly, (4) and (5) give Finally, (3) and (4) give It is clear that (16)–(22) can be satisfied simultaneously only if r = 0; for example, w CC = w CD and w DD = w DC require (16) and (21), that is Since 0 < r < 1, this proves the Lemma.

Proof of Theorem 2

Lemma 1 implies that any equilibrium of (1) must have at least one of its coordinates equal to zero. Furthermore, f CC + f CD + f DC + f DD = 1, so the equilibria of (1) can indeed all be classified by the classes listed in Theorem 2. Let us begin with class (i).

Class (i):

The equality f CC + f CD + f DC + f DD = 1 implies that only (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), and (0, 0, 0, 1) are possible candidates for this class. These four points are equilibria of (1) if the equilibrium condition (11) is satisfied for each of their coordinates. Clearly, we only need to check (11) for the single nonzero coordinate f • = 1, for which we require . However, the mean fitness, simply reduces to w • if three of the four frequencies are zero. Hence, (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), and (0, 0, 0, 1) are all equilibria and there are no restrictions on r for their existence.

Class (ii):

This class contains all equilibria with two coordinates equal to zero. Suppose f CC = 0 and consider the case f CD = 0, while f DC, f DD ≠ 0. Then (21) must hold in order to satisfy (11), but f CC + f CD = 0, so there is no (generic) solution. Similarly, if we assume f CD ≠ 0 and f DC = 0, then (20) implies which is not generic. At the special value a two-dimensional continuum of equilibria (0, 0, f DC, f DD) and another two-dimensional continuum of equilibria (0, f CD, 0, f DD) exist that are both not persistent under variations in r. Hence, a generic equilibrium from class (ii) with f CC = 0 must have f DD = 0. Then (22) holds, which gives the candidate . Since w DC = w CD, the mean fitness becomes so is indeed an equilibrium. Note that this equilibrium exists without restrictions on r.

The only other option for equilibria in this class are equilibria with two zero coordinates and f CC ≠ 0. If we assume that the other nonzero coordinate is f CD ≠ 0, then (16) implies because f DC = f DD = 0, so that f CC + f CD = 1. This is again not generic. The same applies to the case f CD = f DD = 0, using (17), and the only remaining candidate is an equilibrium with f CD = f DC = 0. For this case (19) applies and we find The value for f DD follows from the remainder f DD = 1 − f CC. The equality of the inclusive fitnesses for all nonzero frequencies again implies . Hence, the candidate E 23 as given in (12) is indeed an equilibrium.

The existence interval of E 23 is determined by the fact that all coordinates of E 23 must lie in the interval [0, 1]; it suffices to check this for the f CC-coordinate of E 23, since f CC + f DD = 1 then implies 0 ≤ f DD ≤ 1 as well. Let us first consider the case with d > 0; we have: Hence, the existence interval is , which only makes sense if The bounds and defined in Theorem 2 take into account that one could have c − d < 0, in which case .

The case d < 0 is analogous, with ‘≤’ replaced by ‘≥’ and vice versa as soon as the inequality is multiplied by (1 − r)d. Note that we must consider the possibility b + d < 0, but this leads to r < 0, which is not acceptable. Hence, we have b + d > 0 and the bounds and simply swap places. We now need , which leads to the same condition c − b < d as derived in (23) for d > 0; note that c − b < d ⇒ b + d > 0.

Class (iii):

The final class consists of equilibria with three nonzero coordinates. We obtain E 1 given by (13) if we assume f CC = 0. Indeed, for this case, (20), (21) and (22) must hold, which requires and the value for f DD follows from the fact that all frequencies sum up to one. As before, the equality w CD = w DC = w DD implies that is equal to each of these inclusive fitnesses and E 1 is, indeed, an equilibrium.

The existence interval of E 1 is then determined by the values of r for which ; this automatically implies f DD ∈ [0, 1]. Let us first consider the case d > 0. If we assume 2b − d > 0 then we have: These bounds lead to an r-interval if , which holds if b > c; note that the additional condition 0 < r < 1 defines the bounds and given in Theorem 2. If d is large and 2b − d < 0 then E 1 exists for , if b > c, and for , if b < c.

The case d < 0 is again analogous, and we get . The condition leads to the requirement b > c, which automatically ensures that this r-interval is contained in [0, 1].

Let us now consider the possible existence of an equilibrium with f CD = 0 and all other coordinates nonzero. This means that (17), (19) and (21) must hold. Since f CD = 0, equation (21) defines f CC, and combined with (17), this gives here, we used the fact that 0 < r < 1. Hence, there is no admissible equilibrium in class (iii) that satisfies f CD = 0. A similar argument holds for the case with f DC = 0.

The only other possibility is an equilibrium with all nonzero coordinates except for f DD = 0. We must satisfy (16), (17) and (22), which implies Furthermore, f CC + f CD + f DC = 1, so which fixes f CC as well. Hence, E 4 as defined in (14) is an equilibrium of system (1).

As before, we find the existence interval of the equilibrium E 4 using the condition . Let us first consider the case d > 0, which leads to: As for the other equilibria, we must show that these bounds lead to a nontrivial r-interval. We have so E 4 can only exist for d > 0 if ; the bounds and defined in Theorem 2 take into account that 0 < r < 1 as well.

For the case d < 0 we have (1 − r)d < 0 and we find the existence interval , provided the same bound from (24) is satisfied; here we assume b + d > 0 and 2b + 3d > 0. The case b + d < 0 leads to an interval with r < 0, which is not admissible; the case b + d > 0, but 2b + 3d < 0 also requires r < 0. Note that the condition implies and This concludes the investigation of all possible equilibria for system (1). In total, we found the eight equilibria listed in Theorem 2 and there are no other equilibria.

Note that if d > 0 and these values are positive, which means that E 1 and E 4 do not both exist at the same time. Similarly, the opposite equality is satisfied for d < 0, provided these values are positive, so that the same conclusion holds. The equilibrium E 23 can coexist with either E 1 or E 4 in certain regions of parameter space. Therefore, at most seven equilibria coexist.

Discussion

For a simple yet general evolutionary game theory model of social evolution, in which behaviour is conditioned on social role occupied, recent analysis has shown that the stable equilibria under selection are the same regardless of whether one considers selection acting on the entire strategy [2, 6], or acting on independent ‘genes’ for each role [9]. Thus, it may be tempting to assume that, for certain kinds of sufficiently simple models, strategy-level and independent-gene-level approaches yield equivalent answers. Here, we have presented an analysis using the tools of dynamical systems theory, for the simple case of asymmetric non-additive donation games played between relatives. Although we present our analyses in terms of this game, the payoff matrix we use is equivalent to a fully general 2-player payoff matrix [9]. This analysis reveals the following main points: first, the gene dynamics and the genotype dynamics cannot be made topologically equivalent in a dynamical systems sense, since the dimensions of the respective phase spaces are different. It is also not possible to ‘slave’ the dynamics of the higher-dimensional gene dynamics model to the genotype model, because the two models differ in their number of equilibria and in the locations of some of these equilibria. Second, we find additional equilibria for the genotype model to those previously found using techniques from evolutionary game theory, since we find unstable equilibria as well as the previously-discovered stable equilibria. Third, by observing that the unstable equilibria under the gene and the genotype dynamics are different, we show that although the stable equilibria are the same in the two systems, initial conditions often exist in which the population equilibria that result under selection in each system are different. That is, for the same starting population the two different model analyses predict different evolutionary outcomes.

Our results provide an interesting contrast to earlier results from the population genetics literature, that the genetic bases of traits under selection affect whether population equilibria will maximise population mean fitness. These population genetics approaches, briefly reviewed in [24], rest on analysis of sexual models. In particular, analysis of population genetics models shows that the concept of an ‘adaptive landscape’ independent of genetic details is incorrect; while a single-locus model may predict that population equilibria maximise population-mean-fitness, moving to a two-locus model can break the correspondence between equilibria under selection and fitness maxima [10]. Our analysis is an evolutionary game theory one, which is inherently asexual; strategies, or strategy components (‘genes’) reproduce directly. Our analysis of the particular social game presented here also demonstrates a different effect, since here the stable equilibria are the same and correspond with fitness-maximising equilibria, but the selected equilibria can differ. The fact that it is only equilibrium selection, rather than the stable equilibria themselves, that is affected by using the analytically simpler model of this general game may be of interest to some researchers. For applications in which the only result of interest is the stable equilibria under selection, using the simpler two-dimensional systems of equations is safe; for applications in which equilibrium selection does matter, our results describe the areas of disagreement between the predictions of the two-dimensional and three-dimensional representations of selection, enabling an informed choice to be made about which is appropriate to use. The fact that, as described above, the game analysed is actually equivalent to a fully general game matrix for 2-player interactions where players occupy distinct behavioural roles should reinforce this potential interest.

Topological non-equivalence of the two models.

One may be tempted to believe that the higher-dimensional system (1) implies the behaviour of system (6), because ϕ C1 and ϕ C2 should evolve in the same way as f CC + f CD and f CC + f DC, respectively. While it is straightforward to show that the two systems are not topologically equivalent, also not in this sense, the proof is rather tedious and not very insightful. Therefore, we provide this proof in the Supporting Information for completeness.

(PDF)

Click here for additional data file.