Evolutionary dynamics in the two-locus two-allele model with weak selection.

Two-locus two-allele models are among the most studied models in population genetics. The reason is that they are the simplest models to explore the role of epistasis for a variety of important evolutionary problems, including the maintenance of polymorphism and the evolution of genetic incompatibilities. Many specific types of models have been explored. However, due to the mathematical complexity arising from the fact that epistasis generates linkage disequilibrium, few general insights have emerged. Here, we study a simpler problem by assuming that linkage disequilibrium can be ignored. This is a valid approximation if selection is sufficiently weak relative to recombination. The goal of our paper is to characterize all possible equilibrium structures, or more precisely and general, all robust phase portraits or evolutionary flows arising from this weak-selection dynamics. For general fitness matrices, we have not fully accomplished this goal, because some cases remain undecided. However, for many specific classes of fitness schemes, including additive fitnesses, purely additive-by-additive epistasis, haploid selection, multilinear epistasis, marginal overdominance or underdominance, and the symmetric viability model, we obtain complete characterizations of the possible equilibrium structures and, in several cases, even of all possible phase portraits. A central point in our analysis is the inference of the number and stability of fully polymorphic equilibria from the boundary flow, i.e., from the dynamics at the four marginal single-locus subsystems. The key mathematical ingredient for this is index theory. The specific form of epistasis has both a big influence on the possible boundary flows as well as on the internal equilibrium structure admitted by a given boundary flow.

Introduction

One of the central goals of the pioneers of population genetics was to demonstrate that the inheritance and evolution of continuously varying traits could be explained on the basis of Mendelian genetics (Fisher 1918, 1930). Haldane (1931) and Wright (1935) were apparently the first who formulated explicit dynamical models for the evolution of gene frequencies if selection acts on more than one locus. Under various assumptions about dominance, Haldane considered two loci at each of which a wild type and a deleterious variant segregate. Motivated by empirical examples, he assumed that if both variants occur in the same genotype, they have a selective advantage. Wright investigated a model in which finitely many loci contribute to a quantitative trait that is under quadratic selection toward an intermediate optimum. Both Haldane and Wright assumed that gene frequencies at the loci are probabilistically independent, i.e., they are in linkage equilibrium, and derived the stable equilibrium states as well as other properties of their models.

In an investigation designed to show that selection can lead to tighter linkage, Kimura (1956) derived and studied a full two-locus two-allele model, i.e., one that takes into account linkage disequilibrium. The general (deterministic) two-locus two-allele model for the evolutionary dynamics under selection and recombination was derived and investigated by Lewontin and Kojima (1960). They deduced both general properties as well as properties of special cases, such as additive gene action. Among others, they showed that strong epistasis together with linkage disequilibrium can lead to significantly different outcomes than would occur for independent loci.

Extensive analyses of a special class of fitness patterns, the so-called symmetric viability model (which originated from Wright’s and from Kimura’s work), were performed by Bodmer and Felsenstein (1967) and Feldman and Karlin (1970). The latter authors derived all fifteen possible equilibria and determined their stability for several special cases. Later, Feldman and Liberman (1979) showed that as many as four boundary equilibria and two polymorphic equilibria can be simultaneously stable, and Hastings (1985) demonstrated that up to four stable internal equilibria may coexist. The complexity of this model is also underlined by the finding of Ewens (1968) that there is a gap in the range of recombination rates for which a pair of internal equilibria is stable. A comprehensive review of this model and its extension to multiple loci can be found in Christiansen (1999).

Other important special classes of fitness patterns are the additive and the multiplicative model, in which fitnesses of multilocus genotypes are obtained by adding or multiplying the fitnesses of the constituent single-locus genotypes. In the former, additive epistasis is absent, in the latter, multiplicative epistasis is absent. For additive fitnesses, mean fitness is a Lyapunov function (Ewens 1969), all equilibria are in linkage equilibrium, and, generically, every trajectory converges to an equilibrium point (Karlin and Liberman 1978, 1990; Nagylaki et al. 1999). The multiplicative model is much more complicated, although the linkage-equilibrium manifold is invariant. In this model and away from the linkage-equilibrium manifold, mean fitness may decrease and, for intermediate recombination rates, asymptotically stable equilibria may exist that are in linkage disequilibrium (e.g., Moran 1964, 1968; Bodmer and Felsenstein 1967; Nagylaki 1977; Karlin and Feldman 1978; Hastings 1981a). A detailed review of the theory of two-locus and multilocus models is given in Bürger (2000, Chap. 2).

With general fitnesses, the dynamics in the two-locus two-allele model can be complex. The existence of stable limit cycles has been demonstrated both for the continuous-time model (Akin 1979, 1982) and the discrete-time model (Hastings 1981b; Hofbauer and Iooss 1984). Such complex behavior cannot occur if loci are assumed to be independent, i.e., if linkage equilibrium is imposed. Then the dynamics is gradient-like and mean fitness is a global Lyapunov function (Nagylaki 1989).

However, even the case of two independent, diallelic loci has never been analyzed systematically. Although a Lyapunov function exists, the equilibrium structure can still be quite complicated. For instance, Moran (1963) showed that, apart from degenerate cases, the maximum number of internal (polymorphic) equilibria is five, and up to three can be asymptotically stable.

In this study we perform a systematic analysis of the two-locus two-allele model with constant fitnesses under the assumption of linkage equilibrium. The goal is to determine and classify all possible equilibrium structures and phase portraits (Sects. 3, 4). We assume continuous time for reasons outlined below. We have not fully accomplished our goal, however, we identified all 42 possible (equivalence classes of) boundary flows and 190 potentially possible extended boundary flows, i.e., flows at or close to the boundary. Of these 190 extended boundary flows, the existence of 185 could be proved; the other cases remain undecided. In Section S1 of the Supplementary Information (SI), we present corresponding phase portraits. A large number of extended boundary flows admits not only several non-equivalent phase portraits, but also more than one equilibrium structure, as characterized by the number and stability of boundary and internal equilibria.

We use this general analysis to obtain a detailed classification of equilibrium structures and phase portraits for a number of important special cases that have received considerable attention in the literature. These include the case of marginal overdominance or underdominance (Sect. 5); linear isoclines, which turn out to exist if and only if fitnesses among loci are additive or the only epistatic interactions are additive-by-additive (Sect. 6); multilinear epistasis (Sect. 7); equivalent loci (Sect. 8); and the symmetric viability model (Sect. 9). These results provide considerable insight into the interplay of dominance and epistasis in maintaining genetic polymorphism.

The analysis of this simplified model has immediate implications for the full two-locus two-allele model. This is a consequence of a general theorem by Nagylaki et al. (1999), which applies to multilocus systems. These authors proved under weak technical assumptions that if selection is much weaker than recombination, then after an evolutionarily short period, in which linkage disequilibrium decays to close to zero, the dynamics of the full model (either in discrete or in continuous time) is governed by this weak-selection limit. The model investigated in this paper is the weak-selection limit of the two-locus two-allele model with selection and recombination.

Model

We start with the standard two-locus two-allele model with viability selection and discrete time. Thus, we assume a randomly mating, diploid population with discrete and non-overlapping generations in which viability selection acts on two diallelic, recombining loci. Therefore, gametes are in Hardy–Weinberg proportions. Mutation, random drift, and other evolutionary forces are absent.

Let and be the alleles at locus A, and and those at locus B. Let the frequencies of the four gametes , , , and be denoted by , , , and , respectively, where , and let denote the (constant) viability of an individual with genotype ij. In addition to positing absence of sex effects, i.e., , we posit absence of position effects, i.e., . The recombination probability is denoted by r. The frequencies of alleles and are denoted by and , respectively. Letting be the classical measure of linkage disequilibrium, we obtainUnder the above symmetry assumptions, the fitnesses of genotypes are completely specified by the following matrix:Then evolution of gamete frequencies is given by (Lewontin and Kojima 1960)where , is the marginal fitness of gamete i, and is the mean fitness of the population. This is a dynamical system on the simplex which has received much attention in the literature but is well understood only in special cases (see Sect. 1). For a review consult Chapter 2 in Bürger (2000).

If the assumption of linkage equilibrium, i.e., , is imposed, the dynamics (2.3) simplifies to the following system of difference equations defined on the unit square , where by (2.1) mean fitness is only a function of p and q; cf. Haldane (1931) and Wright (1935, 1942). As shown more generally for multiple multiallelic loci by Nagylaki (1989), is monotone increasing along trajectories of (2.4) and constant only at equilibria. Therefore, if all equilibria are isolated points, every trajectory converges to an equilibrium.

In general, the manifold is not invariant under (2.3). However, assuming weak selection, i.e., setting , rescaling time according to , and letting , the so called weak-selection limit of (2.3) is obtained: Here,is the mean (Malthusian) fitness of the population andthe marginal (Malthusian) fitness of gamete i. We note that the dynamics (2.5) remains unchanged if the same constant is added to every . If every is multiplied by the same positive constant, only a change in time scale results. Therefore, in (2.5) we could substitute for without changing the phase portrait. In particular, two of the nine parameters in the fitness scheme (2.2) could be set to fixed, different values.

Nagylaki et al. (1999) proved (for multiple multiallelic loci) that if is given, selection is sufficiently weak, i.e., is sufficiently small, and if all equilibria of (2.5) are hyperbolic, then every trajectory of (2.3) converges to an equilibrium point on an invariant manifold, , which is contained in an O(s) neighborhood of the linkage-equilibrium manifold . The dynamics on is a small perturbation of the time-s map of the weak-selection limit (2.5), which is gradient-like. In particular, it is easy to show directly for (2.5), but also follows from a result by Nagylaki (1989) for (2.4), that for every and if and only if (p, q) is an equilibrium. Therefore, is a strict Lyapunov function for (2.5).

Because , the equilibria of (2.4) and (2.5) are the same, and so are their stability properties (since ). Therefore, if selection is sufficiently weak and after some (usually short) time has passed (Nagylaki 1993), the dynamics of the full two-locus system (2.3) is closely approximated by the dynamics of the weak-selection limit (2.5). Therefore, each solution of (2.3) converges to an equilibrium point, and this equilibrium point is in an O(s) neighborhood of an equilibrium of (2.5) (Theorem 3.1 in Nagylaki et al. 1999).

Equilibria and their stability

The four monomorphic equilibria of (2.5), corresponding to fixation of gamete i, exist always. They represent the corners of the state space (Fig. 1). The eigenvalues of the Jacobian at are easily calculated and are as follows:An equilibrium is called linearly stable, or a sink, if all eigenvalues (of its Jacobian) have negative real part. If at least one eigenvalue has positive real part, it is linearly unstable. It is called a source if all eigenvalues have positive real part, and it is a saddle if eigenvalues with positive and negative real part occur. Obviously, the corner equilibrium is linearly stable if and only if the fitness of the homozygous genotype ii is higher than each of the two ‘neighboring’ single-locus heterozygous genotypes ij, where j differs from i by a single allele.

Fig. 1

The state space with all possible boundary equilibria. The orientation is such that it corresponds to the fitness matrix in (2.2)

Next there may exist up to four equilibria at which one locus is polymorphic, so-called single-locus polymorphisms (SLPs). They are located on the edges of the state space. We denote the equilibrium on the edge connecting with by (Fig. 1). The coordinates of these edge equilibria are easily calculated and are given by Therefore, the edge equilibrium exists, i.e., is in the interior of the edge, if and only if . Obviously, this is just the well-known condition of either overdominance or underdominance applied to a one-locus boundary system. To distinguish this notion of overdominance (underdominance) from that of marginal overdominance (underdominance) introduced below, we call this overdominance (underdominance) on an edge.

The eigenvalues of the edge equilibria, hence their stability, can also be determined quite straightforwardly, because the Jacobian is in triangular form, so that the eigenvalues appear on the diagonal. One eigenvalue, called internal, determines stability of this equilibrium within its edge. For , it is given byTherefore, as is well known, the SLP is linearly stable (unstable) within its edge if and only if there is overdominance (underdominance) on this edge.

The other eigenvalue, called external, determines stability of the edge equilibrium transversal to the boundary. If the external eigenvalue is negative, the equilibrium is called strictly saturated or externally stable (Hofbauer and Sigmund 1998; Karlin 1980). This has the interpretation that the allele missing on this edge cannot invade the population near this boundary equilibrium. The external eigenvalues are of more complicated form than the internal eigenvalues and are presented in “Appendix A.1”. The distinction between external and internal eigenvalues will be essential for our analysis.

Thus, there are at least four boundary equilibria, but there may be up to eight. We ignore the degenerate cases, in which every point on an edge is an equilibrium. This occurs if and only if , where i and j are gametes differing by one allele.

Finally, we turn to the internal, or fully polymorphic, equilibria. By (2.5), they are the solutions of the two equations and that satisfy and . Equivalently, these equations can be written as where and are the marginal fitnesses of the one-locus genotypes at A and B, respectively. (In the presence of linkage disequilibrium, the expressions in (3.5) and (3.6) need to be normalized; see Ewens and Thomson 1977.) Thus, the internal equilibria are the intersection points of the isoclines and .

The following results were proved by Moran (1963). We will give a slightly different proof of the first statement.

Theorem 3.1

If all equilibria are isolated, then (2.5) has at most five internal equilibria.

Five internal equilibria can be realized, and up to three can be sinks.

Sinks correspond to local maxima of .

Proof

To determine the intersection points of the isoclines (3.4a) and (3.4b), we have to solve the fixed point equation Because the numerators and denominators of f(q) and g(p) are polynomials of degree two or less, numerator and denominator of the rational function f(g(p)) are polynomials of degree four or less. Therefore the intersection points, hence the internal equilibria, are the zeros of a polynomial of degree five or less.

See Moran (1963) and panel 6 in Fig. S1d of the SI, which shows the phase portrait of Moran’s example. (c) is also shown by Moran (1963). It follows immediately from the fact that is a strict Lyapunov function.

To determine the stability of an internal equilibrium, we need the Jacobian.

Lemma 3.2

The Jacobian at an internal equilibrium (p, q) is given bywhere , , and . The eigenvalues of J are real.

A simple calculation shows thatTherefore, (2.5a) yieldsbecause at an internal equilibrium and . The other derivatives in J are calculated in a similar way.

Finally, it is straightforward to compute the discriminant of J, which is the square of the trace minus four times the determinant:The inequality holds because . Therefore, all eigenvalues are real.

For the rest of this paper, we impose the assumption

Therefore, eigenvalues at equilibria are negative or positive, but not zero. This assumptions excludes not only curves of equilibria, but also complete dominance or recessivity of an allele.

The stability of an internal equilibrium is most easily determined by employing the planar Routh–Hurwitz criterion. It states that an equilibrium is (i) a saddle point if , (ii) a sink if and , and (iii) a source if and .

Motivated by Lewontin and Kojima (1960), we say that a locus (e.g., A) exhibits marginal, or induced, overdominance at (p, q) if the inequalitieshold. If both inequality signs are reversed, one obtains marginal, or induced, underdominance. The following result was proved by Kojima (1959). We will give a similar, but more direct, proof.

Corollary 3.3

An internal equilibrium (p, q) is linearly stable, i.e., a sink, if and only if both loci exhibit marginal overdominance at (p, q) and .

By the Routh–Hurwitz criterion, the equilibrium (p, q) is a sink if and only if andTherefore, both and must be negative. Because must hold, (3.4) implies that , , and all have the same sign. This together with an analogous argument for locus B proves that marginal overdominance is necessary. The sufficiency condition follows immediately from (3.12).

Remark 3.4

As pointed out by Kojima (1959), (3.12) is equivalent to the condition that the geometric mean of the dominance variances of each locus exceeds the additive-by-additive variance of the two-locus system.

The proof of Corollary 3.3 shows that, at equilibrium, locus A displays marginal overdominance if and only if (cf. Hastings 1982)

Equilibrium structure and flows

Since the internal equilibria are obtained from the solutions of the quintic polynomial equation (3.7), it is generally difficult to determine their number, position, or stability. Therefore, we first study the possible flows on the boundary of the state space. Subsequently, we extend the flows to a neighborhood of the boundary. Using index theory, we will be able to shed light on the equilibrium structure, i.e., the number and stability properties of equilibria (for a precise definition, see Sect. 4.3), and on the possible phase portraits, i.e., the topological structures of the flow (2.5). Finally, we exclude several potentially possible equilibrium structures and generate phase portraits for most of the remaining cases.

To facilitate the characterization of equilibrium structures and phase portraits, we identify flows that are topologically equivalent or obtained by symmetry operations corresponding to a relabeling of alleles at a locus ( or , or ) or of the loci (A or B). More precisely, we call two systems of the form (2.5), or the corresponding fitness schemes M and in (2.2), or the flows, equivalent if there exists an edge-preserving homeomorphism h of onto itself that maps orbits of (2.5) generated by M onto orbits generated by (preserving the arrow of time), i.e., the phase portraits are topologically equivalent. Here, edge-preserving means that each of the four edges is mapped onto an edge, not necessarily onto itself.

An edge-preserving homeomorphism is a composition of one of the eight symmetry operations of the square with a ‘proper’ or ‘pinned’ homeomorphism of the square, i.e., a homeomorphism that leaves the four corners fixed and maps each edge onto itself. The symmetry group of the square, the dihedral group , consists of four reflections and four rotations (including the identity). Therefore, each equivalence class is invariant under rotations by multiples of , reflections about the diagonal or antidiagonal, and reflections about the middle vertical or horizontal axis of the matrix (2.2), or of the unit square (Fig. 1).

Finally, we call a fitness scheme M [or its induced flow (2.5), or its phase portrait] robust if it has a neighborhood of equivalents in . This is essentially the concept of structural stability, adapted to the selection Eq. (2.5). For the single-locus two-allele model, i.e., on every edge, there are three robust equivalence classes: they correspond to the classical selection patterns of overdominance, underdominance, or intermediate dominance (directional selection).

The goal of our paper is to find all robust (equivalence classes of) phase portraits arising from (2.5). A necessary condition for robustness is that all equilibria are hyperbolic, i.e., condition () holds. A classical characterization of structural stability in two-dimensional systems due to Andronov and Pontryagin (1937) implies:

If ( 2.5 ) satisfies ( ) and there is no saddle connection in the interior of , then this system is robust.

A saddle connection is an orbit which has one saddle as its -limit and an other as its -limit. Most phase portraits shown in Section S1 have no saddle connection. Exceptions occur for the symmetry classes , , and (defined below), where the phase portraits were generated by matrices satisfying the corresponding symmetry condition. In all these cases, however, breaking the symmetry yields phase portraits that are members of the same class. Thus, they are robust in this sense. We now explain how to obtain all phase portraits. This requires three steps.

Flows on the boundary

As shown in Sect. 3, on each edge there is either no or one equilibrium, and each edge equilibrium can be internally stable or unstable. This gives different types of (nondegenerate) flows on the boundary. This type can be easily determined from the selection scheme (2.2), by observing the order relations in each boundary row and boundary column, i.e., overdominance, underdominance, or intermediate dominance. Applying the symmetry operations, their number drops to 42 different boundary flows or, more precisely, boundary-flow equivalence classes. They are displayed in Fig. 2.

Fig. 2

The 42 possible boundary flows, or boundary-flow equivalence classes, for (2.5). The 16 flow-reversal pairs are arranged vertically. For instance, is the reversed boundary flow of , whereas is self inverse. Internal stability is indicated by a solid dot (at a corner) or a cross (at an edge equilibrium). A rhombus indicates an internally unstable edge equilibrium

We recall that a matrix is called centrosymmetric if it is invariant under reflections about the central entry, and bisymmetric if it is symmetric and centrosymmetric. A centrosymmetric fitness matrix (2.2) gives rise to the well-studied symmetric viability model (Sect. 9). A symmetric fitness matrix (2.2) represents a model in which the loci A and B are equivalent (Sect. 8).

The following code was used to label boundary flows: represents a flow of symmetry class with n edge equilibria, m of them internally stable, and k linearly stable corner equilibria. Here, , , or means that at least one flow in this class is generated by a bisymmetric, centrosymmetric, or symmetric matrix, respectively. This does not imply that all matrices generating a flow in such a class have the respective symmetry property. The letter refers to a flow for which the flows at (at least) one pair of opposite edges belong to the same single-locus class and have the same direction if there is intermediate dominance. This is the case if the entries in the two boundary columns (or rows) of the fitness scheme have the same order relation. The letter indicates asymmetry, i.e., the boundary flow cannot be generated by a matrix with one of the properties , , , or . In one case, , the above code is not sufficient to identify the equivalence class; therefore, we used , , .

By reversing all arrows in a flow, the reversed flow is obtained. Equivalently, the signs of all entries in the fitness matrix are reversed. It turns out that 10 of the 42 cases are invariant under flow reversal. In Fig. 2, the other 32 cases are placed such that pairs obtained by reversal are arranged vertically.

Extended boundary flows

Now we consider flows not only on the boundary, but in a sufficiently small neighborhood of the boundary. Here, external stability of the edge equilibria plays a central role (Fig. 3). We need some preparation and recall from below Eq. (3.3) when an equilibrium in our two-dimensional model is strictly saturated. Because we impose the hyperbolicity assumption () throughout, we simply use the term saturated equilibrium.

Fig. 3

The boundary-flow class with positive (left) and negative (right) external eigenvalues at the edge equilibria. In the left panel, the edge equilibria are not saturated. In the right panel, is asymptotically stable and saturated with index , and is a saturated saddle with index . Therefore, in both panels. If the flow in the left panel is reversed, i.e., all arrows are reversed, the flow in the right panel is obtained after a rotation by . Extended boundary flows, where the flow reversal has the same boundary flow and but a different external stability, can occur only for boundary-flow classes that are invariant under flow reversal

Therefore, a corner equilibrium is saturated if and only if both eigenvalues are negative. An edge equilibrium is saturated if and only if the external eigenvalue is negative. An internal equilibrium is, by definition, always saturated. For a hyperbolic equilibrium , the index is defined bywhere k denotes the number of positive eigenvalues of the Jacobian at . It follows, that in a planar system an equilibrium with index is a saddle point; sources and sinks have index . The sum of the indices of all saturated boundary equilibria is called the boundary index sum and is denoted by . If there are no saturated boundary equilibria, then .

In general, is not uniquely determined by the boundary flow because external eigenvalues of edge equilibria may be positive or negative for a given boundary flow (cf. Fig. 3). However, is uniquely determined by the extended boundary flow, which we use as a short hand for the equivalence class (in the above sense) of flows with a given boundary flow (class) together with the signs of the external eigenvalues of all edge equilibria. Different extended boundary flows of the same boundary-flow class may have the same .

Lemma 4.1

For a boundary flow of type , the boundary index sum can assume at most the (integer) values satisfyingEvery boundary flow in Fig. 2 satisfies and .

Proof

We have to sum the indices of the saturated equilibria. Obviously, there are k stable, hence saturated, corner equilibria and at most m saturated, internally stable edge equilibria. This yields the upper bound . Because the number of saturated, internally unstable edge equilibria is at most , the lower bound results.

The second statement follows easily by checking Fig. 2.

By systematic construction of the extended boundary flows for each of the 42 boundary flows, we obtain 200 potential extended boundary flows (see Table S1). It is important to recall that for given more than one extended boundary flow may exist. In the following we show that not all of these 200 extended boundary flows can be realized for (2.5).

Theorem 4.2

Assume (2.5) and ().

The boundary index sum cannot occur.

The boundary index sum cannot occur for any of the boundary flows , , , , or .

The boundary index sum can neither occur for nor if the edge equilibria and are saturated and the others are not.

The boundary index sum cannot occur for if the edge equilibrium is saturated and the others are not.

The proof of (a) is given in “Appendix A.2”. The other statements can be proved in a similar manner.

Remark 4.3

The extended boundary flow described in Theorem 4.2(d) is the reversed flow of the first extended boundary flow in Theorem 4.2(c).

Table S1 informs us that there is exactly one potential extended boundary flow for each of the cases excluded by Theorem 4.2, except for and . Then the two potentially possible extended boundary flows are the reversed flows of each other and thus both are excluded by Theorem 4.2.

We could not exclude the boundary index sum for the cases , , , , and . However, we conjecture that does not occur (see Sect. 4.3).

As a consequence of Theorem 4.2 and statement (b) of this remark, the number of potentially possible extended boundary flows reduces to 190.

It is important to note that statements (c) and (d) of Theorem 4.2 exclude the existence of specific extended boundary flows and do not imply that the boundary index sum does not occur in these cases.

The above theorem constrains the range (4.2) of possible values of to

An important tool for drawing conclusions about internal equilibria is the following index theorem:

Theorem 4.4

(Hofbauer 1990) Assume (2.5) and (). Thenwhere the sum runs over all saturated equilibria.

This theorem has a number of important consequences. The simple proofs are left to the reader.

Corollary 4.5

Assume (2.5) and ().

The number of saturated equilibria is odd;

If is odd (even), the number of internal equilibria is even (odd);

The number of internal equilibria is at least and at most five;

If , there are at least internal saddle points;

If , the total number of sinks and sources is at least ;

The boundary index sum is invariant under flow reversal.

Established equilibrium structures and phase portraits

We show that 185 of the 190 theoretically possible extended boundary flows exist. This task is simplified by taking flow reversals into account.

Importantly, flow reversal neither changes the number of internal equilibria nor their index. Therefore, the boundary index sum is invariant under flow reversal; cf. Corollary 4.5(f). If we identify flow-reversal pairs, it is sufficient to consider boundary flow classes. Here, 16 is the number of representatives of boundary flows in Fig. 2 that are not invariant under flow reversal, and 10 is the number of boundary flows that are invariant. These 16 flows give rise to 79 theoretically possible extended boundary flows (SI, Table S1). The 10 self-inverse boundary flows give rise to 42 extended boundary flows, 22 of which form flow-reversal pairs (Table S1). Thus, up to flow reversal, there are different extended boundary flows. By Theorem 4.2, of these boundary flows cannot exist. Overall, we present phase portraits for extended boundary flows (Section S1 of the SI). Applying flow reversal, this proves existence of different extended boundary flows. All these flows have boundary index sum .

We could not realize potential extended boundary flows (five, including flow reversal), all of them having boundary index sum . Therefore, their existence remains undecided. For the following reasons, we conjecture that they do not exist. (i) Theorem 4.2(b) excludes the other five extended boundary flows with . (ii) A boundary index sum of requires that the number of saddles in the interior is less than the total number of sinks and sources in the interior minus one. Despite extensive search no such equilibrium pattern could be found (see Section S1 of the SI).

Two flows yielding the same extended boundary flow (up to equivalence in the sense defined above) may still differ in their equilibrium structure because they differ in the number of internal equilibria or, if they have the same number of internal equilibria, in the number of sinks, sources, or saddles. In addition, some equilibrium structures can be realized by non-equivalent phase portraits; see below.

We show existence of an extended boundary flow or equilibrium structure by presenting a numerical fitness matrix () yielding a flow on the entire state space, , that has the given extended boundary flow (Section S1 of the SI). We represent each equivalence class of flows on by a phase portrait. In general, the phase portrait is not uniquely determined by the extended boundary flow, because already number and stability of internal equilibria are often not uniquely determined.

We present more than one phase portrait for several extended boundary flows. We do this for cases for which we found more than one equilibrium structure. In addition, we give examples of fitness matrices that yield the same equilibrium structure but non-equivalent phase portraits. One particularly interesting case is with , which admits three extended boundary flows and seven different phase portraits. One of these extended boundary flows has two equilibrium structures, each with two phase portraits (Fig. S2i, panels 1, 2, 3, 4). Another remarkable case is with , which admits four different extended boundary flows (Fig. S2j, panels 2, 3, 4, 5). However, there is also a considerable number of boundary flows for which the value of uniquely determines the extended boundary and also the equilibrium structure and the phase portrait (up to equivalence). See Table S1 and the phase portraits in Section S1 for detailed information.

Permanence

In mathematical modeling of biological systems, the notion of permanence is very important. The dynamical system (2.5) is permanent if there exists a compact subset of such that every solution starting at enters this subset and remains there (consult Hofbauer and Sigmund 1998 for an account on permanence theory). Permanence is equivalent to the boundary being a repeller. It can be shown that it is sufficient that all stationary points on the boundary are repelling.

We want to identify which extended boundary flows can be permanent. First, implies that at least one boundary equilibrium is saturated and, therefore, attracts at least one orbit from the interior. Second, not all extended boundary flows with can be permanent, because there are extended boundary flows for which one equilibrium on the boundary is saturated with index and an other has index . Therefore, an extended boundary flow is permanent if and only if no boundary equilibrium is saturated. Thus, we obtain:

Theorem 4.6

A full flow consistent with a boundary flow of type is permanent if and only if no edge equilibrium is saturated.

Corollary 4.7

Assume (2.5) is permanent. Then

the extended boundary flow is of one of the 14 classes , , , , , , , , , , , , , or , and has no saturated edge equilibrium;

there exists an internal equilibrium that is a sink, for example the global maximum of the mean fitness.

follows from Theorem 4.6 upon comparison with Fig. 2.

follows because is a (strict) Lyapunov function.

We note that each of the extended boundary flows in Corollary 4.7 could give rise to one, three, or five internal equilibria. However, the boundary flow is the only one for which we found a fitness matrix such that (2.5) has five internal equilibria (see Fig. S1d, panel 6).

Continuous isoclines: marginal overdominance or underdominance

In this and the next section we study which equilibrium structures and phase portraits are obtained by imposing specific properties on the isoclines. The structure of the two isoclines is important because the internal equilibria are their intersection points. As special cases, we will encounter well known models.

We start by investigating the equilibrium structure if the isoclines are continuous and map [0, 1] into (0, 1). From (3.11), we recall the definitions of marginal overdominance and underdominance.

Lemma 5.1

The isocline f(q) (g(p)), defined in (3.4), is continuous and maps [0, 1] into (0, 1) if and only if locus A (B) exhibits either marginal overdominance or underdominance for every .

From the definition of the isocline , we obtain that f(q) is continuous if and only if its denominator () does not change sign on [0, 1]. The requirement is then satisfied if and only if both and have the same sign as the denominator. The argument for g(p) is analogous.

This lemma implies that the graph of f(q) intersects the boundary of the state space (precisely) at the equilibria and , and the graph of g(p) does so at and . If both isoclines are continuous and map [0, 1] into (0, 1), we conclude from Lemma 5.1 and Theorem 3.1 that the number of interior equilibria is between one and five, and both bounds can be assumed.

Next we want to derive the possible equilibrium structures if both loci exhibit marginal overdominance everywhere, i.e., at every . Obviously, this is a much stronger assumption than marginal overdominance at an equilibrium, as used in Corollary 3.3. We will need the following lemma:

Lemma 5.2

Locus A exhibits marginal overdominance everywhere if and only ifi.e., overdominance holds on the edges and , and

We observe thatif and only if , , and . Therefore, the conditions in the lemma follow immediately because and can be written in the form of (5.3) and need to be positive.

Theorem 5.3

Assume that both loci exhibit marginal overdominance everywhere. Then the following holds:

The boundary flow is of type .

No boundary equilibrium is saturated, whence and the system is permanent.

There may exist one, three, or five internal equilibria.

If there is a unique internal equilibrium, it is globally asymptotically stable.

If there are three internal equilibria, two are sinks and one is a saddle.

If there are five internal equilibria, three are sinks and two are saddles.

Statement (a) is clear from Lemma 5.2. Combining (3.11) with (A.1) and Theorem 4.6 yields (b). Statement (b) together with Corollary 4.5(b) gives claim (c). Theorem 4.4 yields the desired numbers of saddles for (d), (e), and (f). Because the trace of J (3.2) is negative due to (3.11), the internal equilibria with index are sinks.

Phase portraits with one or three internal equilibria are shown in Fig. S1d, panels 4 and 5. We could neither find an example for case (f) in Theorem 5.3 nor could we exclude it.

Theorem 5.4

Assume that both loci exhibit marginal underdominance everywhere. Then the following holds:

The boundary flow is of type .

All boundary equilibria are saturated, whence .

There may exist one, three, or five internal equilibria.

If there is a unique internal equilibrium, it is a source.

If there are three internal equilibria, two are sources and one is a saddle.

If there are five internal equilibria, three are sources and two are saddles.

Remark 5.5

The possible flows on for Theorem 5.4 are obtained by flow reversal from flows occurring by Theorem of 5.3.

Finally, the following can be shown. We leave the simple proof to the reader.

Theorem 5.6

Assume one locus exhibits marginal overdominance everywhere and the other marginal underdominance everywhere.

The boundary flow is of type .

The two internally stable edge equilibria are saturated, whence .

There exists a unique internal equilibrium, and it is a saddle.

Remark 5.7

It seems interesting that there exist fitness matrices that do not satisfy the assumption of the above theorem and generate a phase portrait with three internal equilibria and the boundary flow . Thus, the assumptions of marginal overdominance and underdominance everywhere do not only constrain the potential boundary flows substantially, but also the phase portraits that can occur for a given boundary flow.

Linear isoclines

Here, we study the dynamics of systems obtained by fitness matrices yielding linear isoclines. This leads to a system of differential equations that has been investigated by Schuster et al. (1981) to study the evolution of two strategies in asymmetric animal contests, for example between two species. It also turns out to be equivalent to a model used by Zhivotovsky and Gavrilets (1992) to study quantitative genetic variation under epistatic selection. A number of well known models emerge as special cases, for instance the additive model or the haploid model.

It will be convenient to reparameterize the fitnesses as follows:Here, are allelic effects, dominance effects, and epistatic effects. According to the remark below (2.7), one additional parameter could be specified (e.g., by setting ), but we refrain from doing so. The general assumption () implies that we exclude complete dominance or recessivity of an allele, i.e., for , and analogous restrictions apply to and .

A straightforward calculation shows that the dynamics (2.5) becomes We note that the right-hand side of (6.2a) is simply the additive effect of allele multiplied by p, and analogously the right-hand side of (6.2b) yields the additive effect of allele .

Remark 6.1

A standard decomposition of the total genetic variance (e.g., Kempthorne 1955) shows that dominance-by-dominance interactions are absent if and only if , additive-by-dominance interactions are absent if and only ifadditive-by-additive interactions are absent if and only if , and dominance is absent if and only if and (6.3) holds.

We infer from (6.2) that both isoclines are linear if and only if (6.3) holds. With this assumption, (6.2) simplifies to

Interestingly, this is a model that was studied independently and in different contexts by Schuster et al. (1981) and Zhivotovsky and Gavrilets (1992). The latter authors introduced an n-locus version of the fitness scheme (6.1) with a generalized version of constraint (6.3), as ‘the simplest generalization of the additive model to include dominance and pairwise additive-by-additive epistasis’. Schuster et al. (1981) identified all possible (non-degenerate) phase portraits for (6.4). Because they were interested in a game-theoretical context, they studied a more general model in which one of the coefficients in (6.4) was substituted by a sixth, independent coefficient. Then limit cycles can occur because the corresponding fitness matrix is not symmetric. Here, we present the results of Schuster et al. (1981) that apply to our model (6.4) in our terminology and complement them. They also follow directly from our results in Sect. 4.

Since the isoclines are linear, there is either no or one internal equilibrium. Therefore, by Corollary 4.5(c), the boundary index sum can assume only the values , or 2. Furthermore, if two edge equilibria exist on opposite edges, they must have the same internal stability. This rules out the following 10 of the 42 boundary flows in Fig. 2: , , , , , , , , , and . In addition, the following boundary flows are easily ruled out: and Finally, the planar Routh–Hurwitz criterion implies that the internal equilibrium resulting from (6.4) isWe summarize the results:

a saddle if ;

a sink if and ; and

a source if and .

Corollary 6.2

Assume (6.4).Therefore, a phase portrait is uniquely determined by its boundary flow and boundary index sum . However, for 12 boundary flows, is not uniquely determined by the boundary flow. For , , , and , can assume all three possible values, 0, 1, and 2.

The following eight boundary flows can have : , , , (in each case the internal equilibrium is globally attracting); , , , (in each case the internal equilibrium is a source). The last four cases are the flow reversals of the first four cases.

The following ten boundary flows and their flow reversals can have : , , , , , (if only is saturated), , (if both edge equilibria are not saturated), (if only is saturated), (if only is saturated). In addition, and have . In all these 22 cases, an internal equilibrium does not exist.

The following five boundary flows and their flow reversals can have : , , , , . In addition, and have . In each of these twelve cases, the internal equilibrium is a saddle.

Zhivotovsky and Gavrilets (1992) derived results about the maintenance of a multilocus polymorphism, i.e., an internal equilibrium, and showed that several models investigated previously, mainly in a quantitative-genetic context, can be obtained as special cases of their model. Among others, they deduced the conditions for the existence and for the stability of a fully polymorphic equilibrium if all loci have equal effects on fitness and selection is sufficiently weak that linkage disequilibrium can be ignored. The assumption of equal effects is equivalent to a symmetric fitness matrix (2.2), i.e., and for every i.

A more detailed study of this symmetric model was performed by Gavrilets (1993). If the fitness matrix is symmetric ( and ), the conditions for the existence and stability of an internal equilibrium become very simple. An internal (nondegenerate) equilibrium exists if and only if . It exists and is linearly stable if and only if . This condition can be satisfied only if . It also implies that there are internally stable equilibria on the edges and . The corresponding boundary flow is with . We will briefly return to this model in Sect. 8, where we treat general symmetric fitness matrices.

Now we turn to two important special cases of the model (6.4).

The additive fitness model

This classical model (Bodmer and Felsenstein 1967; Ewens 1969) assumes that loci contribute additively to fitness. Thus there is no epistasis, but dominance is admitted. It is obtained from (6.1) by assumingIn fact, this model has been analyzed without ignoring linkage disequilibrium, i.e., the full four-gamete system (2.3), and indeed for any number of loci. Ewens (1969) proved for an arbitrary number of multiallelic loci that mean fitness is a Lyapunov function. Moreover, every equilibrium is in linkage equilibrium (for non-zero recombination, as we assume here). For diallelic loci, an internal equilibrium exists if and only if at every locus there is either overdominance or underdominance. An internal equilibrium is unique if it is isolated; it is globally asymptotically stable if there is overdominance on both loci (Karlin and Liberman 1978, 1990). All trajectories converge exponentially to an equilibrium on the linkage-equilibrium manifold (Lyubich 1992). For a more detailed review, see Bürger 2000 (pp. 48–50, 76–78). Therefore, the flows derived by assuming are representative for the full dynamics after a sufficiently long time has passed and also for every trajectory that starts close to .

The isocline simplifies to a horizontal straight line (in a representation as in Fig. 1), and the isocline to a vertical straight line. It follows that locus A (B) exhibits marginal overdominance if and only if it exhibits overdominance on one of the respective edges, i.e., if and only if (); analogously for underdominance. Therefore, in accordance with the above mentioned results, we obtain:

Corollary 6.3

For the additive fitness model the following boundary flows with boundary index sum occur:An internal equilibrium exists and is unique for the boundary flows , , and . It is globally attracting, a source, and a saddle, respectively. Each of the above boundary flows admits only one phase portrait.

with if dominance is intermediate at both loci, i.e., ;

with if dominance is intermediate at one locus and the other locus is overdominant (underdominant);

with if both loci are overdominant (underdominant);

with if one locus is overdominant, the other underdominant.

The haploid model

If selection acts on haploids instead of diploids, fitnesses can be assigned directly to gametes. Denoting the fitness of gamete i by , straightforward calculations show that under the assumption the following dynamics is obtained (cf. Haldane 1931): Comparison with (6.4) reveals that this is obtained from the fitness scheme (6.1) if one sets , (without loss of generality), , , , and assumes (6.3). Therefore, there is additive-by-additive epistasis but no dominance.

The isocline is given by the vertical line , and the isocline is given by the horizontal line . Hence, there are no edge equilibria, but there may exist one internal equilibrium given by these values. The only possible boundary flows are , , and . Their boundary index sum is uniquely determined and , 1, and 2, respectively. It is easy to show that all three cases can be realized. An internal equilibrium exists only in the third case, and it is a saddle. These results are in accordance with results obtained previously for the following more general models.

In game theory, (6.6) is known as the replicator dynamics for partnership games, which is generalized by the replicator dynamics for bimatrix games. There, the coefficients of p and q in (6.6) may differ. Schuster and Sigmund (1981) proved that periodic orbits occur if these coefficients have opposite signs (see Hofbauer and Sigmund 1998, Sections 10 and 11, for a treatment of bimatrix games and the replicator dynamics).

The complete haploid selection model, i.e., without the assumption of linkage equilibrium, was investigated by Felsenstein (1965), Feldman (1971), Rutschman (1994), and Bank et al. (2012). The latter authors proved that there exists at most one internal equilibrium and, if it exists, it is unstable. Then two vertex equilibria are asymptotically stable. This corresponds to the boundary flow . Otherwise, one vertex is globally asymptotically stable and the boundary flow is either or . Bank et al. (2012) also identified the parameter combinations that lead to the respective equilibrium structures.

The multilinear epistasis model

Hansen and Wagner (2001) introduced a model of gene interaction that assumes that the effect of gene substitutions due to changes in the genetic background (the other loci) can be described by a linear transformation. Although their model is formulated in terms of genetic effects on quantitative traits, it can be applied to our context if we consider fitness as the trait. A two-locus version of this model, and its applications to the maintenance of genetic variation, was studied by Hermisson et al. (2003). Following their formulation, the fitness of the two-locus genotype can be written asComparison with the fitness scheme (6.1) shows, after some calculation, that their model is the special case of (6.1) obtained by assumingIt follows that the isoclines and take the form where Therefore, they are in product form. It is not difficult to show that the isoclines of (6.2) are in product form if and only if (7.2) holds. In this case, which we assume below, the dynamics has separated variables.

We defineand observe that () if and only if there is overdominance or underdominance at locus A (B). If and , we call the equilibrium the central equilibrium. We note that and are invariant lines of the dynamics (6.2) with (7.2). As a consequence, saddle connections occur robustly in this class, e.g., for the boundary flow with (Fig. S3h, panel 6).

The isoclines may be composed of up to three straight lines. For (and in a representation as in Fig. 1), these are the horizontal line and the two vertical lines that are given by the solutions of . If , then the two edge equilibria and exist and have the same value, . Thus, edge equilibria occur always in pairs on opposite edges. If in addition , the other two edge equilibria and the central equilibrium exist. Furthermore, the intersection points of and may yield up to four additional internal equilibria. Thus, in total there may be up to five internal equilibria.

The case of five internal equilibria can be realized for overdominance at both loci as well as for underdominance at both loci, but not if one locus exhibits overdominance and the other underdominance (Theorem 7.1 and Table S3). The corresponding boundary flows are with and its flow reversal .

Because opposite edge equilibria occur pairwise, flows of type or cannot occur. The following theorem lists all possible equilibrium structures. Only 16 of the 42 boundary flows in Fig. 2 can occur. In particular, it states the extent to which the equilibrium structure can be inferred from the boundary flow. The proof is given in “Appendix A.3”. Parameter combinations that yield all possible equilibrium structures are given in Table S3.

Theorem 7.1

Assume (6.2) and (7.2).

All possible equilibrium structures in the interior are given in Table 1.

Table 1

Internal equilibrium structures for the multilinear epistasis model

Boundary flow	δ=0	δ=1	δ=2	δ=3	δ=4
E00C_1s	−	0	−	−	−
E00C_1a	−	0	−	−	−
E00C_2b	−	−	1 saddle	−	−
E22C_0c	−	0	−	−	−
E02C_2c	−	0	−	−	−
E22C_0e	−	0^∗	−	−	−
E02C_2e	−	0^∗	−	−	−
E12C_0e	−	0	−	−	−
E12C_1a^′′′	−	−	1 saddle	−	−
E12C_2e	−	0	−	2 saddles	−
E44C_0b	1 sink	−	1 saddle	−	1 source, 4 saddles
E04C_4b	1 source	−	1 saddle	−	1 sink, 4 saddles
E34C_0e	1 sink	−	1 saddle or 1 source and 2 saddles	−	−
E14C_2e	1 source	−	1 saddle or 1 sink and 2 saddles	−	−
E24C_0c	1 sink or 1 source	−	1 saddle	−	−
E24C_1s	−	1 sink or source, 1 saddle	−	−	−

A ‘−’ indicates that this value of does not occur, and ‘0’ indicates that the number of internal equilibria is zero. A comma means ‘and’. Indicates that this case can be realized by matrices with different extended boundary flows; thus there are two different equilibrium structures. For the four boundary flows at the bottom, which admit two different internal equilibrium structures for one value of , each of the internal equilibrium structures is generated by a different, but unique, extended boundary flow

For the top twelve boundary-flow classes in Table 1, the number of asymptotically stable edge equilibria is . For the other four classes, this is in general wrong.

Remark 7.2

Because mean fitness is a strict Lyapunov function (Sect. 2), there is a globally asymptotically stable equilibrium if and only if there is precisely one sink. This occurs in the following cases: and (a corner equilibrium is globally attracting); , , and (an edge equilibrium is globally attracting); and , each with (the central equilibrium is globally attracting); with (in this case the central equilibrium is a sink, hence globally attracting).

We observe from (7.2) that the multilinear model reduces to that with linear isoclines in Sect. 6 if and only if or . Thus, the multilinear model of epistasis coincides with the epistatic model of Zhivotovsky and Gavrilets (1992) if and only if dominance is absent, when it simplifies to the haploid model (Sect. 6.2), or if epistasis is absent, when it reduces to the additive model (Sect. 6.1). It is easy to show that if dominance and epistasis are present, the multilinear model and that of Zhivotovsky and Gavrilets are different and none is a special case of the other. The reason is that in the model of Zhivotovsky and Gavrilets there are no additive-by-dominance or dominance-by-dominance interactions, whereas they occur in the multilinear model if dominance is present (as follows from Eq. 7.2 and Remark 6.1).

Remark 7.3

The following (extended) boundary flows occur in the multilinear epistasis model but not in the model of Zhivotovsky and Gavrilets (1992): , , , , , ; and (each with and ) and with (see Table S2).

The extended boundary flows that occur both in the multilinear model and in the model of Zhivotovsky and Gavrilets are those occuring in the haploid model (no edge equilibria), those occuring in the absence of dominance and epistasis (listed in Corollary 6.3), and and , each with (see Table S2). However, in the latter case the fitness matrices are different.

Finally, we note that the above model was investigated by Bomze et al. (1983) in the context of Mendelian game dynamics of asymmetric contests. Because then the fitness matrix (2.2) is in general not symmetric, periodic solutions may occur. Bomze (1982) provided a classification of all possible phase portraits.

Equal locus effects

If both loci contribute to fitness equally, the matrix in (6.1) is symmetric, i.e.,Under this assumption, which we impose throughout this section, only the nine boundary flows in Fig. 2 occur that have a code ending by an or . However, additional restrictions on the possible phase portraits and equilibrium structures occur as we will demonstrate now.

The symmetry implies that edge equilibria occur pairwise, i.e., () and (). Equilibria of a pair have the same internal and the same external stability. Therefore, boundary flows with an odd number of stable corner equilibria give rise to odd values ; with an even number of stable corner equilibria is even.

Using the symmetry (8.1), we infer from (6.2) that the isoclines take the formwhere for the isocline we have , and for the isocline we have . Therefore, we obtain up to three symmetric internal equilibria. They satisfy and are the solutions of the cubic equation . In addition, there may exist a pair of internal equilibria that are mirror images of each other with respect to the diagonal . Since the internal equilibria are the zeros of the quintic polynomial (3.7), this pair of equilibria is given by the zeros of a quadratic polynomial.

The following theorem determines the number of internal equilibria for all possible boundary flows. We recall from Corollary 4.5 that the number of internal equilibria is odd if and only if is even.

Theorem 8.1

Assume (6.2) and (8.1). Possible boundary flows are of class , where , , and . The boundary index sum is odd for boundary-flow classes with an odd number of stable corner equilibria, and is even otherwise. In addition, the following holds:

For the boundary-flow classes and , each with or , and for with , the number of internal equilibria can be zero (only for ), two, or four.

For the boundary-flow classes and , each with or , the number of internal equilibria can be one (only for ), three, or five.

For the boundary-flow classes (which has ), (which has ), and (which both have or ), and with , and with , the respective maximum number of internal equilibria (four or five) cannot be assumed. Any smaller number of internal equilibria admitted by Corollary 4.5 (1 or 3 if or , 0 or 2 if ) can be assumed.

The proof is given in “Appendix A.4”. Without proof (which is simple), we note that in this symmetric case, can be excluded for . Therefore, the only possible values are and .

Finally, we briefly treat two special cases. For the model of Gavrilets (1993) mentioned in Sect. 6, i.e., Eq. (6.4) with (8.1), the following boundary flows and values can occur: All eight boundary flows from Corollary 6.2 ending with an or occur. For , , , , , and , the boundary index sum is already uniquely determined (, 2, 1, 1, 0, and 0, respectively). For and the case is easily excluded, whence only or are possible.

If the symmetry assumption (8.1) is imposed on the multilinear epistasis model treated in Sect. 7, the following boundary flows and equilibrium structures occur.

() and (); in both cases, the equilibrium structure is unique.

(): both cases (sink and saddle, source and saddle) can be realized.

and : each with and . cannot occur because this would require neighboring edge equilibria to differ in their external stability.

Interestingly, there are no equilibrium structures with two edge equilibria in the multilinear model if loci have equal effects.

The symmetric viability model

As already outlined in the Introduction, the so-called symmetric viability model has received much attention in the literature. One reason is that special cases of it arise naturally when two diallelic loci determine a quantitative character that is under stabilizing selection toward an intermediate optimum (e.g., Wright 1935, 1952; Hastings 1987; Nagylaki 1989; Gavrilets and Hastings 1993; Bürger and Gimelfarb 1999; Willensdorfer and Bürger 2003). For a detailed review, consult Bürger (2000, Chap. 6.2).

It has the property that fitnesses of genotypes are invariant under the simultaneous exchange of with and with . Therefore the resulting fitness matrix (2.2) is centro-symmetric and depends only on four parameters. We shall use the following parametrization (cf. Nagylaki 1989):This is equivalent to the most general form, as first introduced by Bodmer and Felsenstein (1967). The model studied by Lewontin and Kojima (1960) corresponds to the special case . Interestingly, this is also a special case of the multilinear epistasis model treated in Sect. 7 [by setting in (6.1)]. Bodmer and Parson (1962) assumed , which makes the matrix symmetric, in addition to being centro-symmetric. Thus, only boundary flows ending with can be realized. The most general models of stabilizing selection (among those referred to above) require all four parameters in (9.1). However, because the double heterozygote has the highest fitness and fitness of trait values decays symmetrically with distance from the optimum, the four parameters have to satisfy certain inequalities (see Nagylaki 1989).

We analyze the general model (9.1). Instead of , we use the coordinates defined byThis transforms the system (2.5) into which will form the basis for the subsequent analysis.

It is immediate that the origin is an equilibrium and the dynamics is point-symmetric with respect to it. In particular, if is an equilibrium of (9.3), so is , and both have the same eigenvalues because the respective Jacobian matrices are equal. We observe that the isoclines of (9.3) are in product form if , whence a special case of the multilinear epistasis model emerges. The isoclines are linear if .

If the isoclines are not linear, they are given by

Theorem 9.1

In particular, the following holds:

The central, or symmetric, equilibrium O exists always.

If and , the following pairs of unsymmetric internal equilibria may exist:

If , then at most one of these pairs is admissible.

If , then both pairs may be admissible.

If , then O is the only internal equilibrium.

The proof is given in “Appendix A.5.1”.

It is straightforward to compute the Jacobian of (9.3). At an internal equilibrium it simplifies tocf. Lemma 3.2. Evaluation at the central equilibrium O yieldsFrom the Routh–Hurwitz criterion we infer immediately:

Lemma 9.2

Let Then O is

a saddle with index if ,

a sink with index if and ,

a source with index if and .

Note that the sign of can be determined immediately from the fitness scheme (9.1). If the fitness of the double heterozygous genotpye exceeds the arithmetic mean fitness of, for instance, the homozygous genotypes and ( or 2), then .

Remark 9.3

Lemma 9.2 settles the internal equilibrium structure for the case , in which O is the only internal equilibrium.

Let . If , the equilibrium is a saddle, and there exists either no other internal equilibrium (if ) or the curve of equilibria , yielding a degenerate flow. If , then is a line of equilibria, hence again degenerate. If and and , then O is the unique internal equilibrium, and it is a saddle.

For the remainder of this section, we assume and . As a corollary to the above theorem, we obtain

Corollary 9.4

If , then O is the unique internal equilibrium. It is a saddle point and .

We already know from Theorem 9.1 that O is unique. If , then and O is a saddle. Since the index of a saddle point is , Theorem 4.4 shows that the boundary index sum is .

Here is the main result of this section.

Theorem 9.5

The internal equilibrium structures that can occur in the symmetric viability model are given in Table 2.

Table 2

Internal equilibrium configurations for the symmetric viability model

Boundary flow	δ=0		δ=2		δ=4
	ρ<0	ρ>0	ρ<0	ρ>0	ρ<0	ρ>0
E00C_2b	−	−	1 saddle	1 sink or source^b, 2 saddles	−	−
E22C_0c	1 saddle, 2 sinks	1 sink	1 saddle	1 sink or source^a, 2 saddles	−	−
E02C_2c	1 saddle, 2 sources	1 source	1 saddle	1 sink or source^a, 2 saddles	−	−
E44C_0b	1 saddle, 2 sinks	1 sink	1 saddle	1 source, 2 saddles	3 saddles	1 source, 4 saddles
E04C_4b	1 saddle, 2 sources	1 source	1 saddle	1 sink, 2 saddles	3 saddles	1 sink, 4 saddles
E24C_0c	×	1 sink or 1 source^b	1 saddle	1 sink or source^b, 2 saddles	−	−

A ‘−’ indicates that this value of does not occur (see Table S1) and ‘’ indicates that this combination of and cannot occur. A comma means ‘and’. Whether O is a sink or a source needs to be determined from the sign of . The stability of O switches under flow reversal, which (in these classes) does not alter the extended boundary-flow class

The stability of the boundary equilibria can be inferred from the boundary flow type as follows:

The number of asymptotically stable corner equilibria is k.

The number of asymptotically stable edge equilibria is .

The proof is given in “Appendix A.5.2”. Examples of all possible flows can be found in Section S1 of the SI.

From Table 2, we infer immediately that O is the unique internal equilibrium if either and or if and . In fact, the table shows that the number of internal equilibria can always be determined from the extended boundary flow and .

Nagylaki (1989) identified all equilibrium structures and phase portraits for the special case of the symmetric viability model that arises from stabilizing selection toward an optimum situated on the genotypic value of the double heterozygote. He admitted arbitrary functions decaying from the optimum monotonically and symmetrically, and he assumed absence of linkage equilibrium. Thus, his model is a special cases of (9.3) (after transformation between (p, q) and (x, y) coordinates). He proved that only the two phase portraits for can occur and four of the five listed for (if and , then O is a sink in his model).

Discussion

The analysis of the classical two-locus two-allele selection-recombination model is notoriously difficult. As outlined in the Introduction, despite considerable efforts only special cases are well understood, for instance, the models with additive or multiplicative fitnesses. We investigated a simplification of the full two-locus two-allele model (2.3) by assuming that the two loci are independent, i.e., in linkage equilibrium. This simplified model is given by (2.5) and, essentially, goes back to Wright (1942), although special cases had been studied earlier (Haldane 1931; Wright 1935).

The model (2.5) is not only more accessible to mathematical analysis than (2.3), it is also of biological relevance because it has been derived as the weak-selection limit of the full model. Therefore, it provides a good approximation if selection is not too strong and the two loci are unlinked or only weakly linked. Indeed, under the non-degeneracy assumption (), a theorem by Nagylaki et al. (1999) demonstrates that in the full model (2.3) and for sufficiently weak selection relative to recombination, every trajectory converges to an equilibrium point, and every such equilibrium point is a perturbation of an equilibrium point of the weak-selection limit (2.5). Throughout this paper, we imposed condition ().

Importantly, mean fitness is a strict Lyapunov function for (2.5). Hence, every solution converges to an equilibrium point (Sect. 2). This is not always the case in the full model, neither for (2.3) nor its continuous-time analog. For each the existence of stable limit cycles was demonstrated (Akin 1979, 1982; Hastings 1981b; Hofbauer and Iooss 1984) if selection and recombination are of similar strength. For the weak-selection limit (2.5), Moran (1963) showed that, in addition to the eight possible boundary equilibria, there may exist up to five internal equilibria, and three can be simultaneously stable (Theorem 3.1). For the full model (2.3), the maximum number of equilibria is 15, seven of them being internal equilibria. This is an immediate consequence of a result of Altenberg (2010), who proved the conjecture of Feldman and Karlin (1970) that the maximum number of equilibria in a selection-recombination model with n gametes is .

Boundary flows, extended boundary flows, and phase portraits

Although we could not fully accomplish our goal of deriving and classifying all possible equilibrium structures and (equivalence classes of) phase portraits of the weak-selection limit (2.5), we identified the extended boundary flows and determined all equilibrium structures for several important types of fitness patterns. These results yield interesting insights into the role of epistasis and dominance in generating equilibrium structures.

For general fitnesses, we identified all possible boundary flows, i.e., flows on the boundary of the state space (Fig. 1). The four corners correspond to the monomorphic equilibria. The dynamics on the edges correspond to the single-locus dynamics when the other locus is fixed for one or the other allele. There are 42 (topologically) different boundary flows or, more precisely, boundary-flow classes because boundary flows that are obtained by relabeling loci or alleles are identified (Sect. 4, Fig. 2). Of these 42 boundary flows, there are 16 pairs for which a member of a pair is obtained by reversing the flow of the other member. The other 10 boundary-flow classes are self inverse under flow reversal.

These boundary flows are by far not sufficient to describe all possible phase portraits, i.e., the topological structure of the flow on . As an intermediate step, we studied the extended boundary flows (Sect. 4.2). They describe the dynamics in a small neighborhood of the boundary, in particular, the external stability of the edge equilibria. A key ingredient for deriving the possible equilibrium structures and phase portaits on the full state space from a given extended boundary flow is the boundary index sum (defined below Eq. 4.1). The possible values of are constrained by Lemma 4.1, and Theorem 4.2 shows that for specific boundary flows additional values can be excluded. Thus, we are still left with 190 potentially possible extended boundary flows. The most important tools for drawing conclusions about the internal equilibrium structure are Theorem 4.4, which is a special case of a more general index theorem (Hofbauer 1990), and Corollary 4.5. Overall, we showed existence of 185 extended boundary flows by providing a fitness matrix and a phase portrait generating such an extended boundary flow (SI, Figs. S1–S5, Table S1).

For a given boundary-flow class and a value of , there may still exist more than one extended boundary flow yielding this boundary-flow class and this . Moreover an extended boundary flow may be compatible with more than one equilibrium structure, and an equilibrium structure may be generated by non-equivalent phase portraits (see Sect. 4.3 and Table S1, as well as the phase portraits in Section S1). Apart from characterizing all possible equilibrium structures or even all possible phase portraits, also more specific problems remain unresolved. For instance, can a sink and a source in the interior coexist? There are also five extended boundary flows, all with , whose existence we could not exclude. We conjecture that they do not exist (see Sect. 4.3).

An important notion in modeling biological systems is permanence. This is a generalization of the notion of a protected polymorphism, which is mainly used for one-locus two-allele models in (spatially) structured populations. Loosely speaking, permanence means that no type or species will be lost because its frequency will remain above a certain threshold. For our model (2.5), this implies that a permanent system can exhibit only 14 types of extended boundary flows (Corollary 4.7), all having . At least one of these extended boundary flows () can be generated by permanent flows with one, three, or five internal equilibria (Fig. S1d, panels 4, 5, 6). Thus, the phase portrait is not uniquely determined by the extended boundary flow.

For every given fitness matrix generating one of the 14 boundary-flow classes of Corollary 4.7 and parametrized as in (6.1), a sufficiently strong increase of the epistasis parameter yields a permanent flow. This is a consequence of Corollaries 4.7 and A.1, since in the parametrization (6.1), is the only parameter in that is independent of the given boundary flow.

Marginal overdominance or underdominance

In Sect. 5, we assumed that each locus exhibits either marginal overdominance at every point (p, q) or marginal underdominance; see (3.11) and Lewontin and Kojima (1960). This is equivalent to having continuous isoclines (3.4) that map [0, 1] into (0, 1). If there is marginal overdominance at both loci, then is the only possible boundary-flow class and no boundary equilibrium is saturated. Therefore, the extended boundary flow is uniquely determined, , and the system is permanent. This is compatible with having one, three, or five internal equilibria, of which one, two, or three, respectively, are sinks (Theorem 5.3). We could find phase portraits with one or three internal equilibria, but not with five. However, we could not prove that five internal equilibria cannot be realized. Although Moran’s (1963) example (Figure S1d, panel 6) has boundary-flow class and , it does not satisfy the assumptions of Theorem 5.3.

A result analogous to Theorem 5.3 holds if both loci exhibit marginal underdominance because then every flow can be obtained by flow reversal. If one locus exhibits marginal overdominance and the other marginal underdominance, then the boundary-flow class is , there exists a unique internal equilibrium, which is a saddle, and the two internally stable edge equilibria are linearly stable (Theorem 5.6). Hence, the phase portrait is uniquely determined by the boundary-flow class.

Table S2 lists the extended boundary flows and indicates by which of the special fitness patterns they can be generated.

Linear isoclines, or additive-by-additive epistasis

A particularly interesting special class arises if linear isoclines are posited (Sect. 6). With the fitness parameterization (6.1) and the assumption (6.3), the dynamics (2.5) simplifies to (6.4). This type of model has been studied independently, in different generality and in different contexts, by Schuster et al. (1981) and by Zhivotovsky and Gavrilets (1992) and Gavrilets (1993); see Sect. 6 for a more detailed appraisal. In population genetics terms, (2.5) has linear isoclines if and only if epistatic interactions are absent or there are only additive-by-additive interactions. The assumption of linear isoclines rules out 16 of the 42 possible boundary flows and greatly reduces the number of possible phase portraits because there can be at most one internal equilibrium. Also the stability of the internal equilibrium, if it exists, is easily determined. Corollary 6.2 lists all possible flows. In particular, the phase portraits of (6.4) are uniquely determined by the boundary flow and by . For 12 boundary flows, however, may assume more than one value.

Corollary 6.2 shows also that a stable internal equilibrium cannot exist unless there is at least one internally stable edge equilibrium, i.e., there is overdominance in at least one single-locus boundary system. Analogously, an internal source can occur only if there is underdominance in at least one single-locus boundary system.

As discussed in Sect. 6.1 on the additive model, in the absence of epistasis all trajectories of the full model (2.3) converge to the linkage-equilibrium manifold if the recombination rate is positive. Therefore, the phase portraits derived for (2.5) are representative for the full model. In addition, every possible boundary flow determines the phase portrait uniquely (Corollary 6.3).

Another important special case with linear isoclines is the haploid selection model (6.6). Since there is no dominance, the only possible boundary flows are , , and . As in the additive case, the phase portraits of (6.6) are uniquely determined by the boundary flows: in the first two cases there is no internal equilibrium, in the third case there is a saddle. The weak-selection limit (6.6) of the haploid model captures all possible equilibrium structures and phase portraits of the full haploid model with linkage disequilibrium, at least for continuous time (Bank et al. 2012 and Sect. 6.2).

Multilinear epistasis

The so-called multilinear model of epistasis was introduced by Hansen and Wagner (2001). It assumes that the effects of gene substitutions due to changes in the genetic background can be described by a linear transformation (Sect. 7). In this model, all types of epistatic interactions can occur. In the two-locus case, these are additive-by-additive, additive-by-dominance, and dominance-by-dominance interactions (cf. Remark 6.1). Interestingly, the multilinearity assumption (7.1) turns out to be equivalent to assuming that the isoclines are in product form; see (7.3). Therefore, they are composed of vertical or horizontal straight lines. In the absence of dominance, the multilinear model and that of Zhivotovsky and Gavrilets (1992) coincide. In this case, both models reduce formally to the haploid model (6.6). Otherwise, they differ. In fact, up to five internal equilibria may occur in the multilinear model. It is also possible to have overdominance on every edge, i.e., in each marginal one-locus system, but no stable internal equilibrium (Table 1, boundary flow with ). All possible equilibrium structures could be identified (Theorem 7.1). Seven of the 16 possible boundary flows determine the phase portrait uniquely.

Equivalent loci

In Sect. 8, we briefly treat the equilibrium structures generated by symmetric fitness matrices. This is equivalent to assuming that both loci are equivalent, an assumption made in many investigations. Only the nine boundary flows of Fig. 2 ending with an or can occur. The symmetry properties of this model greatly simplify its analysis, so that for every possible boundary flow the possible values of and the possible number of internal equilibria can be determined (Theorem 8.1). Moran’s (1963) example of a flow with three stable internal equilibria is symmetric. Thus, the assumption of symmetry reduces the complexity of the two-locus model only in certain aspects.

The symmetric viability model

As indicated in the Introduction and in Sect. 9, the symmetric viability model may be one of the best studied dynamical systems in population genetics. Although it depends on only four parameters, its complexity seems to preclude a comprehensive mathematical analysis. Even the subclass arising from models of stabilizing selection on a quantitative trait toward an intermediate optimum is well understood only for special cases, such as quadratic or Gaussian stabilizing selection (Bürger 2000, Chap. 6.2; Willensdorfer and Bürger 2003).

The weak-selection limit, Eq.  (9.3) in our parameterization (9.1) and in the transformed coordinates, is still sufficiently complex to admit a wide variety of equilibrium structures and phase portraits. Nevertheless, it is simple enough to admit the identification of all possible equilibrium structures (Theorem 9.5 and Table 2). This requires knowledge of the boundary-flow class (by definition, the types ending with or occur), the boundary index sum, and the sign of two simple compound parameters ( and ).

Theorem 9.5 shows that the maximum number of internal equilibria of (9.3) is five and the maximum number of stable internal equilibria is two. This is different in the full two-locus symmetric viability model with linkage disequilibrium. Then the maximum number of internal equilibria is seven (Feldman and Karlin 1970), and four can be simultaneously stable (Hastings 1985). Interestingly, the same maximum numbers of internal and of stable internal equilibria occur in the special case arising from Gaussian stabilizing selection on a quantitative trait (Willensdorfer and Bürger 2003). In the full two-locus model, three of the seven internal equilibria are so called symmetric equilibria whose gamete frequencies satisfy and . In the weak-selection limit, these symmetry conditions collapse to and or, in the (x, y) coordinates of Sect. 9, to . Thus, in the weak-selection limit, the manifold defined by and collapses to the single point . The four unsymmetric equilibria determined by Feldman and Karlin (1970) for the full model correspond to the four unsymmetric internal equilibria determined by Theorem 9.1.

Multiplicative fitnesses

The reader may have noticed that we did not consider multiplicative fitnesses, although the multiplicative viability model has received great attention in the literature. This model has always been treated in discrete time because then multiplicate fitnesses have an immediate biological meaning. In contrast to the additive viability model, it has the property that the linkage-equilibrium manifold () is invariant under the map (2.3). Also in contrast to the additive model, solutions do not necessarily approach linkage equilibrium, and stable equilibria with may exist. By performing the weak-selection limit, multiplicative fitnesses become additive because . Thus, the weak-selection limits of the additive and the multiplicative model coincide.

Inferring stable two-locus polymorphisms

An interesting, old question is whether the maintenance of a stable two-locus polymorphism requires some form of overdominance at the individual loci. Kojima (1959) showed for independent loci, i.e., our model (2.5), that marginal overdominance (3.11) on both loci at the equilibrium point is necessary for the maintenance of a stable two-locus polymorphism. However, this condition is not sufficient (Corollary 3.3). For the full two-locus model, Hastings (1982) showed that for a small range of recombination rates stable internal equilibria (sinks) may display marginal underdominance. Theorem 5.3 and Fig. S1d (panel 5) show that unstable internal equilibria may exist even if both loci display marginal overdominance on the whole state space.

A simpler and more general question is whether inferences about existence and stability of internal equilibria can be drawn from knowledge of the boundary flow. In general, the answer is negative. Without further restrictions on the fitness scheme, (2.2) or (6.1), there is no boundary flow that determines the phase portrait uniquely. The boundary-flow class is the only one for which an internal equilibrium exists always. All other boundary-flow classes can have and all can be realized without an internal equilibrium (Figs. S2a – S2k). It is even possible to maintain an internal sink if no edge equilibrium exists, i.e., if all four marginal one-locus systems display intermediate dominance. Indeed, all three boundary-flow classes admit an internal sink (Fig. S2a, panels 2 and 4; Fig. S3a, panel 2), but none of our special fitness patterns does so.

In the model with linear isoclines, a stable internal equilibrium can occur only if there is overdominance at least at two edges (boundary flows , , ). Overdominance at all four edges () is sufficient for the existence of an internal sink (Corollary 6.2). For multilinear epistasis, which admits not only additive-by-additive but also additive-by-dominance and dominance-by-dominance interactions (though not in their most general form), an internal sink can exist only if there are four edge equilibria. However, for each of the six resulting boundary flows, there may exist either an internal sink or an internal source (Theorem 7.1). For the symmetric viability model, at least two edge equilibria are necessary for the existence of an internal sink. However, for each of the five resulting boundary flows, there may exist either an internal sink or an internal source (Theorem 9.5). Nevertheless, for the models with linear isoclines (i.e., only additive-by-additive epistasis) and with multilinear epistasis, many boundary-flow classes determine the equilibrium structure uniquely (Corollary 6.2 and Theorem 7.1).

Our treatment of special fitness patterns is not exhaustive. For instance, Hastings and Hom (1990) determined the equilibrium structure of a model with stabilizing selection toward an optimum with arbitrary position by assuming absence of linkage disequilibrium. More recently, Feldman and Puniyani (2006) examined an extension of the Lewontin–Kojima version of the full two-locus symmetric viability model, in which substitutions at locus A do not depend on the background alleles at locus B, but substitutions at locus B depend on the background locus A. The boundary flow classes that arise in this framework are , and

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary material 1 (pdf 9323 KB)

Table 3

Each cell lists the conditions for the validity of the inequalities on top under the respective assumptions on and l on the left

	0<y_+<1	0<y_-<1
r_1>0∧l>0	ρ<0∧ϵ_1<0	−
r_1<0∧l<0	−	ρ<0∧ϵ_1<0
r_1>0>l	(ρ>0∧h_1<0)∨(ρ>0∧h_1>0∧ϵ_1>0)	h_1<0∧ϵ_1<0
r_1<0<l	h_1>0∧ϵ_1>0	(ρ>0∧h_1>0)∨(ρ>0∧h_1<0∧ϵ_1>0)

The symbols and mean ‘and’ and ‘or’