A Survey on Migration-Selection Models in Population Genetics

This survey focuses on the most important aspects of the mathematical theory of population genetic models of selection and migration between discrete niches. Such models are most appropriate if the dispersal distance is short compared to the scale at which the environment changes, or if the habitat is fragmented. The general goal of such models is to study the influence of population subdivision and gene flow among subpopulations on the amount and pattern of genetic variation maintained. Only deterministic models are treated. Because space is discrete, they are formulated in terms of systems of nonlinear difference or differential equations. A central topic is the exploration of the equilibrium and stability structure under various assumptions on the patterns of selection and migration. Another important, closely related topic concerns conditions (necessary or sufficient) for fully polymorphic (internal) equilibria. First, the theory of one-locus models with two or multiple alleles is laid out. Then, mostly very recent, developments about multilocus models are presented. Finally, as an application, analysis and results of an explicit two-locus model emerging from speciation theory are highlighted.

populations, and the study of how this genetic variation, shaped by environmental influences, leads to evolutionary change, adaptation, and speciation. Therefore, population genetics provides the basis for understanding the evolutionary processes that have led to the diversity of life we encounter and admire.
Mathematical models and methods have a long history in population genetics, tracing back to Gregor Mendel, who used his education in mathematics and physics to draw his conclusions. Francis Galton and the biometricians, notably Karl Pearson, developed new statistical methods to describe the distribution of trait values in populations and to predict their change between generations. Yule [103], Hardy [44], and Weinberg [99] worked out simple, but important, consequences of the particulate mode of inheritance proposed by Mendel in 1866 that contrasted and challenged the then prevailing blending theory of inheritance. However, it was not before 1918 that the synthesis between genetics and the theory of evolution through natural selection began to take shape through Fisher's [35] work. By the early 1930s, the foundations of modern population genetics had been laid by the work of Ronald A. Fisher, J.B.S. Haldane, and Sewall Wright. They had demonstrated that the theory of evolution by natural selection, proposed by Charles Darwin in 1859, can be justified on the basis of genetics as governed by Mendel's laws. A detailed account of the history of population genetics is given in Provine [87].
In the following, we explain some basic facts and mechanisms that are needed throughout this survey. Mendel's prime achievement was the recognition of the particulate nature of the hereditary determinants, now called genes. Its position along the DNA is called the locus, and a particular sequence there is called an allele. In most higher organisms, genes are present in pairs, one being inherited from the mother, the other from the father. Such organisms are called diploid. The allelic composition is called the genotype, and the set of observable properties derived from the genotype is the phenotype.
Meiosis is the process of formation of reproductive cells, or gametes (in animals, sperm and eggs) from somatic cells. Under Mendelian segregation, each gamete contains precisely one of the two alleles of the diploid somatic cell and each gamete is equally likely to contain either one. The separation of the paired alleles from one another and their distribution to the gametes is called segregation and occurs during meiosis. At mating, two reproductive cells fuse and form a zygote (fertilized egg), which contains the full (diploid) genetic information.
Any heritable change in the genetic material is called a mutation. Mutations are the ultimate source of genetic variability and form the raw material upon which selection acts. Although the term mutation includes changes in chromosome structure and number, the vast majority of genetic variation is caused by changes in the DNA sequence. Such mutations occur in many different ways, for instance as base substitutions, in which one nucleotide is replaced by another, as insertions or deletions of DNA, as inversions of sequences of nucleotides, or as transpositions. For the population-genetic models treated in this text the molecular origin of a mutant is of no relevance because they assume that the relevant alleles are initially present.
During meiosis, different chromosomes assort independently and crossing over between two homologous chromosomes may occur. Consequently, the newly formed gamete contains maternal alleles at one set of loci and paternal alleles at the complementary set. This process is called recombination. Since it leads to random association between alleles at different loci, recombination has the potential to combine favorable alleles of different ancestry in one gamete and to break up combinations of deleterious alleles. These properties are generally considered to confer a substantial evolutionary advantage to sexual species relative to asexuals.
The mating pattern may have a substantial influence on the evolution of gene frequencies. The simplest and most important mode is random mating. This means that matings take place without regard to ancestry or the genotype under consideration. It seems to occur frequently in nature. For example, among humans, matings within a population appear to be random with respect to blood groups and allozyme phenotypes, but are nonrandom with respect to other traits, for example, height.
Selection occurs when individuals of different genotype leave different numbers of progeny because they differ in their probability to survive to reproductive age (viability), in their mating success, or in their average number of produced offspring (fertility). Darwin recognized and documented the central importance of selection as the driving force for adaptation and evolution. Since selection affects the entire genome, its consequences for the genetic composition of a population may be complex. Selection is measured in terms of fitness of individuals, i.e., by the number of progeny contributed to the next generation. There are different measures of fitness, and it consists of several components because selection may act on each stage of the life cycle.
Because many natural populations are geographically structured and selection varies spatially due to heterogeneity in the environment, it is important to study the consequences of spatial structure for the evolution of populations. Dispersal of individuals is usually modeled in one of two alternative ways, either by diffusion in space or by migration between discrete niches, or demes. If the population size is sufficiently large, so that random genetic drift can be ignored, then the first kind of model leads to partial differential equations ( [37], [59]). This is a natural choice if genotype frequencies change continuously along an environmental gradient, as it occurs in a cline [43]. Here we will not be concerned with this wide and fruitful area and refer instead to [7], [9], [67], and [82] for recent developments and references.
Instead, this survey focuses on models of selection and migration between discrete demes. Such models are most appropriate if the dispersal distance is short compared to the scale at which the environment changes, or if the habitat is fragmented. They originated from the work of Haldane [41] and Wright [102]. Most of the existing theory has been devoted to study selection on a single locus in populations with discrete, nonoverlapping generations that mate randomly within demes. However, advances in the theory of multilocus models have been made recently. The general goal is to study the influence of population subdivision and of gene flow among subpopulations on the amount and pattern of genetic variation maintained. The models are formulated in terms of systems of nonlinear difference or differential equations. The core purpose of this survey is the presentation of the methods of their analysis, of the main results, and of their biological implications.
For mathematically oriented introductions to the much broader field of population genetics, we refer to the books of Nagylaki [73], Bürger [13], Ewens [30], and Wakeley [98]. The two latter texts treat stochastic models in detail, an important and topical area ignored in this survey. As an introduction to evolutionary genetics, we recommend Charlesworth and Charlesworth [21].
2. Selection on a multiallelic locus. Darwinian evolution is based on selection and inheritance. In this section, we summarize the essential properties of simple selection models. It will prepare the ground for the subsequent study of the joint action of spatially varying selection and migration. Proofs and a detailed treatment may be found in Chapter I of [13]. Our focus is on the evolution of the genetic composition of the population, but not on its size. Therefore, we always deal with relative frequencies of genes or genotypes within a given population.
Unless stated otherwise, we consider a population with discrete, nonoverlapping generations, such as annual plants or insects. We assume two sexes that need not be distinguished because gene or genotype frequencies are the same in both sexes (as is always the case in monoecious species). Individuals mate at random with respect to the locus under consideration, i.e., in proportion to their frequency. We also suppose that the population is large enough that gene and genotype frequencies can be treated as deterministic, and relative frequency can be identified with probability. Then the evolution of gene or genotype frequencies can be described by difference or recurrence equations. These assumptions reflect an idealized situation which will model evolution at many loci in many populations or species, but which is by no means universal.
2.1. The Hardy-Weinberg law. With the blending theory of inheritance variation in a population declines rapidly, and this was one of the arguments against Darwin's theory of evolution. With Mendelian inheritance there is no such dilution of variation, as was shown independently by the famous British mathematician Hardy [44] and, in much greater generality, by the German physician Weinberg ([99], [100]).
We consider a single locus with I possible alleles A i and write I = {1, . . . , I} for the set of all alleles. We denote the frequency of the ordered genotype A i A j by P ij , so that the frequency of the unordered genotype A i A j is P ij + P ji = 2P ij . Subscripts i and j always refer to alleles. Then the frequency of allele A i in the population is After one generation of random mating the zygotic proportions satisfy 2 P ij = p i p j for every i and j .
A mathematically trivial, but biologically very important, consequence is that (in the absence of other forces) gene frequencies remain constant across generations, i.e., p i = p i for every i .
(2.1) In other words, in a (sufficiently large) randomly mating population reproduction does not change allele frequencies. A population is said to be in Hardy-Weinberg equilibrium if P ij = p i p j .
(2.2) In a (sufficiently large) randomly mating population, this relation is always satisfied among zygotes.
Evolutionary mechanisms such as selection, migration, mutation, or random genetic drift distort Hardy-Weinberg proportions, but Mendelian inheritance restores them if mating is random.

2.2.
Evolutionary dynamics under selection. Selection occurs when genotypes in a population differ in their fitnesses, i.e., in their viability, mating success, or fertility and, therefore, leave different numbers of progeny. The basic mathematical models of selection were developed and investigated in the 1920s and early 1930s by Fisher [36], Wright [102], and Haldane [42].
We will be concerned with the evolutionary consequences of selection caused by differential viabilities, which leads to simpler models than (general) fertility selection (e.g. [50], [73]). Suppose that at an autosomal locus the alleles A 1 , . . . , A I occur.
We count individuals at the zygote stage and denote the (relative) frequency of the ordered genotype A i A j by P ij (= P ji ). Since mating is at random, the genotype frequencies P ij are in Hardy-Weinberg proportions. We assume that the fitness (viability) w ij of an A i A j individual is nonnegative and constant, i.e., independent of time, population size, or genotype frequencies. In addition, we suppose w ij = w ji , as is usually the case. Then the frequency of A i A j genotypes among adults that have survived selection is where we have used (2.2). Here, is the mean fitness of the population and is the marginal fitness of allele A i . Both are functions of p = (p 1 , . . . , p I ) . 3 Therefore, the frequency of A i after selection is (2.5) Because of random mating, the allele frequency p i among zygotes of the next generation is also p * i (2.1), so that allele frequencies evolve according to the selection dynamics p i = p i w ī w , i ∈ I . and describes the evolution of allele frequencies at a single autosomal locus in a diploid population. We view the selection dynamics (2.6) as a (discrete) dynamical system on the simplex S I = p = (p 1 , . . . , p I ) ∈ R I : p i ≥ 0 for every i ∈ I , Although selection destroys Hardy-Weinberg proportions, random mating re-establishes them. Therefore, (2.6) is sufficient to study the evolutionary dynamics.
The right-hand side of (2.6) remains unchanged if every w ij is multiplied by the same constant. This is very useful because it allows to rescale the fitness parameters according to convenience (also their number is reduced by one). Therefore, we will usually consider relative fitnesses and not absolute fitnesses.
Fitnesses are said to be multiplicative if constants v i exist such that for every i, j. Then w i = v iv , wherev = i v i p i , andw =v 2 . Therefore, (2.6) simplifies to which can be solved explicitly because it is equivalent to the linear system It is easy to show that (2.9) also describes the dynamics of a haploid population if the fitness v i is assigned to allele A i . Fitnesses are said to be additive if constants v i exist such that for every i, j. Then w i = v i +v, wherev = i v i p i , andw = 2v. Although this assumption is important (it means absence of dominance; see Sect. 2.3), it does not yield an explicit solution of the selection dynamics.
Example 2.1. Selection is very efficient. We assume (2.8). Then the solution of (2.9) is (2.11) Suppose that there are only two alleles, A 1 and A 2 . If A 1 is the wild type and A 2 is a new beneficial mutant, we may set (without loss of generality!) v 1 = 1 and v 2 = 1 + s. Then we obtain from (2.11): (2.12) Thus, A 2 increases exponentially relative to A 1 . For instance, if s = 0.5, then after 10 generations the frequency of A 2 has increased by a factor of (1 + s) t = 1.5 10 ≈ 57.7 relative to A 1 . If s = 0.05 and t = 100, this factor is (1 + s) t = 1.05 100 ≈ 131.5. Therefore, slight fitness differences may have a big long-term effect, in particular, since 100 generations are short on an evolutionary time scale.

REINHARD BÜRGER
The statement (2.13) is closely related to Fisher's Fundamental Theorem of Natural Selection, which Fisher (1930) formulated as follows: "The rate of increase in fitness of any organism at any time is equal to its genetic variance in fitness at that time." For recent discussion, see [31] and [18]. In mathematical terms, (2.13) shows thatw is a Lyapunov function. This has a number of important consequences. For instance, complex dynamical behavior such as limit cycles or chaos can be excluded. From (2.6) it is obvious that the equilibria are precisely the solutions of (2.14) We call an equilibrium internal, or fully polymorphic, if p i > 0 for every i (all alleles are present). The I equilibria defined by p i = 1 are called monomorphic because only one allele is present.
The following result summarizes a number of important properties of the selection dynamics. Proofs and references to the original literature may be found in [13], Chap. I.9; see also [69], Chap. 9.
Theorem 2.2. 1. If an isolated internal equilibrium exists, then it is uniquely determined.
2.p is an equilibrium if and only ifp is a critical point of the restriction of mean fitnessw(p) to the minimal subsimplex of S I that contains the positive components ofp.
3. If the number of equilibria is finite, then it is bounded above by 2 I − 1.
4. An internal equilibrium is asymptotically stable if and only if it is an isolated local maximum ofw. Moreover, it is isolated if and only if it is hyperbolic (i.e., the Jacobian has no eigenvalues of modulus 1).
5. An equilibrium point is stable if and only if it is a local, not necessarily isolated, maximum ofw.
6. If an asymptotically stable internal equilibrium exists, then every orbit starting in the interior of S I converges to that equilibrium.
7. If an internal equilibrium exists, it is stable if and only if, counting multiplicities, the fitness matrix W = (w ij ) has exactly one positive eigenvalue.
8. If the matrix W has i positive eigenvalues, at least (i − 1) alleles will be absent at a stable equilibrium.
9. Every orbit converges to one of the equilibrium points (even if they are not isolated).

2.3.
Two alleles and the role of dominance. For the purpose of illustration, we work out the special case of two alleles. We write p and 1 − p instead of p 1 and p 2 . Further, we use relative fitnesses and assume

Schematic selection dynamics with two alleles
It is easily verified that the allele-frequency change from one generation to the next can be written as There exists an internal equilibrium if and only if h < 0 (overdominance) or h > 1 (underdominance). It is given bŷ (2.16) shows that ∆p > 0 if 0 < p < 1, hence p = 1 is globally asymptotically stable.
If h > 1, then the monomorphic equilibria p = 0 and p = 1 each are asymptotically stable andp is unstable. The three possible convergence patterns are shown in Figure 2.1. Figure 2.2 demonstrates that the degree of (intermediate) dominance strongly affects the rate of spread of an advantageous allele. The initial frequency is p 0 = 0.005 and the selective advantage is s = 0.05. If the advantageous allele is recessive, its initial rate of increase is vanishingly small because the frequency p 2 of homozygotes is extremely low when p is small. However, only homozygotes are 'visible' to selection.
2.4. The continuous-time selection model. Most higher animal species have overlapping generations because birth and death occur continuously in time. This, however, may lead to substantial complications if one wishes to derive a continuoustime model from biological principles. By contrast, discrete-time models can frequently be derived straightforwardly from simple biological assumptions. If evolutionary forces are weak, a continuous-time version can usually be obtained as an approximation to the discrete-time model.
A rigorous derivation of the differential equations that describe gene-frequency change under selection in a diploid population with overlapping generations is a formidable task and requires a complex model involving age structure (see [73], Chap. 4.10). Here, we simply state the system of differential equations and justify it in an alternative way.
In a continuous-time model, the (Malthusian) fitness m ij of a genotype A i A j is defined as its birth rate minus its death rate. Then the marginal fitness of allele A i is the mean fitness of the population is and the dynamics of allele frequencies becomeṡ This is the analogue of the discrete-time selection dynamics (2.6). Its state space is again the simplex S I . The equilibria are obtained from the conditionṗ i = 0 for every i. The dynamics (2.19) remains unchanged if the same constant is added to every m ij . We note that (2.19) is a so-called replicator equation [51]. If we set w ij = 1 + sm ij for every i, j ∈ I , (2.20) where s > 0 is (sufficiently) small, the difference equation (2.6) and the differential equation (2.19) have the same equilibria. This is obvious upon noting that (2.20) implies w i = 1 + sm i andw = 1 + sm.
Following [73], p. 99, we approximate the discrete model (2.6) by the continuous model (2.19) under the assumption of weak selection, i.e., small s in (2.20). We rescale time according to t = τ /s , where denotes the closest smaller integer. Then s may be interpreted as generation length and, for p i (t) satisfying the difference equation (2.6), we write π i (τ ) = p i (t). Then we obtain formally From (2.6) and (2.20), we obtain . We note that (2.6) is essentially the Euler scheme for (2.19).
The exact continuous-time model reduces to (2.19) only if the mathematically inconsistent assumption is imposed that Hardy-Weinberg proportions apply for every t which is generally not true. Under weak selection, however, deviations from Hardy-Weinberg decay to order O(s) after a short period of time [73].
One of the advantages of models in continuous time is that they lead to differential equations, and usually these are easier to analyze because the formalism of calculus is available. An example for this is that, in continuous time, (2.13) simplifies tȱ which is much easier to prove than (2.13): Remark 2.3. The allele-frequency dynamics (2.19) can be written as a (generalized) gradient system ( [96], [91]): Here, G p = (g ij ) is a quadratic (covariance) matrix, where 1 2 if k = l and k = i or l = i , 0 otherwise .

(2.24)
An equivalent formulation is the following covariance form: where m is interpreted as the random variable m(A k A l ) = m kl (Li 1967). It holds under much more general circumstances ( [85], [62]).
3. The general migration-selection model. We assume a population of diploid organisms with discrete, nonoverlapping generations. This population is subdivided into Γ demes (niches). Viability selection acts within each deme and is followed by adult migration (dispersal). After migration random mating occurs within each deme. We assume that the genotype frequencies are the same in both sexes (e.g., because the population is monoecious). We also assume that, in every deme, the population is so large that gene and genotype frequencies may be treated as deterministic, i.e., we ignore random genetic drift.
3.1. The recurrence equations. As before, we consider a single locus with I alleles A i (i ∈ I). Throughout, we use letters i, j to denote alleles, and greek letters α, β to denote demes. We write G = {1, . . . , Γ} for the set of all demes. The presentation below is based on Chapter 6.2 of [73]. We denote the frequency of allele A i in deme α by p i,α . Therefore, we have for every α ∈ G. Because selection may vary among demes, the fitness (viability) w ij,α of an A i A j individual in deme α may depend on α. The marginal fitness of allele A i in deme α and the mean fitness of the population in deme α are respectively. Next, we describe migration. Letm αβ denote the probability that an individual in deme α migrates to deme β, and let m αβ denote the probability that an (adult) individual in deme α immigrated from deme β. The Γ × Γ matrices M = (m αβ ) and M = (m αβ ) (3.3) are called the forward and backward migration matrices, respectively. Both matrices are stochastic, i.e., they are nonnegative and satisfy βm αβ = 1 and β m αβ = 1 for every α . (3.4) Given the backward migration matrix and the fact that random mating within each demes does not change the allele frequencies, the allele frequencies in the next generation are describes the change due to selection alone; cf. (2.6). These recurrence equations define a dynamical system on the Γ-fold Cartesian product S Γ I of the simplex S I . The investigation of this dynamical system, along with biological motivation and interpretation of results, is one of the main purposes of this survey.
The difference equations (3.5) require that the backward migration rates are known. In the following, we derive their relation to the forward migration rates and discuss conditions when selection or migration do not change the deme proportions.
3.2. The relation between forward and backward migration rates. To derive this relation, we describe the life cycle explicitly. It starts with zygotes on which selection acts (possibly including population regulation). After selection adults migrate and usually there is population regulation after migration (for instance because the number of nesting places is limited). By assumption, population regulation does not change genotype frequencies. Finally, there is random mating and reproduction, which neither changes gene frequencies (Section 2.1) nor deme proportions. The respective proportions of zygotes, pre-migration adults, post-migration adults, and post-regulation adults in deme α are c α , c * α , c * * α , and c α : Because no individuals are lost during migration, the following must hold: The (joint) probability that an adult is in deme α and migrates to deme β can be expressed in terms of the forward and backward migration rates as follows: Inserting (3.6a) into (3.7), we obtain the desired connection between the forward and the backward migration rates: Therefore, ifM is given, an ansatz for the vector c * = (c * 1 , . . . , c * Γ ) in terms of c = (c 1 , . . . , c Γ ) is needed to compute M (as well as a hypothesis for the variation, if any, of c).
Two frequently used assumptions are the following [23]. 1) Soft selection. This assumes that the fraction of adults in every deme is fixed, i.e., c * α = c α for every α ∈ G .
This may be a good approximation if the population is regulated within each deme, e.g., because individuals compete for resources locally [28].
2) Hard selection. Following [28], the fraction of adults will be proportional to mean fitness in the deme if the total population size is regulated. This has been called hard selection and is defined by is the mean fitness of the total population. Essentially, these two assumptions are at the extremes of a broad spectrum of possibilities. Soft selection will apply to plants; for animals many schemes are possible.

REINHARD BÜRGER
Under soft selection, (3.8) becomes As a consequence, if c is constant (c = c), M is constant if and only ifM is constant. If there is no population regulation after migration, then c will generally depend on time because (3.6a) yields c = c * * =M c. Therefore, the assumption of constant deme proportions, c = c, will usually require that population control occurs after migration.
A migration pattern that does not change deme proportions (c * * α = c * α ) is called conservative. Under this assumption, (3.7) yields and, by stochasticity of M andM , we obtain If there is soft selection and the deme sizes are equal (c * α = c α ≡ constant), then m αβ =m βα .
Remark 3.1. Conservative migration has two interesting special cases.
1) Dispersal is called reciprocal if the number of individuals that migrate from deme α to deme β equals the number that migrate from β to α: If this holds for all pairs of demes, then (3.6a) and (3.4) immediately yield c * * β = c * β . From (3.7), we infer m αβ =m αβ , i.e., the forward and backward migration matrices are identical.
2) A migration scheme is called doubly stochastic if αm αβ = 1 for every α . (3.16) If demes are of equal size, then (3.6a) shows that c * * α = c * α . Hence, with equal deme sizes a doubly stochastic migration pattern is conservative. Under soft selection, deme sizes remain constant without further population regulation. Hence, m αβ = m βα and M is also doubly stochastic.
Doubly stochastic migration patterns arise naturally if there is a periodicity, e.g., because the demes are arranged in a circular way. If we posit equal deme sizes and homogeneous migration, i.e.,m αβ =m β−α so that migration rates depend only on distance, then the backward migration pattern is also homogeneous because m αβ =m βα =m α−β and, hence, depends only on β − α. If migration is symmetric, m αβ =m βα , and the deme sizes are equal, then dispersion is both reciprocal and doubly stochastic.
3.3. Important special migration patterns. We introduce three migration patterns that play an important role in the population genetics and ecological literature.  [26]. This model assumes that a proportion µ ∈ [0, 1] of individuals in each deme leaves their deme and is dispersed randomly across all demes. Thus, they perform outbreeding whereas a proportion 1−µ remains at their home site. If c * * α is the proportion (of the total population) of post-migration adults in deme α, then the forward migration rates are defined bym αβ = µc * * β if α = β, andm αα = 1 − µ + µc * * α . If µ = 0, migration is absent; if µ = 1, the Levene model is obtained (see below). Because this migration pattern is reciprocal, M =M holds.
Example 3.4. In the linear stepping-stone model the demes are arranged in a linear order and individuals can reach only one of the neighboring demes. It is an extreme case among migration patterns exhibiting isolation by distance, i.e., patterns in which migration diminishes with the distance from the parental deme. In the classical homogeneous version, the forward migration matrix is We leave it to the reader to derive the backward migration matrix using (3.8). It is a special case of the following general tridiagonal form: where n α ≥ 0 and q α + n α + r α = 1 for every α, q α > 0 for α ≥ 2, r α > 0 for α ≤ Γ − 1, and q 1 = r Γ = 0. This matrix admits variable migration rates between neighboring demes. If all deme sizes are equal, the homogeneous matrix (3.21) satisfies M =M , and each deme exchanges a fraction m of the population with each of its neighboring demes. The stepping-stone model has been used as a starting point to derive the partial differential equations for selection and dispersal in continuous space [72]. Also circular and infinite variants have been investigated.
Juvenile migration is of importance for many marine organisms and plants, where seeds disperse. It can be treated in a similar way as adult migration. Also models with both juvenile and adult migration have been studied. Some authors investigated migration and selection in dioecious populations, as well as selection on X-linked loci (e.g. [73], pp. 143, 144).
Unless stated otherwise, throughout this survey we assume that the backward migration matrix M is constant, as is the case for soft selection if deme proportions and the forward migration matrix are constant. Then the recurrence equations (3.5) provide a self-contained description of the migration-selection dynamics. Hence, they are sufficient to study evolution for an arbitrary number of generations. 4. Two alleles. Of central interest is the identification of conditions that guarantee the maintenance of genetic diversity. Often it is impossible to determine the equilibrium structure in detail because establishing existence and, even more so, stability or location of polymorphic equilibria is unfeasible. Below we introduce an important concept that is particularly useful to establish maintenance of genetic variation at diallelic loci. Throughout this section we consider a single locus with two alleles. The number of demes, Γ, can be arbitrary.

Protected polymorphism.
There is a protected polymorphism [86] if, independently of the initial conditions are, a polymorphic population cannot become monomorphic. Essentially, this requires that if an allele becomes very rare, its frequency must increase. In general, a protected polymorphism is neither necessary nor sufficient for the existence of a stable polymorphic equilibrium. For instance, on the one hand, if there is an internal limit cycle that attracts all solutions, then there is a protected polymorphism. On the other hand, if there are two internal equilibria, one asymptotically stable, the other unstable, then selection may remove one of the alleles if sufficiently rare. A generalization of this concept to multiple alleles would correspond to the concept of permanence often used in ecological models (e.g. [51]).
Because we consider only two alleles, we can simplify the notation. We write p α = p 1,α for the frequency of allele A 1 in deme α (and 1 − p α for that of A 2 in deme α). Let p = (p 1 , . . . , p Γ ) denote the vector of allele frequencies. Instead of using the fitness assignments w 11,α , w 12,α , and w 22,α , it will be convenient to scale the fitness of the three genotypes in deme α as follows . This can be achieved by setting x α = w 11,α /w 12,α and y α = w 22,α /w 12,α , provided w 12,α > 0. With these fitness assignments, one obtains and the migration-selection dynamics (3.5) becomes We consider this as a (discrete) dynamical system on [0, 1] Γ . We call allele A 1 protected if it cannot be lost. Thus, it has to increase in frequency if rare. In mathematical terms this means that the monomorphic equilibrium p = 0 must be unstable. To derive a sufficient condition for instability of p = 0, we linearize (4.3) at p = 0. If y α > 0 for every α (which means that A 2 A 2 is nowhere lethal), a simple calculation shows that the Jacobian of (4.3a), is a diagonal matrix with (nonzero) entries d αα = y −1 α . Because (4.3b) is linear, the linearization of (4.3) is i.e., q αβ = m αβ /y β . To obtain a simple criterion for protection, we assume that the descendants of individuals in every deme be able eventually to reach every other deme. Mathematically, the appropriate assumption is that M is irreducible. Then Q is also irreducible and it is nonnegative. Therefore, the Theorem of Perron and Frobenius (e.g. [90]) implies the existence of a uniquely determined eigenvalue λ 0 > 0 of Q such that |λ| ≤ λ 0 holds for all eigenvalues of Q. In addition, there exists a strictly positive eigenvector pertaining to λ 0 which, up to multiplicity, is uniquely determined. As a consequence, (if λ 0 = 1, then stability cannot be decided upon linearization). This maximal eigenvalue satisfies with equality if and only if all the row sums are the same.
Example 4.1. Suppose that A 2 A 2 is at least as fit as A 1 A 2 in every deme and more fit in at least one deme, i.e., y α ≥ 1 for every α and y β > 1 for some β. Then q αβ = m αβ /y β ≤ m αβ for every β. Because M is irreducible, there is no β such that m αβ = 0 for every α. Therefore, the row sums β q αβ = β m αβ /y β in (4. 7) are not all equal to one, and we obtain Thus, A 1 is not protected, and this holds independently of the choice of the x α , or w 11,α . It can be shown similarly that A 1 is protected if A 1 A 2 is favored over A 2 A 2 in at least one deme and is nowhere less fit than A 2 A 2 .
One obtains the condition for protection of A 2 if, in (4.6), A 1 is replaced by A 2 and λ 0 is the maximal eigenvalue of the matrix with entries m αβ /x β . Clearly, there is a protected polymorphism if both alleles are protected.
In the case of complete dominance the eigenvalue condition (4.6) cannot be satisfied. Consider, for instance, protection of A 1 if A 2 is dominant, i.e., y α = 1 for every α. Then q αβ = m αβ , β q αβ = β m αβ = 1, and λ 0 = 1. This case is treated in Section 6.2 of [73].

4.2.
Two demes. It will be convenient to set where we assume r α < 1 and s α < 1 for every α ∈ {1, 2}. We write the backward migration matrix as where 0 < m α < 1 for every α ∈ {1, 2}. Now we derive the condition for protection of A 1 . The characteristic polynomial of Q is proportional to It is convex and satisfies (4.13) By Example 4.1, A 1 is not protected if A 1 A 2 is less fit than A 2 A 2 in both demes (more generally, if s 1 ≤ 0, s 2 ≤ 0, and s 1 + s 2 < 0). Of course, A 1 will be protected if A 1 A 2 is fitter than A 2 A 2 in both demes (more generally, if s 1 ≥ 0, s 2 ≥ 0, and s 1 + s 2 > 0). Hence, we restrict attention to the most interesting case when A 1 A 2 is fitter than A 2 A 2 in one deme and less fit in the other, i.e., s 1 s 2 < 0.
If there is no dominance (r α = −s α and 0 < |s α | < 1 for α = 1, 2), then further simplification can be achieved. From the preceding paragraph the results depicted in Figure 4.2 are obtained. The region of a protected polymorphism is In a panmictic population, a stable polymorphism can not occur in the absence of overdominance. Protection of both alleles in a subdivided population requires that selection in the two demes is in opposite direction and sufficiently strong relative to migration. Therefore, the study of the maintenance of polymorphism is of most interest if selection acts in opposite direction and dominance is intermediate, i.e., r α s α < 0 for α = 1, 2 and s 1 s 2 < 0. In the first, the fitnesses are given by thus, there is deme independent degree of dominance [75]. In the second scenario, the fitnesses are given by In both cases, we assume 0 < s < 1, 0 < a < 1 (a is a measure of the selection intensity in deme 2 relative to deme 1), −1 ≤ h ≤ 1 (intermediate dominance), and m 1 = m 2 = m.
If, in (4.17), selection is sufficiently symmetric, i.e., a > 1/(1 + 2s), there exists a protected polymorphism for every h ≤ 1 and every m ≤ 1. Otherwise, for given a, the critical migration rate m admitting a protected polymorphism decreases with increasing h because this increases the average (invasion) fitness of A 1 . This can be shown analytically by studying the principal eigenvalue.
For (4.18), increasing h greatly facilitates protected polymorphism because it leads to an increase of the average fitness of heterozygotes relative to the homozygotes. The precise argument is as follows. If we rescale fitnesses in (4.18) according to (4.1), the matrix Q in (4.5) has the entries q α1 = m α1 (1 + hs)/(1 − s) and q α2 = m α2 (1 + has)/(1 − as), which are increasing in h. Therefore, the principal eigenvalue λ 0 of Q increases in h (e.g. [11], Chapter 1.3), and (4.6) implies that protection of A 1 is facilitated by increasing h. Because an analogous reasoning applies to A 2 , the result is proved.
Indeed, the above finding on the role of dominance for the fitness scheme (4.18) is a special case of the following result of Nagylaki (personal communication). Assume fitnesses of A 1 A 1 , A 1 A 2 , and A 2 A 2 in deme α are 1 + s α , 1 + h α s α , and 1, respectively, where s α > −1 and 0 ≤ h α ≤ 1. If the homozygote fitnesses are fixed, increasing the heterozygote fitness in each deme aids the existence of a protected polymorphism. The proof follows by essentially the same argument as above.
where s 1 s 2 < 0. Therefore, for given s 1 , s 2 , and c 1 , there is a critical value µ 0 such that allele A 1 is protected if and only if µ < µ 0 . This implies that for two diallelic demes a protected polymorphism is favored by a smaller migration rate.
Example 4.4. In the Levene model, the condition for a protected polymorphism is We close this subsection with an example showing that already with two alleles and two demes the equilibrium structure can be quite complicated.
Example 4.5. In the absence of migration, the recurrence equations for the allele frequencies p 1 , p 2 in the two demes are two decoupled one-locus selection dynamics of the form (2.16). Therefore, if there is underdominance in each deme, the top convergence pattern in Figure 2.1 applies to each deme. As a consequence, in the absence of migration, the complete two-deme system has nine equilibria, four of which are asymptotically stable and the others are unstable. Under sufficiently weak migration all nine equilibria are admissible and the four stable ones remain stable, whereas the other five are unstable. Two of the stable equilibria are internal. For increasing migration rate, several of these equilibria are extinguished in a sequence of bifurcations [56].

Arbitrary number of demes.
The following result is a useful tool to study protection of an allele. Let I (n) and Q (n) respectively designate the n × n unit matrix and the square matrix formed from the first n rows and columns of Q. for some n, where 1 ≤ n ≤ Γ, then A 1 is protected.
This theorem is sharp in the sense that if the inequality in (4.21) is reversed for every n ≤ Γ, then A 1 is not protected.
The simplest condition for protection is obtained from Theorem 4.6 by setting n = 1. Hence, A 1 is protected if an α exists such that [27] m αα /y α > 1 . This condition ensures that, when rare, the allele A 1 increases in frequency in deme α even if this subpopulation is the only one containing the allele. The general condition in the above theorem ensures that A 1 increases in the n subpopulations if they are the only ones that contain it. Therefore, we get the following simple sufficient condition for a protected polymorphism: There exist α and β such that m αα /y α > 1 and m ββ /x β > 1. If we apply these results to the Deakin model with an arbitrary number of demes, condition (4.23) becomes There exist α and β such that A more elaborate and less stringent condition follows from Theorem 4.6 by nice matrix algebra: 22]). For the Deakin model with Γ ≥ 2 demes, the following is a sufficient condition for protection of A 1 . There exists a deme α such that If (4.25) is violated for every α and the inequality in (4.26) is reversed, then A 1 is not protected.
Corollary 4.7 can be extended to the inhomogeneous Deakin model, which allows for different homing rates ( [22]; [53], pp. 85, 127). Using Corollary 4.7, we can generalize the finding from two diallelic demes that a lower degree of outbreeding is favorable for protection of one or both alleles. More precisely, we show: Proof. If condition (4.25) holds for µ 2 , it clearly holds for µ 1 . Now suppose that 1 − y α < µ 1 for every α (hence 1 − y α < µ 2 ) and that µ 2 satisfies (4.26). This is equivalent to which proves our assertion.
Karlin ([53], Theorem 5.2 and p. 128) extended this result to migration matrices of the form where M is irreducible, by proving that for any diagonal matrix D with positive entries on the diagonal the spectral radius of the matrix M (µ) D is decreasing as µ increases (strictly if D = dI). From the criterion (4.6) of protection, one obtains Corollary 4.9. Assume the family of backward migration matrices (4.29). If 0 < µ 1 < µ 2 ≤ 1, then allele A 1 is protected for µ 1 if it is protected for µ 2 . Therefore, increased mixing reduces the potential for a protected polymorphism.
As noted by Altenberg [5], Karlin's theorem on the monotonicity of the spectral radius is related to a result of Kato on the convexity of the spectral bound of so-called resolvent-positive operators. Altenberg provides a generalization of both results implying that increased mixing leads to a certain "reduction" phenomenon in a large class of models. In general, however, it is not true that less migration favors the maintenance of polymorphism. For instance, if the amount of homing, 1 − µ α , varies among demes, a protected polymorphism may be destroyed by decreasing one µ α ( [53], p. 128).
Example 4.10. We apply Corollary 4.7 to the Levene model. Because (4.25) can never be satisfied if µ = 1, A 1 is protected from loss if the harmonic mean of the y α is less than one, i.e., if Jointly, (4.30a) and (4.30b) provide a sufficient condition for a protected polymorphism. If y * > 1 or x * > 1, then [86]. Therefore, a sufficient condition for a protected polymorphism is Whereas in the Deakin model and in its special case, the Levene model, dispersal does not depend on the geographic distribution of niches, in the stepping-stone model it occurs between neighboring niches. How does this affect the maintenance of a protected polymorphism? Here is the answer: Example 4.11. For the general stepping-stone model (3.22) with fitnesses given by (4.1), Karlin [53] provided the following explicit criterion for protection of A 1 : where π 1 = 1 and π α = r α−1 r α−2 · · · r 1 q α q α−1 · · · q 2 (4.34) if α ≥ 2. This result is a simple consequence of Lemma 4.14.
For the homogeneous stepping-stone model with equal demes sizes, i.e., M given by (3.21), (4.33) simplifies to which is the same as condition (4.30) in the Levene model with equal deme sizes.
Thus, we obtain the surprising result that the conditions for protection are the same in the Levene model and in the homogeneous stepping stone model provided all demes have equal size.

4.4.
The continent-island model. We consider a population living on an island. Initially, we admit an arbitrary number of alleles. Each generation, a proportion a of adults is removed by mortality or emigration and a fraction b of migrants with constant allele frequenciesq i is added. A simple interpretation is that of one-way migration from a continent to an island. The continental population is in equilibrium with allele frequenciesq i > 0. Sometimes,q i is interpreted as the average frequency over (infinitely) many islands and the model is simply called island model.
Assuming random mating on the island, the dynamics of allele frequencies p i on the island becomes where w i is the fitness of allele A i on the island and m = b/(1 − a + b). Thus, m is the fraction of zygotes with immigrant parents. Throughout we assume 0 < m < 1.
If we define u ij = mq j and consider u ij as the mutation rate from A i to A j , a special case of the (multiallelic) mutation-selection model is obtained (the socalled house-of-cards model). Therefore ( [13], pp. 102-103), (4.36) has the Lyapunov function (4.37) It follows that all trajectories are attracted by the set of equilibria.
In the sequel we determine the conditions under which an 'island allele' persists in the population despite immigration of other alleles from the continent. Clearly, no allele carried to the island recurrently by immigrants can be lost.
We investigate the diallelic case and consider alleles A 1 and A 2 with frequencies p and 1 − p on the island, andq 2 = 1 on the continent. Thus, all immigrants are A 2 A 2 . For the fitnesses on the island, we assume where 0 < s ≤ 1 and −1 ≤ h ≤ 1. Therefore, A 1 evolves according to We outline the analysis of (4.39). Since f (1) = 1 − m < 1, this confirms that A 2 cannot be lost (p = 1 is not an equilibrium). Because Obviously, p = 0 is an equilibrium. If there is no other equilibrium, then it must be globally asymptotically stable (as is also obvious from f (1) < 1, which implies p < p).
The equilibria with p = 0 satisfȳ which is quadratic in p. With the fitnesses (4.38), the solutions arê These solutions give rise to feasible equilibria if 0 ≤p ± ≤ 1. As h → 0, (4.43) gives the correct limit,p We define Then µ 2 > µ 1 if and only if h < h 0 . If µ 2 < µ 1 , thenp + is not admissible. As m increases from 0 to 1, h 0 decreases from − 1 3 to −1. Hence, if there is no dominance or the fitter allele For fixed h and m it is straightforward, but tedious, to study the dependence ofp + andp − on s. The results of this analysis can be summarized as follows (for an illustration, see (b) Let m/s > µ 2 and −1 ≤ h < h 0 , or m/s ≥ µ 1 and h 0 ≤ h ≤ 1, i.e., migration is strong relative to selection. Then there exists no internal equilibrium and p(t) → 0 as t → 0, i.e., the continental allele A 2 becomes fixed.
(c) Let µ 1 < m/s ≤ µ 2 and −1 ≤ h < h 0 , i.e., migration is moderately strong relative to selection and the fitter allele is (almost) recessive. Then there are the three equilibria 0,p + , andp − , where 0 <p + <p − , and For a detailed proof, see [73], Chapter 6.1. The global dynamics of the diallelic continent-island model admitting evolution on the continent is derived in [75].
Sometimes immigration from two continents is considered [24]. Some authors call our general migration model with n demes the n-island model (e.g. [24]). The following (symmetric) island model is a standard model in investigations of the consequences of random drift and mutation in finite populations: Clearly, migration is doubly stochastic in this model and there is no isolation by distance. In the special case when all deme sizes are equal, i.e., c α = 1/Γ, (4.47) reduces to a special case of the Deakin model, (3.19), with m = µ(1 − Γ −1 ). In this, and only this case, the island model is conservative and the stationary distribution of M is c = (1/Γ, . . . , 1/Γ).

Submultiplicative fitnesses.
Here, we admit arbitrary migration patterns but assume that fitnesses are submultiplicative, i.e., the fitnesses in (4.1) satisfy for every α [54]. Submultiplicative fitnesses include a number of important cases: multiplicative fitnesses, no dominance, partial or complete dominance of the fitter allele, and overdominance. We recall that with multiplicative fitnesses (x α y α = 1), the diploid selection dynamics (2.6) reduces to the haploid dynamics (2.9).

Theorem 4.13 ([54], Result I).
If Γ ≥ 2 and fitnesses are submultiplicative in each deme, at most one monomorphic equilibrium can be asymptotically stable and have a geometric rate of convergence.
The proof of this theorem applies the following very useful result to the conditions for protection of A 1 or A 2 (Section 4.1). Throughout, we denote the left principal eigenvector of M by ξ ∈ S Γ ; cf. (7.23). Without restrictions on the migration matrix, the sufficient condition for protection of A 1 is sharp. This can be verified for a circulant stepping stone model, when ξ α = 1/Γ. The conclusion of Theorem 4.13 remains valid under some other conditions (see [54]). For instance, if migration is doubly stochastic, then ξ α = 1/Γ for every α and implies the conclusion of the theorem.
Here is another interesting result: Because the proof given in [54] is erroneous, we present a corrected proof (utilizing their main idea).
Proof of Theorem 4.15. Because fitnesses are multiplicative, we can treat the haploid model. Without loss of generality we assume that the fitnesses of alleles A 1 and A 2 in deme 1 are s 1 ≥ 1 and 1, respectively, and in deme 2 they are s 2 (0 < s 2 < 1) and 1. (We exclude the case when one allele is favored in both demes, hence goes to fixation.) Let the frequency of A 1 in demes 1 and 2 be p 1 and p 2 , respectively. Then the recursion (4.3) simplifies to Solving p 1 = p 1 for p 2 , we obtain . (4.53) Substituting this into p 2 = p 2 , we find that every equilibrium must satisfy We want to show that A(p 1 ) has always two zeros, because then zeros in (0, 1) can emerge only by bifurcations through either p 1 = 0 or p 1 = 1. Therefore, we calculate the discriminant where we set s 1 = 1 + σ and s 2 = 1 − τ with σ ≥ 0 and 0 < τ < 1. Now we consider D as a (quadratic) function in and We distinguish two cases. 1. If m 2 (2 + σ − τ ) ≥ στ , then D (0) ≥ 0 and there can be no zero in (0, 1).
Therefore, D has no zero in (0, 1). Thus, we have shown that A(p 1 ) = 0 has always two real solutions. For the rest of the proof we can follow Karlin and Campbell: If s 1 = 1 and s 2 < 1, there are no polymorphic equilibria, p 1 = 0 is stable and p 1 = 1 is unstable. As s 1 increases from 1, a bifurcation event at p 1 = 0 occurs before p 1 = 1 becomes stable. At this bifurcation event, p 1 = 0 becomes unstable and a stable polymorphic equilibrium bifurcates off is linear in s 1 , no further polymorphic equilibria can appear as s 1 increases until p = 1 becomes stable, which then remains stable for larger s 1 .
Karlin and Campbell [54] pointed out that the dynamics (4.52) has the following monotonicity property: If p 1 < q 1 and p 2 < q 2 , then p 1 < q 1 and p 2 < q 2 . By considering 1 − p 2 and 1 − q 2 , p 1 < q 1 and p 2 > q 2 implies p 1 < q 1 and p 2 > q 2 . In particular, if p 1 < p 1 and p 2 < p 2 , then p 1 < p 1 and p 2 < p 2 . Therefore, this monotonicity property implies global monotone convergence to the asymptotically stable equilibrium.
Remark 4.16. (a) Because the monotonicity property used in the above proof holds for an arbitrary number of demes and for selection on diploids, Karlin and Campbell [54] stated that "periodic or oscillating trajectories seem not to occur" (in diallelic models). However, the argument in the above proof cannot be extended to more than two demes, and Akin (personal communication) has shown that for three diallelic demes unstable periodic orbits can occur in the corresponding continuoustime model (cf. Section 6.2). For two diallelic demes, however, every trajectory converges to an equilibrium.
(b) It would be of interest to prove that attracting periodic orbits exist. In view of the results derived or discussed above, the simplest possible scenarios seem to be three alleles and two demes or two alleles and four demes.
In [54] it was also shown that for some migration patterns, Theorem 4.15 remains true for an arbitrary number of demes. The Levene model is one such example. Theorem 5.16 provides a stronger result. It is an open problem whether Theorem 4.15 holds for arbitrary migration patterns if there are more than two demes. A further interesting result is the following:   5. The Levene model. This is a particularly simple model to examine the consequences of spatially varying selection for the maintenance of genetic variability. It can also be interpreted as a model of frequency-dependent selection (e.g. [101]). From the definition (3.20) of the Levene model and equation (3.5a), we infer which is independent of α. Therefore, after one round of migration allele frequencies are the same in all demes, whence it is sufficient to study the dynamics of the vector of allele frequencies The migration-selection dynamics (3.5) simplifies drastically and yields the recurrence equation for the evolution of allele frequencies in the Levene model. Therefore, there is no population structure in the Levene model although the population experiences spatially varying selection. In particular, distance does not play any role. The dynamics (5.3) is also obtained if, instead of the life cycle in Section 3.2, it is assumed that adults form a common mating pool and zygotes are distributed randomly to the demes.
Remark 5.1. 1. In the limit of weak selection, the Levene model becomes equivalent to panmixia. Indeed, with w ij,α = 1 + sr ij,α , an argument analogous to that below (2.20) shows that the resulting dynamics is of the forṁ [78]. Therefore, the Levene model is of genuine interest only if selection is strong. For arbitrary migration, the weak-selection limit generally does not lead to a panmictic dynamics (Section 6.3).
2. The weak-selection limit for juvenile migration is panmixia [82]. 3. For hard selection (3.10), the dynamics in the Levene model becomes which is again equivalent to the panmictic dynamics with fitnesses averaged over all demes. Therefore, no conceptually new behavior occurs. With two alleles, there exists a protected polymorphism if and only if the averaged fitnesses display heterozygous advantage. Interestingly, it is not true that the condition for protection of an allele under hard selection is always more stringent than under soft selection, although under a range of assumptions this is the case ( [73], Sect. 6.3).

5.1.
General results about equilibria and dynamics. We definẽ which is the geometric mean of the mean fitnesses in single demes. Furthermore, we define Both functions will play a crucial role in our study of the Levene model. From we obtain Because i p i = 1, we can write (5.9) in the form A simple exercise using Lagrange multipliers shows that the equilibria of (5.10), hence of (5.3), are exactly the critical points ofw, or F . Our aim is to prove the following theorem ( [64], [20], [73]): (a) Geometric-mean fitness satisfies ∆w(p) ≥ 0 for every p ∈ S I , and ∆w(p) = 0 if and only if p is an equilibrium of (5.3). The same conclusion holds for F .
(b) The set Λ of equilibria is globally attracting, i.e., p(t) → Λ as t → ∞. If every point in Λ is isolated, as is generic 6 , then p(t) converges to somep ∈ Λ. Generically, p(t) converges to a local maximum ofw.
This is an important theoretical result because it states that in the Levene model no complex dynamical behavior, such as limit cycles or chaos, can occur. One immediate consequence of this theorem is Outline of the proof of Theorem 5.2. (a) We assume c α ∈ Q for every α. Then there exists n ∈ N such that W =w n is a homogeneous polynomial in p with nonnegative coefficients. From From the inequality of Baum and Eagon [10], we infer immediately W (p ) > W (p) unless p = p. Therefore, W andw are (strict) Lyapunov functions. By a straightforward but tedious approximation argument, which uses compactness of S I , it follows thatw is also a Lyapunov function if c α ∈ R (see [73], Sect. 6.3). This proves the first assertion in (a). The second follows because the logarithm is strictly monotone increasing.
(b) The first statement is an immediate consequence of LaSalle's invariance principle ( [60], p. 10). The second is obvious, and the third follows becausew is nondecreasing. (Note, however, that convergence to a maximum holds only generically because some trajectories may converge to saddle points. Moreover, becausew may also have minima and saddle points, the globally attracting set Λ may not be stable; only a subset of Λ is stable.) Next, we derive a useful criterion for the existence of a unique stable equilibrium.
Lemma 5.4. Ifw α is concave for every α, then F is concave.
Proof. Because the logarithm is strictly monotone increasing and concave, a simple estimate shows that lnw α is concave. Hence, F is concave. 78]). If F is concave, there exists exactly one stable equilibrium (point or manifold) and it is globally attracting. If there exists an internal equilibrium, it is the global attractor.
Proof. Concavity and the fact that ∆F ≥ 0 imply that if an internal equilibrium exists, it must be the global attractor.
Suppose now there exist two stable equilibrium points on the boundary of the simplex S I . Because F is concave, F is constant on the line that joins them. If that line is on the boundary of S I , then these two equilibrium points are elements of the same manifold of equilibria (on the boundary of S I ). If the line connecting them is an internal equilibrium, then the case treated above applies.
If dominance is absent in every deme, there exist constants v i,α such that for every i, j ∈ I and α ∈ G, whence (3.2) yields A simple calculation shows that without dominance the dynamics (5.3) becomes Fitnesses are called multiplicative if there exist constants v i,α such that for every i, j ∈ I and α ∈ G. Then (3.2) yields With multiplicative fitnesses, one obtains the haploid Levene model, i.e., (c) In every deme a globally attracting internal equilibrium exists without migration (i.e., with the dynamics p i,α = p i,α w i,α /w α ).
Therefore, the conclusions of Theorem 5.5 apply in each of the cases.
(c) If there is a globally attracting internal equilibrium in deme α when it is isolated, then the fact that ∆w α ≥ 0 implies that the quadraticw α must be concave.
In all three cases, the assertion follows from Lemma 5.4.
Remark 5.7. Theorems 5.5 and 5.6 hold for hard selection if we replace F bȳ w. For haploids, Strobeck [95] proved uniqueness and asymptotic stability of an internal equilibrium.
In [80] sufficient conditions are provided for the nonexistence of an internal equilibrium.
The following theorem shows that a globally stable internal equilibrium can be maintained on an open set of parameters, even in the absence of overdominance. Clearly, this requires that there are at least two demes; cf. Section 2.3. We say there is partial dominance in every deme if holds for every pair i, j ∈ I with i = j and for every α ∈ Γ.  [80]. An interesting question to ask is: What determines the number of alleles that can be maintained at an equilibrium? If there is no dominance, there is a simple answer. Proof. From (5.14) we conclude that, at an internal equilibriump, holds. The substitution The statement about neutral demes follows because in a neutral deme v i,α = v α = 1 for every i, hence x α = c α holds.  Remark 5.11. Theorem 5.9 holds for hard selection. A slight modification of the proof shows that it holds also for multiplicative fitnesses. In fact, in [95] an analog of Theorem 5.9 is proved for haploid species. Let v i,α = u i δ iα , where u i > 0 and δ iα is the Kronecker delta. This means that allele A i has fitness u i in deme i and fitness 0 elsewhere. Hence, every allele is the best in one deme. Then (5.14) simplifies to In the formulation of Theorem 5.9, the word 'generic' is essential. If Γ and I are arbitrary and one assumes w ij,α = 1 + r ij g α (5.23) for (sufficiently small) constants r ij and g α , it can be shown that an internal (hence globally attracting) manifold of equilibria exists for an open set of parameter combinations. This holds also for the additive case, when r ij = s i + s j [78]. Another interesting problem is to determine conditions when a specific allele will go to fixation. A simple and intuitive result is the following: for every j = i. Then p i (t) → 1 as t → ∞.
A proof as well as additional results on the loss or fixation of alleles may be found in [80]. The following example is based on a nice application of Theorem 5.6 and shows that migration may eliminate genetic variability.
Example 5.14 ( [78]). We suppose two demes of equal size (c 1 = c 2 = 1 2 ) and three alleles without dominance. The alleles A 1 and A 2 have extreme fitnesses and A 3 is intermediate in both demes. More precisely, we assume v 1,1 = 1, v 2,1 = 0, v 3,1 = u, and v 1,2 = 0, v 2,2 = 1, v 3,2 = u, where 0 < u < 1. Without migration, A 1 is ultimately fixed in deme 1 and A 2 is fixed in deme 2. We now establish that the Levene model evolves as sketched in Figure 5.1, i.e., Because there is no dominance, Theorem 5.6 informs us that there can be only one stable equilibrium. Therefore, it is sufficient to establish asymptotic stability of the limiting equilibrium, which we leave as an easy exercise.
Following [75], we say there is deme-independent degree of intermediate domi- holds for constants ϑ ij such that for every α and every pair i, j. In particular, ϑ ii = 1 2 . Obviously, DIDID covers complete dominance or recessiveness (ϑ ij = 0 or ϑ ij = 1 if i = j), and no dominance (ϑ ij = 1 2 ), but not multiplicativity. We also note that DIDID includes the biologically important case of absence of genotype-byenvironment interaction. In general, F is not concave under DIDID. However, Nagylaki ([75], Theorem 3.2) proved that with DIDID the evolution of p(t) is qualitatively the same as without dominance. Therefore, the conclusion of Theorem 5.5 holds, i.e., under DIDID there exists exactly one stable equilibrium (point or manifold) and it is globally attracting. If there exists an internal equilibrium, it is the global attractor. Moreover, the number of demes is a generic upper bound on the number of alleles that can segregate at an equilibrium (generalization of Theorem 5.9). This contrasts sharply with Theorem 5.8. Finally, the condition for loss of an allele in Theorem 5.13 has a simple generalization to DIDID. For proofs and further results about DIDID consult [75].
Remark 5.15. For general migration, the dynamics under DIDID may differ qualitatively from that under no dominance. For instance, in a diallelic model with one way migration, two asymptotically stable equilibria may coexist if there is DIDID [75]. Moreover, in two triallelic demes, a globally asymptotically stable internal equilibrium exists for an open set of parameters [84]. Therefore, with arbitrary migration and DIDID, the number of alleles maintained at a stable equilibrium can be higher than the number of demes. This is not the case in the absence of dominance when the number of demes is a generic upper bound (Theorem 6.1).

5.3.
Two alleles with dominance. Here, we specialize to two alleles but admit arbitrary dominance. We assume that fitnesses are given by (4.1). Our first result yields a large class of fitness schemes for which a unique stable equilibrium exists.
Then F is concave on [0, 1]. Hence, there exists at most one internal equilibrium. If an internal equilibrium exists, it is globally asymptotically stable. If a monomorphic equilibrium is stable, then it is globally asymptotically stable.
In the proof it is shown that lnw α is concave if and only if (5.27) is fulfilled. Then the conclusion follows from Theorem 5.5.
Theorem 5.16 shows that the protection conditions (4.30) imply the existence of a globally asymptotically stable internal equilibrium if (5.27) holds. It generalizes Result I in [52], where uniqueness of an internal equilibrium and global convergence were proved under the more restrictive assumption of submultiplicative fitnesses (4.48). The proof there is quite different.
where κ is as above. The unique, asymptotically stable equilibrium is given bŷ p = 1 − κ.   30) and there exists at least one internal equilibrium. Then both monomorphic equilibria are asymptotically stable and the internal equilibrium is unique.  Table 5.1 presents numerical results that demonstrate that stronger selection, submultiplicativity, and a higher number of demes all facilitate protected polymorphism in the Levene model. These findings agree with intuition. For instance, if an allele has sufficiently low fitness in just one deme, i.e., c α /y α > 1, the other allele is protected.
So far, we derived sufficient conditions for a unique (internal) equilibrium, but we have not yet considered the question of how many (stable) internal equilibria can coexist. This turns out to be a difficult question which was solved only recently for diallelic loci.
With fitnesses given by (4.1), a simple calculation shows that the dynamics (5.3) can be written as (5.31) Expressing this with a common denominator, we see that the internal equilibria are the solutions of a polynomial in p of degree 2Γ − 1. Thus, in principle, there can be up to 2Γ − 1 internal equilibria. For two demes, Karlin [52] provided a construction principle for obtaining three internal equilibria. It requires quite extreme fitness differences and that in both demes the less fit allele is close to dominant. By a clever procedure, Novak [83] proved the following result: Theorem 5.20. The diallelic Levene model allows for any number j ∈ {1, . . . , 2Γ− 1} of hyperbolic internal equilibria with any feasible stability configuration.
His numerical results, which admit arbitrary dominance, show that the proportion of parameter space supporting more than Γ equilibria becomes extremely small if Γ > 2.
6. Multiple alleles and arbitrary migration. Because for multiple alleles and arbitrary migration few general results are available, we focus on three limiting cases that are biologically important and amenable to mathematical analysis. These are weak selection and weak migration, weak migration (relative to selection), and weak selection (relative to migration). The first leads to the continuous-time migration model, the second to the so-called weak-migration limit, and the third to the strongmigration limit. The latter two are based on a separation of time scales. With the help of perturbation theory, results about existence and stability of equilibria, but also about global convergence, can be derived for an open subset of parameters of the full model. In addition, we report results about the case of no dominance and about uniform selection, i.e., when selection is the same in every deme. We start with no dominance.  This theorem also holds for hard selection. As shown by Peischl [84], it can not be extended to DIDID; cf. Remark 5.15. An example showing that the upper bound can be achieved if I = Γ is obtained by setting v i,α = u i δ iα , where u i > 0. The internal equilibrium can be written asp = 1 2 (I − 1 2 M ) −1 M , and convergence is geometric (at a rate ≤ 1 2 ). As we shall see below, with dominance the number of alleles that can be maintained at equilibrium may depend on the strength and pattern of migration.
Example 6.2. Eyland [32] provided a global analysis of (6.9) for the special case of two diallelic demes without dominance. As in Sect. 4.2, we assume that the fitnesses of A 1 A 1 , A 1 A 2 , and A 2 A 2 in deme α are 1 + s α , 1, and 1 − s α , respectively, where s α = 0 (α = 1, 2). Moreover, we set µ 1 = µ 12 > 0, µ 2 = µ 21 > 0, and write p α for the frequency of A 1 in deme α. Then (6.9) becomeṡ The equilibria can be calculated explicitly. At equilibrium, p 1 = 0 if and only if p 2 = 0, and p 1 = 1 if and only if p 2 = 1. In addition, there may be an internal equilibrium point. We set and 13) The internal equilibrium exists if and only if s 1 s 2 < 0 and |κ| < 1; cf. (4.15). If s 2 < 0 < s 1 , it is given bŷ 14) It is straightforward to determine the local stability properties of the three possible equilibria. Gobal asymptotic stability follows from the results cited above about quasimonotone systems. Let p = (p 1 , p 2 ) . Then allele A 1 is eliminated in the region Ω 0 in Figure 4.2, i.e., p(t) → (0, 0) as t → ∞, whereas A 1 is ultimately fixed in the region Ω 1 . In Ω + , p(t) converges globally to the internal equilibrium pointp given by (6.14).
For the discrete-time dynamics (3.5) such a detailed analysis is not available. However, for important special cases, results about existence, uniqueness, and stability of equilibria were derived in [54]; see Section 4.5.
Finally, we present sufficient conditions for global loss of an allele. In discrete time, such conditions are available only for the Levene model. For (6.9), however, general conditions were derived in [81]. Suppose that there exists i ∈ I and constants γ ij such that γ ij ≥ 0 , γ ii = 0 , j γ ij = 1 , (6.15a) and j γ ij r jk,α > r ik,α (6.15b) for every α ∈ G and every k ∈ I. . If the matrix (µ αβ ) is irreducible and the conditions (6.15) are satisfied, then p i (t) → 0 as t → ∞ whenever p i (0) > 0 and p j (0) > 0 for every j such that γ ij > 0.
If there is no dominance, constants s i,α exist such that r ij,α = s i,α + s j,α for every i, j, α. Then condition (6.15b) simplifies to j γ ij s j,α > s i,α , (6.16) and Theorem 6.3 applies.
If, in addition to (6.17), we assume that min{s 1,α , s I,α } < s i,α < max{s 1,α , s I,α } (6.18) for i = 2, . . . , I − 1 and every α, then every allele A i with 1 < i < I is intermediate in every deme, and A 1 and A I are extreme. Thus, Theorem 6.3 implies that all intermediate alleles are eliminated. This conclusion can be interpreted as the elimination of generalists by specialists, and it can yield the increasing phenotypic differentiation required for parapatric speciation (cf. [66] for an analogous result and discussion in the context of diffusion models). If s 1,α − s I,α changes sign among demes, so that every allele is the fittest in some deme(s) and the least fit in the other(s), then both alleles may be maintained. Another application of Theorem 6.3 is the following. If in every deme the homozygotes have the same fitness order and there is strict heterozygote intermediacy, then the allele with the greatest homozygous fitness is ultimately fixed (Remark 3.20 in [81]. Remark 6.4. The slow-evolution approximation of the exact juvenile-migration model is also (6.9) ( [82], [73], pp. 143-144).
6.3. Weak migration. If migration is sufficiently weak relative to selection, properties of the dynamics can be inferred by perturbation techniques from the wellunderstood case of a finite number of isolated demes in which there is selection and random mating.
To study weak migration, we assume where µ αβ is fixed for every α, β, and > 0 is sufficiently small. Because M is stochastic, we obtain for every α ∈ Γ (cf. Section 6.2) µ αβ ≥ 0 for every β = α and β µ αβ = 0 . (6.20) If there is no migration ( = 0), the dynamics in each deme reduces to the pure selection dynamics Because (6.21) is defined on the Cartesian product S Γ I , it may exhibit a richer equilibrium and stability structure than the panmictic selection dynamics (2.6). This was already illustrated by Example 4.5, in which two asymptotically stable internal equilibria may coexist. By Theorem 2.2, such an equilibrium configuration does not occur for (2.6).
The following is a central perturbation result that has a number of important consequences. Among others, it excludes complex dynamics in the full system (3.5) provided migration is sufficiently weak. We note that hyperbolicity is a generic property for (6.21) (Appendix A in [77]), and an internal equilibrium is hyperbolic if and only if it is isolated (Lemma 3.2 in [79]). . Suppose that every equilibrium of (6.21) is hyperbolic, that (6.19) holds, and that > 0 is sufficiently small.
(a) The set of equilibria Σ 0 ⊂ S Γ I of (6.21) contains only isolated points, as does the set of equilibria Σ ⊂ S Γ I of (3.5). As → 0, each equilibrium in Σ converges to the corresponding equilibrium in Σ 0 .
(b) In the neighborhood of each asymptotically stable equilibrium in Σ 0 , there exists exactly one equilibrium in Σ , and it is asymptotically stable. In the neighborhood of each unstable internal equilibrium in Σ 0 , there exists exactly one equilibrium in Σ , and it is unstable. In the neighborhood of each unstable boundary equilibrium in Σ 0 , there exists at most one equilibrium in Σ , and if it exists, it is unstable.
(c) Every solution p(t) of (3.5) converges to one of the equilibrium points in Σ .
The perturbation results in (a) and (b) are essentially due to Karlin and McGregor ( [56], [57]). The proof of (c), which also yields (a) and (b), is a simplification of that of Theorem 2.3 in [77] and is based on quite deep results. To outline the proof, we need some preparation.
Consider a family of maps (difference equations) that depends on a parameter , where x ∈ X ⊆ R n (X compact and convex) and ∈ E ⊆ R k (E open). We assume thatx is a fixed point of (6.22) if = 0 and that the Jacobian f (x, 0 ) of f evaluated at (x, 0 ) exists and is continuous. Furthermore, we posit thatx is a hyperbolic equilibrium. If we define the function Hence, for every ∈ U, i.e., for close to 0 , (6.22) has a uniquely determined fixed pointx = φ( ) close tox. With the help of the Hartman-Grobman theorem, it can be shown that the stability properties of the perturbed fixed pointsx are the same as those of the unperturbed,x. The reason is that if an equilibrium is hyperbolic, this property persists under small perturbations (the Jacobian changes continuously if parameters change continuously); see also Theorem 4.4 in [57]. The extension of this argument to finitely many hyperbolic equilibria is evident.
Although hyperbolic equilibria change continuously under small perturbations, limit sets of trajectories do not: they can explode. Thus, perturbations could introduce 'new' limit sets away from the hyperbolic equilibria. What has good properties under perturbations is the set of chain-recurrent points [25].
Let X be a compact set with metric d and let f : X → X be a continuous map. A point x ∈ X is called chain recurrent (with respect to f ) if, for every δ > 0, there exists a finite sequence x 0 = x, x 1 , . . . , x r−1 , x r = x (often called a δ-pseudo-orbit) such that d(f (x m ), x m+1 ) < δ for m = 0, 1, . . . , r − 1. The set of chain-recurrent points contains the limit sets of all orbits and behaves well under perturbations ( [4], p. 244).
Outline of the proof of Theorem 6.5. (a) and (b) follow from the implicit function theorem and the Hartman-Grobman theorem. That asymptotically stable boundary equilibria remain in the state space follows from Brouwer's fixed point theorem (for details, see [57], especially Theorem 4.4).
(c) Let F = {p ∈ S Γ I : p i,α (w i,α −w α ) = 0 ∀i ∈ I, ∀α ∈ Γ} (6.25) denote the set of equilibria of (6.21). Then, within each deme, ∆w α ≥ 0 with equality only at equilibrium; cf. (2.13). Hence, satisfies ∆w(p) ≥ 0 with ∆w(p) = 0 if and only if p ∈ F. Because, on F,w takes only finitely many values, Theorem 3.14 in [4] implies that the chain-recurrent points of (6.21) are exactly the equilibria. Now we can follow the proof of Theorem 2.3 in [77] almost verbally. Because the set of chain-recurrent points consists only of hyperbolic equilibria, this is also true for small C 1 perturbations of the dynamics ( [4], p. 244). Indeed, as an immediate consequence of the definition of chain recurrence, it follows that the chain-recurrent set of (6.21) changes in an upper semicontinuous way with . In particular, the chain-recurrent set for > 0 is contained in the union of the δ-neighborhoods of the equilibria for = 0, with δ → 0 for → 0. By the implicit function theorem and the openness of hyperbolicity (Hartman-Grobman theorem), if > 0 is small, then the maximal invariant sets in those neighborhoods are hyperbolic equilibria. Hence, for small , the chain-recurrent set consists only of finitely many equilibria, which implies convergence of all trajectories. Remark 6.6. Boundary equilibria that are unstable in the absence of migration, can disappear under weak migration because they may leave S Γ I . A simple example is overdominance in two diallelic demes: the unstable zero-migration equilibria (p 1,1 , p 1,2 ) = (1, 0) and (p 1,1 , p 1,2 ) = (0, 1) do not survive perturbation. Indeed, if > 0, the perturbation of (1, 0) must have the form (1 − z 1 , z 2 ). If the fitnesses of the genotypes are as in (4.1) and the migration rates are m 12 = m 1 and m 21 = m 2 , then straightforward calculations show that this equilibrium is given by , if x α < 1 and y α < 1.
As another application of Theorem 6.5, we study the number of alleles that can be maintained at an asymptotically stable equilibrium under weak migration. First we prove a simple result for a single deme. For more general results, see [79]. Proposition 6.8. Let Γ = 1 and assume that the alleles are ordered such that w ii ≥ w i+1,i+1 for i = 1, . . . , I − 1. In addition, assume that w 11 > w 22 and there is intermediate dominance, i.e., for every i and every j > i. Then allele A 1 is fixed as t → ∞. The corresponding equilibrium is globally asymptotically stable.
Proof. The assumptions imply w 1i ≥ w ii ≥ w Ii for every i ∈ I and w 11 > w I1 or w I1 > w II . It follows that w 1 = i w 1i p i ≥ i w Ii p i = w I and, if p 1 = 0 and p I = 0, w 1 > w I . (6.28) Therefore, we obtain if p 1 > 0 and p I > 0. Hence, p I (t) → 0 as t → ∞ provided p I < 1 and p 1 > 0. Now we can repeat this argument with p 1 and p I−1 . Proceeding inductively, we obtain p 1 (t) → 1 as t → ∞.
Theorem 6.9. Suppose that migration is sufficiently weak and there is partial dominance in every deme (5.18).
(a) Generically, there is global convergence to an asymptotically stable equilibrium at which at most Γ alleles are present. Thus, the number of demes is a generic upper bound for the number of alleles that can be maintained at a stable equilibrium.
(b) If Γ ≤ I, then there is an open set of parameters such that Γ alleles are segregating at a globally asymptotically stable equilibrium.
Proof. Proposition 6.8 shows that, in the absence of migration, one allele is fixed in every deme. The parameter combinations satisfying the assumptions of the theorem clearly form an open set of all possible parameter combinations. Therefore, without migration, at most Γ alleles can be maintained at an asymptotically stable equilibrium, and this can be achieved on an open set in the parameters space. Moreover, this equilibrium is globally asymptotically stable. Therefore, Theorem 6.5 yields statement (a) for weak migration. Clearly, the same set of alleles as without migration occurs at this equilibrium.
If Γ ≤ I, we still obtain an open set of parameters if we choose fitnesses such that in each deme a different allele has the highest homozygous fitness. Hence, the upper bound Γ can be achieved on an open set, and Theorem 6.5 yields statement (b). Therefore, in contrast to the Levene model (Theorem 5.8), in which migration is strong, for weak migration the number of alleles that can be maintained at a stable equilibrium cannot exceed the number of demes.
6.4. Strong migration. If migration is much stronger than selection, we expect that rapid convergence to spatial quasi-homogeneity occurs. After this short initial phase, evolution should be approximately panmictic with suitably averaged allele frequencies. Because in the absence of selection, there exists a globally attracting manifold of equilibria, so that the dynamics is not gradient like, the derivation of perturbation results is much more delicate than for weak migration. Since the fundamental ideas in the proofs of the most relevant results are essentially the same as if selection acts on many loci, we defer the analysis of strong migration to Section 7.6, where the multilocus case is treated. 6.5. Uniform selection. Selection is called uniform if w ij,α = w ij for every i, j, and every α. In the following we state sufficient conditions under which there is no genetic indication of population structure. We callp ∈ S Γ I a uniform selection equilibrium if everyp ·,α is an equilibrium of the pure selection dynamics (6.21) andp ·,α =p ·,β for every α, β. Theorem 6.10 ([81], Theorem 5.1). Ifp ∈ S Γ I is a uniform selection equilibrium, thenp is an equilibrium of the (full) migration-selection dynamics (3.5), andp is either asymptotically stable for both (6.21) and (3.5), or unstable for both systems.
It can also be shown that the ultimate rate of convergence to equilibrium is determined entirely by selection and is independent of migration.
Next, one may ask for the conditions when the solutions of (3.5) converge globally to a uniform selection equilibrium. One may expect global convergence under migration if the uniform selection equilibrium is globally asymptotically stable without migration. So far, only weaker results could be proved. (a) Migration is weak and, without migration, every equilibrium is hyperbolic. (b) Migration is strong, M is ergodic, and every equilibrium of the strong-migration limit (Section 7.6) is hyperbolic.
7. Multilocus models. Since many phenotypic traits are determined by many genes, many loci are subject to selection. If loci are on the same chromosome, especially if they are within a short physical distance, they can not be treated independently. In order to understand the evolutionary effects of selection on multiple loci, models have to be developed and studied that take linkage and recombination into account. Because of the complexity of the general model, useful analytical results can be obtained essentially only for limiting cases or under special assumptions. Whereas the former can be sometimes extended using perturbation theory to obtain general insight, the latter are mainly useful to study specific biological questions or to demonstrate the kind of complexity that can arise.
Before developing the general model, we introduce the basic model describing the interaction of selection and recombination and point out some of its fundamental properties. We begin by illustrating the effects of recombination in the simplest meaningful setting. 7.1. Recombination between two loci. We consider two loci, A and B, each with two alleles, A 1 , A 2 , and B 1 , B 2 . Therefore, there are four possible gametes, Genes on different chromosomes are separated during meiosis with probability one half (Mendel's Principle of Independent Assortment). If the loci are on the same chromosome, they may become separated by a recombination event (a crossover) between them. We denote this recombination probability by r. The value of r usually depends on the distance between the two loci along the chromosome. Loci with r = 0 are called completely linked (and may be treated as a single locus), and loci with r = 1 2 are called unlinked. The maximum value of r = 1 2 occurs for loci on different chromosomes, because then all four gametes are produced with equal frequency 1 4 . Thus, the recombination rate satisfies 0 ≤ r ≤ 1 2 . If, for instance, in the initial generation only the genotypes A 1 B 1 /A 1 B 1 and A 2 B 2 /A 2 B 2 are present, then in the next generation only these double homozygotes, as well as the two double heterozygotes A 1 B 1 /A 2 B 2 and A 1 B 2 /A 2 B 1 will be present. After further generations of random mating, all other genotypes will occur, but not immediately at their equilibrium frequencies. The formation of gametic types other than A 1 B 1 or A 2 B 2 requires that recombination occurs between the two loci.
We denote the frequency of gamete A i B j by P ij and, at first, admit an arbitrary number of alleles at each locus. Let the frequencies of the alleles A i at the first locus be denoted by p i and those of the alleles B j at the second locus by q j . Then The allele frequencies are no longer sufficient to describe the genetic composition of the population because, in general, they do not evolve independently. Linkage equilibrium (LE) is defined as the state in which holds for every i and j. Otherwise, the population is in linkage disequilibrium (LD). LD is equivalent to probabilistic dependence of allele frequencies accross loci. Given P ij , we want to find the gametic frequencies P ij in the next generation after random mating. The derivation of the recursion equation is based on the following basic fact of Mendelian genetics: an individual with genotype A i B j /A k B produces gametes of parental type if no recombination occurs (with probability 1 − r), and recombinant gametes if recombination occurs (with probability r). Therefore, the fraction of gametes A i B j and A k B is 1 2 (1 − r) each, and that of A i B and A k B j is 1 2 r each. From these considerations, we see that the frequency of gametes of type A i B j in generation t + 1 produced without recombination is (1 − r)P ij , and that produced with recombination is rp i q j because of random mating. Thus, This shows that the gene frequencies are conserved, but the gamete frequencies are not, unless the population is in LE. Commonly, LD between alleles A i and B j is measured by the parameter The D ij are called linkage disequilibria. From (7.3) and (7.4) we infer and Therefore, unless r = 0, linkage disequilibria decay at the geometric rate 1 − r and LE is approached gradually without oscillation.
For two alleles at each locus, it is more convenient to label the frequencies of the gametes A 1 B 1 , A 1 B 2 , A 2 B 1 , and A 2 B 2 by x 1 , x 2 , x 3 , and x 4 , respectively. A simple calculation reveals that Thus, the recurrence equations (7.3) for the gamete frequencies can be rewritten as The two-locus gametic frequencies are the elements of the simplex S 4 and may be represented geometrically by the points in a tetrahedron. The subset where D = 0 forms a two-dimensional manifold and is called the linkage equilibrium, or Wright, manifold. It is displayed in Figure 7.1.
If r > 0, (7.6) implies that all solutions of (7.9) converge to the LE manifold along straight lines, because the allele frequencies x 1 + x 2 and x 1 + x 3 remain constant, where sets such as x 1 + x 2 = const. represent planes in this geometric picture. The LE manifold is invariant under the dynamics (7.9). 7.2. Two diallelic loci under selection. To introduce selection, we assume that viability selection acts on juveniles. Then recombination and random mating occurs. Since selection acts on diploid individuals, we assign fitnesses to two-locus genotypes. We denote the fitness of genotype ij (i, j ∈ {1, 2, 3, 4}) by w ij , where we assume w ij = w ji , because usually it does not matter which gamete is paternally or maternally inherited. The marginal fitness of gamete i is defined by w ij x j , and the mean fitness of the population isw = 4 i,j=1 w ij x i x j . If we assume, as is frequently the case, that there is no position effect, i.e., w 14 = w 23 , simple considerations yield the selection-recombination dynamics [63]: This is a much more complicated dynamical system than either the pure selection dynamics (2.6) or the pure recombination dynamics (7.9), and has been studied extensively (for a review, see [13], Sects. II.2 and VI.2). In general, mean fitness may decrease and is no longer maximized at an equilibrium. In addition, the existence of stable limit cycles has been established for this discrete-time model ( [46], [49]) as well as for the corresponding continuous-time model ( [2], [3]). Essentially, the demonstration of limit cycles requires that selection coefficients and recombination rates are of similar magnitude.
There is a particularly important special case in which the dynamics is simple. This is the case of no epistasis, or additive fitnesses. Then the fitness of every genotype A i B j /A k B can be expressed as u (1) ik + u (2) j ; see (7.49) for the general definition. In the absence of epistasis, i.e., if (7.49) holds, mean fitnessw is a (strict) Lyapunov function [29]. In addition, a point p is an equilibrium point of (7.10) if and only if it is both a selection equilibrium for each locus and it is in LE ( [69], [77]). In particular, the equilibria are the critical points of mean fitness.
The reason for the increased complexity of two-locus (or multilocus) systems lies not so much in the increased dimensionality but arises mainly from the fact that epistatic selection generates nonrandom associations (LD) among the alleles at different loci. Recombination breaks up these associations to a certain extent but changes gamete frequencies in a complex way. Thus, there are different kinds of interacting nonlinearities arising in the dynamics under selection and recombination. 7.3. The general model. We extend the migration-selection model of Section 3 by assuming that selection acts on a finite number of recombining loci. The treatment follows [14], which was inspired by [76]. We consider a diploid population with discrete, nonoverlapping generations, in which the two sexes need not be distinguished. The population is subdivided into Γ ≥ 1 panmictic colonies (demes) that exchange adult migrants independently of genotype. In each of the demes, selection acts through differential viabilities, which are time and frequency independent. Mutation and random genetic drift are ignored.
The genetic system consists of L ≥ 1 loci and I n ≥ 2 alleles, A is I 1 + · · · + I L . We use letters i, j, ∈ I for gametes, k, n ∈ L for loci, and α, β ∈ G for demes.
Let p i,α = p i,α (t) represent the frequency of gamete i among zygotes in deme α in generation t. We define the following column vectors: Here, p i , p ·,α , and p signify the frequency of gamete i in each deme, the gamete frequencies in deme α, and all gamete frequencies, respectively. We will use analogous notation for other quantities, e.g. for D i,α . The frequency of allele A where the sum runs over all multi-indices i with the kth component fixed as i k . We write for the vector of frequencies of allele A (k) i k in each deme. Let x ij,α and w ij,α denote the frequency and fitness of genotype ij in deme α, respectively. We designate the marginal fitness of gamete i in deme α and the mean fitness of the population in deme α by The life cycle starts with zygotes in Hardy-Weinberg proportions. Selection acts in each deme on the newly born offspring. Then recombination occurs followed by adult migration and random mating within in each deme. This life cycle extends that in Section 3.2. To deduce the general multilocus migration-selection dynamics ( [76]), let x * ij,α = p i,α p j,α w ij,α /w α (7.15a) be the frequency of genotype ij in deme α after selection, and its frequency after recombination. Here, R i,j is the probability that during gametogenesis, paternal haplotypes j and produce a gamete i by recombination. Finally, migration occurs and yields the gamete frequencies in the next generation in each deme: The recurrence equations (7.15) describe the evolution of gamete frequencies under selection on multiple recombining loci and migration. We view (7.15) as a dynamical system on S Γ I . We leave it to the reader to check the obvious fact that the processes of migration and recombination commute.
The complications introduced by recombination are disguised by the terms R i,j which depend on the recombination frequencies among all subsets of loci. To obtain an analytically useful representation of (7.15b), more effort is required.
Let {K, N} be a nontrivial decomposition of L, i.e., K and its complement N = L\K are each proper subsets of L and, therefore, contain at least one locus. (The decompositions {K, N} and {N, K} are identified.) We designate by c K the probability of reassociation of the genes at the loci in K, inherited from one parent, with the genes at the loci in N, inherited from the other. Let where K runs over all (different) decompositions {K, N} of L, denote the total recombination frequency. We designate the recombination frequency between loci k and n, such that k < n, by c kn . It is given by We define This is a measure of LD in gamete i in deme α. A considerable generalization of the arguments leading to (7.3) shows that (7.15b) can be expressed in the following form [74]: If there is only one deme, then (7.15b') provides the recurrence equation describing evolution under selection on multiple recombining loci. It seems worth noting that the full dynamics, (7.15), does not depend on linkage disequilibria between demes. Let Λ 0,α = p ·,α : p i,α = p If there is no position effect, i.e., if w ij,α = w i K j N ,j K i N ;α for every i, j, and K, then D i,α = 0 for every p ·,α ∈ Λ 0,α . Hence, where D ·,α is defined in analogy to (7.11b). In the absence of selection, equality holds in (7.22). We will often need the following assumption: The backward migration matrix M is ergodic, i.e., irreducible and aperiodic. (E) Given irreducibility, the biologically trivial condition that individuals have positive probability of remaining in some deme, i.e., m αα > 0 for some α, suffices for aperiodicity ( [34], p. 426). Because M is a finite matrix, ergodicity is equivalent to primitivity.

Migration and recombination.
We study migration and recombination in the absence of selection, i.e., if w ij,α = 1 for every i, j, α. Then the dynamics (7.15) reduces to where and {K, N} as above (7.16). Here, where i|i K runs over all multi-indices i with the components in K fixed as i K , denotes the marginal frequency of the gamete with components i k fixed for the loci k ∈ K. In vector form, i.e., with the notation (7.11a), (7.26) becomes 29) The following theorem shows that in the absence of selection trajectories quickly approach LE and spatial homogeneity. This theorem generalizes the well known fact that in multilocus systems in which recombination is the only evolutionary force, linkage disequilibria decay to zero at a geometric rate ( [39], [69]). The proof is very technical and uses induction on the number of embedded loci of the full L-locus system. Theorem 7.1 does not hold if the migration matrix is reducible or periodic (Remark 3.3 in [14]).
In the following, we investigate several limiting cases in which useful general results can be proved and perturbation theory allows important extensions. 7.5. Weak selection. We assume that selection is much weaker than recombination and migration. The main result will be that all trajectories converge to an invariant manifold Ψ close to Ψ 0 (7.30), on which there is LE and allele frequencies are deme independent. On Ψ , the dynamics can be described by a small perturbation of a gradient system on Ψ 0 . This implies that all trajectories converge, i.e., no cycling can occur, and the equilibrium structure can be inferred from that of the much simpler gradient system.
Throughout this section, we assume (E), i.e., the backward migration matrix is ergodic. We assume that there are constants r ij,α ∈ R such that w ij,α = 1 + r ij,α , (7.31) where ≥ 0 is sufficiently small and |r ij | ≤ 1. Migration and recombination rates, m αβ and c K , are fixed so that fitness differences are small compared with them. From (7.14) and (7.31), we deduce in which r i,α (p ·,α ) = j r ij,α p j,α ,r α (p ·,α ) = i,j r ij,α p i,α p j,α .
When selection is dominated by migration and recombination, we expect that linkage disequilibria within demes as well as gamete-and gene-frequency differences between demes decay rapidly to small quantities. In particular, we expect approximately panmictic evolution of suitably averaged gamete frequencies in 'quasi-linkage equilibrium'. We also show that all trajectories converge to an equilibrium point, i.e., no complicated dynamics, such as cycling, can occur. In the absence of migration, this was proved in [77], Theorem 3.1. For a single locus under selection and strong migration, this is the content of Theorem 4.5 in [81]. Theorem 7.2 and its proof combine and extend these results as well as some of the underlying ideas and methods.
To formulate and prove this theorem, we define the vector 1,α , . . . , p of all averaged allele frequencies at every locus. We note that in the presence of selection the P (k) i k , hence π, are time dependent. Instead of p, we will use π, D, and q to analyze (7.15), and occasionally write p = (π, D, q).
On the LE manifold Λ 0,α (7.20), which is characterized by the ρ α (α ∈ G), the selection coefficients of gamete i, allele i n at locus n, and of the subpopulation are cf. (7.33). As in (7.23), let ξ denote the principal left eigenvector of M . We introduce the average selection coefficients of genotype ij, gamete i, allele i n at locus n, and of the entire population: in,α (π) , (7.37c) Forω, we obtain the alternative representations For reasons that will be justified by the following theorem, we call the differential equationṖ in (π) −ω(π) , (7.39a) D = 0 , q = 0 (7.39b) on S Γ I the weak-selection limit of (7.15). In view of the following theorem, it is more convenient to consider (7.39a) and (7.39b) on S Γ I instead of (7.39a) on S I1 ×· · ·×S I L . The differential equation (7.39a) is a Svirezhev-Shashahani gradient (Remark 2.3) with potential functionω. In particular,ω increases strictly along nonconstant solutions of (7.39a) becausė We will also need the assumption All equilibria of (7.39a) are hyperbolic.
(H) 14]). Suppose that (7.15), (7.31), (7.18), (E) and (H) hold, the backward migration matrix M and all recombination rates c K are fixed, and > 0 is sufficiently small. (a) The set of equilibria Ξ 0 ⊂ S Γ I of (7.39) contains only isolated points, as does the set of equilibria Ξ ⊂ S Γ I of (7.15). As → 0, each equilibrium in Ξ converges to the corresponding equilibrium in Ξ 0 .
(b) In the neighborhood of each equilibrium in Ξ 0 , there exists exactly one equilibrium point in Ξ . The stability of each equilibrium in Ξ is the same as that of the corresponding equilibrium in Ξ 0 ; i.e., each pair is either asymptotically stable or unstable.
(c) Every solution p(t) of (7.15) converges to one of the equilibrium points in Ξ .
The essence of this theorem and its proof is that under weak selection (i) the exact dynamics quickly leads to spatial quasi-homogeneity and quasi-linkage equilibrium, and (ii) after this time, the exact dynamics can be perceived as a perturbation of the weak-selection limit (7.39). The latter is much easier to study because it is formally equivalent to a panmictic one-locus selection dynamics. Theorem 7.2 is a singular perturbation result because in the absence of selection every point on Ψ 0 is an equilibrium.
Parts (a) and (b) of the above theorem follow essentially from Theorem 4.4 in [57] which is an application of the implicit function theorem and the Hartman-Grobman theorem. Part (c) is much stronger and relies, among others, on the notion of chain-recurrent points and their properties under perturbations of the dynamics (see Section 6.3).
Outline of the proof of Theorem 7.2. Theorem 7.1 and the theory of normally hyperbolic manifolds imply that for sufficiently small s, there exists a smooth invariant manifold Ψ close to Ψ 0 , and Ψ is globally attracting for (7.15) at a geometric rate (see [77], and the references there). The manifold Ψ is characterized by an equation of the form (D, q) = ψ(π, ) , (7.41) where ψ is a smooth function of π. Thus, on Ψ , and more generally, for any initial values, after a long time, The next step consists in deriving the recurrence equations in an O( ) neighborhood of Ψ 0 which, in particular, contains Ψ . By applying (7.31) and D(t) = O( ) to (7.15b'), straightforward calculations yield for every α ∈ G (see eq. (3.10) in [77]). By averaging, invoking (7.15), (7.43), and (7.41) one obtains By rescaling time as τ = t and letting → 0, the leading term in (7.44), P in (π) −ω(π) , approximates the gradient system (7.39a). Therefore, we haveω(π ) >ω(π) unless π = π. In particular, the dynamics (7.44) on Ψ 0 is gradient like.
The eigenvalues of the Jacobian of (7.44) have the form 1 + ν + O( 2 ), where ν signifies an eigenvalue of the Jacobian of (7.39a). Therefore, (H) implies that also every equilibrium of (7.44) is hyperbolic.
The rest of the proof is identical to that of Theorem 3.1 in [77]. The heart of its proof is the following. Since solutions of (7.15) are in phase with solutions on the invariant manifold Ψ , it is sufficient to prove convergence of trajectories for initial conditions p ∈ Ψ . With the help of the qualitative theory of numerical approximations, it can be concluded that the chain-recurrent set of the exact dynamics (7.15) on Ψ is a small perturbation of the chain-recurrent set of (7.39a). The latter dynamics, however, is gradient like. Because all equilibria are hyperbolic by assumption, the chain-recurrent set consists exactly of those equilibria. Therefore, the same applies to the chain-recurrent set of (7.15) if > 0 is small.
A final point to check is that unstable boundary equilibria remain in the simplex after perturbation. The reason is that an equilibriump of (7.39) on the boundary of S Γ I satisfiesp i = 0 for every i in some subset of I, and this condition is not altered by migration. The positive components ofp are perturbed within the boundary, where the equilibrium remains.
Our next result concerns the average of the (exact) mean fitnesses over demes: 14]). Suppose the assumptions of Theorem 7.2 apply. If (7.42) holds, π is bounded away from the equilibria of (7.39a), and p is within O( 2 ) of Ψ , then ∆w(p) > 0 .
The time to reach an O( 2 ) neighborhood of Ψ can be estimated (Remark 4.6 in [14] and Remark 4.13 in [81]. Unless migration is very weak or linkage very tight, convergence occurs in an evolutionary short period. It follows that mean fitness is increasing during the long period when most allele-frequency change occurs, i.e., after reaching Ψ and before reaching a small neighborhood of the stable equilibrium.
A key step in the proof is to show that is satisfied if (7.42) holds. Therefore, linkage disequilibria and the measure q of spatial diversity change very slowly on Ψ . This justifies to call states on Ψ spatially quasi-homogeneous and to be in quasi-linkage equilibrium. Theorem 7.2 enables the derivation of the following result about the equilibrium structure of the multilocus migration-selection dynamics (7.15). It establishes that arbitrarily many loci with arbitrarily many alleles can be maintained polymorphic by migration-selection balance: Theorem 7.4. Let L ≥ 1, Γ ≥ 2, I n ≥ 2 for every n ∈ L, and let all recombination rates c K be positive and fixed.
(a) There exists an open set Q of migration and selection parameters, such that for every parameter combination in Q, there is a unique, internal, asymptotically stable equilibrium point. This equilibrium is spatially quasi-homogeneous, is in quasi-linkage equilibrium, and attracts all trajectories with internal initial condition. Furthermore, every trajectory converges to an equilibrium point as t → ∞.
(b) Such an open set, Q , also exists if the set of all fitnesses is restricted to be nonepistatic and to display intermediate dominance.
For a more detailed formulation and a proof, see Theorem 2.2 in [15]. The constructive proof suggests that with increasing number of alleles and loci, the proportion of parameter space that allows for a fully polymorphic equilibrium shrinks rapidly because at every locus some form of overdominance of spatially averaged fitnesses is required. In addition, the proof shows that this set Q exists in the parameter region where migration is strong relative to selection.
Remark 7.5. The set Q in the above theorem does not contain parameter combinations such that all one-locus fitnesses are multiplicative or additive, or more generally satisfy DIDID (5.26). In each of these cases, one gamete becomes fixed as t → ∞ (Proposition 2.6 and Remark 2.7 in [15]). 7.6. Strong migration. Now we assume that selection and recombination are weak relative to migration, i.e., in addition to (7.31), we posit where ≥ 0 is sufficiently small and γ K is defined by this relation. Then every solution p(t) of the full dynamics (7.15) converges to a manifold close to q = 0. On this manifold, the dynamics is approximated by the differential equatioṅ which we augment with q = 0 . (7.48b) This is called the strong-migration limit of (7.15).
In general, it cannot be expected that the asymptotic behavior of solutions of (7.15) under strong migration is governed by (7.48) because its chain-recurrent set does not always consist of finitely many hyperbolic equilibria. The differential equation (7.48a) is the continuous-time version of the discrete-time dynamics (7.15b') which describes evolution in a panmictic population subject to multilocus selection and recombination. Although the recombination term is much simpler than in the corresponding difference equation, the dynamics is not necessarily less complex. Akin ([2], [3]) proved that (7.48) may exhibit stable cycling. Therefore, under strong migration and if selection and recombination are about equally weak, convergence of trajectories of (7.15) will not generally occur, and only local perturbation results can be derived ( [14], Proposition 4.10).
The dynamics (7.48) becomes simple if there is a single locus (then the recombination term in (7.48a) vanishes) or if there is no epistasis, i.e., if fitnesses can be written in the form ω ij = n u (n) injn (7.49) for constants u (n) injn ≥ 0. Essentially, this means that fitnesses can be assigned to each one-locus genotype and there is no nonlinear interaction between loci. Then mean fitness is a Lyapunov function, all equilibria are in LE, and global convergence of trajectories occurs generically ( [29], [69], [77]).
Remark 7.6. If there is only a single locus or there is no epistasis, an analog of Theorem 7.2 holds for strong migration, i.e., for sufficiently strong migration the global dynamics of (7.15) is a perturbation of (7.48). In particular, all trajectories converge to an equilibrium close to one of (7.48). For the one-locus case, this is Theorem 4.5 in [81]. For two loci, see Proposition 4.3 in [1]. The proof given there is general. 7.7. Weak recombination. If recombination is weak relative to migration and selection, the limiting dynamics is formally equivalent to a multiallelic one-locus migration-selection model, as treated in Sections 3 -6. Then local perturbation results can be applied but, in general, global convergence of trajectories cannot be concluded and, in fact, does not always occur. 7.8. Weak migration. In contrast to the one-locus case treated in Section 6.3, in the multilocus case the assumption of weak migration is insufficient to guarantee convergence of trajectories. The reason is that in the absence of migration the dynamics (7.15) reduces to for every i ∈ I and every α ∈ G; cf. (7.15b'). Therefore, we have Γ decoupled multilocus selection dynamics, one for each deme. Because already for a single deme, stable cycling has been established ( [46], [49]), global perturbation results can not be achieved without additional assumptions. Such an assumption is weak epistasis. We say there is weak epistasis if we can assign (constant) fitness components u (n) injn,α > 0 to single-locus genotypes, such that where the numbers s ij,α satisfy |s ij,α | ≤ 1 and η ≥ 0 is a measure of the strength of epistasis. It is always assumed that η is small enough, so that w ij,α > 0. For the rest of this section, we assume weak migration, i.e., (6.19) and (6.20), weak epistasis (7.51), and η = η( ) , (7.52) where η : [0, 1) → [0, ∞) is C 1 and satisfies η(0) = 0. Therefore, migration and epistasis need not be 'equally' weak. In particular, no epistasis (η ≡ 0) is included.
In the absence of epistasis and migration ( = η = 0), p is an equilibrium point of (7.50) if and only if for every α ∈ G, p ·,α is both a selection equilibrium for each locus and is in LE ( [69], [77]). In addition, the only chain-recurrent points of (7.50) are its equilibria (Lemmas 2.1 and 2.2 in [77]). Therefore, one can apply the proof of Theorem 2.3 in [77] to deduce the following result, which simultaneously generalizes Theorem 2.3 in [77] and Theorem 6.5: Theorem 5.4). Suppose that in the absence of epistasis every equilibrium of (7.50) is hyperbolic, and > 0 is sufficiently small.
(a) The set of equilibria Σ 0 ⊂ S Γ I of (6.21) contains only isolated points, as does the set of equilibria Σ ⊂ S Γ I of (7.15). As → 0, each equilibrium in Σ converges to the corresponding equilibrium in Σ 0 .
(b) In the neighborhood of each asymptotically stable equilibrium in Σ 0 , there exists exactly one equilibrium in Σ , and it is asymptotically stable. In the neighborhood of each unstable internal equilibrium in Σ 0 , there exists exactly one equilibrium in Σ , and it is unstable. In the neighborhood of each unstable boundary equilibrium in Σ 0 , there exists at most one equilibrium in Σ , and if it exists, it is unstable.
(c) Every solution p(t) of (7.15) converges to one of the equilibrium points in Σ .
In contrast to the case of weak selection (Theorem 7.2), unstable boundary equilibria can leave the state space under weak migration [56]. For an example in the one-locus setting, see Remark 4.2 in [81]. If = 0, then all equilibria are in Λ 0 , i.e., there is LE within each deme, but not between demes. If > 0 is sufficiently small, then there is weak LD within each deme, i.e., D i,α = O( ) for every i and every α.
Remark 7.8. If in Theorem 7.7 the assumption of weak epistasis is not made, statements (a) and (b) remain valid, but (c) does not necessarily follow. Statement (c) holds for arbitrary epistasis if, in the absence of migration, the chain recurrent set consists precisely of the hyperbolic equilibria. This is the case if, for instance, the unperturbed system is a gradient system or if it has a globally asymptotically stable equilibrium [4].
One of the consequences of the above theorem is the following, which contrasts sharply with Theorem 7.4: Theorem 7.9 ([15], Theorem 3.9). For an arbitrary number of loci, sufficiently weak migration and epistasis, and partial dominance, the number of demes is the generic maximum for the number of alleles that can be maintained at any locus at any equilibrium (stable or not) of (7.15). 7.9. Strong recombination. If recombination is much stronger than selection and migration, we expect evolution in quasi-linkage equilibrium. By scaling fitnesses as in (7.31) and migration rates as in (6.19) and (6.20), a procedure analogous to that in the proof of Theorem 7.2 yields the following limiting dynamics of allele frequencies, which is defined on the LE-manifold Λ 0 , d dt p where r (n) in,α (ρ α ) andr α (ρ α ) are defined in (7.36). If the chain-recurrent points of (7.53) are precisely the hyperbolic equilibria, a perturbation result analogous to Theorem 7.2 applies. Unless epistasis is absent, this system is more complicated than L independent one-locus migration-selection systems because the equations for different loci are coupled by the mean fitnesses. 7.10. Weak evolutionary forces. If all evolutionary forces are weak, i.e., if → 0 in (7.31), (6.19), and in c K = γ K for all subsets K ⊆ L , (7.54) the limiting dynamics of (7.15) on S Γ I becomeṡ where c α denotes a relative deme size (and not a recombination rate). To deduce the main results, it is not necessary to express the probabilities R i,j of gamete production in terms of the recombination probabilities c K among sets of loci. Many of the results on the one-locus Levene model in Sections 5.1 and 5.2 can be generalized to the multilocus Levene model if epistasis is absent or weak. These generalizations are based on two key results. The first is that in the absence of epistasis geometric mean fitness is again a Lyapunov function, and the internal equilibria are its stationary points (Theorems 3.2 and 3.3 in [76]). The second key result is that in the absence of epistasis, generically, every trajectory converges to an equilibrium point that is in LE (Theorem 3.1 in [17]). This is true for the haploid and the diploid model. Then a proof analogous to that of Theorem 7.7 yields a global perturbation result for weak epistasis, in particular, generic global convergence to an equilibrium point in quasi-linkage equilibrium (Theorem 7.2 in [17]).
With the aid of these results, the next two theorems can be derived about the maintenance of multilocus polymorphism in the Levene model.   (c 1 , . . . , c Γ−1 ), such that for every parameter combination in W, there is a unique, internal, asymptotically stable equilibrium point of (7.56). This equilibrium is in quasi-linkage equilibrium and globally attracting. The set W is independent of the choice of the recombination rates.
Whereas Theorem 7.10 suggests that spatially varying selection has the potential to maintain considerable multilocus polymorphism, Theorem 7.11 shows that properties of the genetic architecture of the trait (no dominance or DIDID) may greatly constrain this possibility. Numerical results in [16] quantify how frequent polymorphism is in the two-locus two-allele Levene model, and how this depends on the selection scheme and dominance relations. With increasing number of loci or number of alleles per locus, the volume of parameter space in which a fully polymorphic equilibrium is maintained will certainly decrease rapidly.
Other aspects of the multilocus Levene model have also been investigated. Zhivotovsky et al. [104] employed a multilocus Levene model to study the evolution of phenotypic plasticity. Wiehe and Slatkin [101] explored a haploid Levene model in which LD is caused by epistasis. More recently, van Doorn and Dieckmann [97] performed a numerical study that admits new mutations (of variable effect) and substitutions at many loci. They showed that genetic variation becomes increasingly concentrated on a few loci. Roze and Rousset [88] derived recursions for the allele frequencies and for various types of genetic associations in a multilocus infinite-island model. Barton [8] studied certain issues related to speciation using a generalized haploid multilocus Levene model that admits habitat preferences. 7.12. Some conclusions. Local adaptation of subpopulations to their specific environment occurs only if the alleles or genotypes that perform well are maintained within the population. Similarly, differentiation between subpopulations occurs only if they differ in allele or gamete frequencies. Consequently, maintenance of polymorphism is indispensable for evolving or maintaining local adaptation and differentiation. The more loci and alleles are involved, the higher is the potential for local adaptation and differentation. This is the main reason, why conditions for the maintenance of polymorphism played such an important role in this survey.
For a single randomly mating population, polymorphic equilibria do not exist in the absence of epistasis and of overdominance or underdominance. Therefore, in order to study the potential of migration in a heterogeneous environment for maintaining genetic variation, it is natural to assume that in every subpopulation epistasis is absent (or weak) and dominance is intermediate. In addition, intermediate dominance seems to be by far the most common form of dominance in nature (e.g. [21]).
We briefly recall the most relevant results concerning maintenance of genetic variation and put them into perspective. Particularly relevant results for a single locus are provided in Sections 4.2 -4.5, and by Theorems 5.9, 5.16, 5.19, 6.1, and 6.9. Theorems 7.4, 7.9, 7.10, and 7.11 treat multiple loci. Theorems 7.4 and 7.10 establish that in structured populations and in the Levene model, respectively, spatially varying selection can maintain multiallelic polymorphism at arbitrarily many loci under conditions for which in a panmictic population polymorphism can not be maintained.
Interestingly, for strong migration and only two demes, an arbitrary number of alleles can be maintained at each of arbitrarily many loci, whereas for weak migration, the number of demes is a generic upper bound for the number of alleles that can be maintained (Theorem 7.9). Given the widespread opinion that it is easier to maintain polymorphism under weak migration than under strong migration, this is counterintuitive.
This opinion derives from the fact that for a single diallelic locus and two demes, the parameter region for which a protected polymorphism exists increases with decreasing migration rate, provided there is a single parameter that measures the strength of migration, as is the case in the Deakin model; see condition (4.19). Additional support for this opinion comes from the study of weak migration in homogeneous and heterogeneous environments ( [56], [57], [24]), from the analysis of the general Deakin model and the result on stronger mixing (Section 4.3), as well as from various numerical studies (e.g. [92], [93], [94]). However, as already noted in Section 4.3, this is not generally true. In fact, if only one of the migration rates in (4.13) is reduced, the protection condition (4.14) implies that it depends on the sign of the corresponding selection coefficient if protection is facilitated or not.
Theorem 7.4 is not at variance with Theorem 7.9 or the widespread opinion expressed above. It is complimentary, and its proof yields deeper insight into the conditions under which genetic variation can be maintained by migration-selection balance. With strong migration, a stable multiallelic polymorphism requires some form of overdominance of suitably averaged fitnesses for the loci maintained polymorphic. The reason is that strong migration leads to strong mixing, so that gamete and allele frequencies become similar among demes. Because the fraction of volume of parameter space, in which average overdominance holds for multiple alleles, shrinks rapidly with increasing number of alleles or loci, we expect very stringent constraints on selection and dominance coefficients for maintaining polymorphism at many loci if migration is strong.
The conditions for maintaining loci polymorphic under weak migration are much less restrictive. With arbitrary intermediate dominance, the proof of Theorem 7.9 shows that those alleles will be maintained that are the fittest in at least one niche. This, obviously, limits the number of alleles that can be maintained.
Also in the Levene model arbitrarily many loci can be maintained polymorphic if epistasis is absent and dominance is intermediate in every deme (Theorem 7.10). However, the maximum number of polymorphic loci depends on the pattern of dominance and the number of demes ( [76], [17], and Theorem 7.11). It is an open problem if for every (ergodic) migration scheme, arbitrarily many loci can be maintained polymorphic in the absence of epistasis, overdominance, and underdominance.
Although the fraction of parameter space, in which the conditions for maintenance of multiallelic multilocus polymorphisms are satisfied, decreases rapidly with increasing number of alleles or loci, stable polymorphism involving small or moderate numbers of alleles or loci does not seem unlikely. What is required if dispersal is weak, is a mosaic of directional selection pressures, and different genes or genotypes that are locally adapted. 8. Application: A two-locus model for the evolution of genetic incompatibilities with gene flow. Here we show how to apply some of the above perturbation results in conjunction with Lyapunov functions and an index theorem to derive the complete equilibrium and stability structure for a two-locus continentisland model, in which epistatic selection may be strong. The model was developed and studied by Bank et al. [6] to explore how much gene flow is needed to inhibit speciation by the accumulation of so-called Dobzhansky-Muller incompatibilities (DMI).
The idea underlying the concept of DMIs is the following. If a population splits into two, for instance because at least one of the subpopulations moves to a different habitat, and the two subpopulations remain separated for a sufficiently long time, in each of them new, locally adapted alleles may emerge that substitute previous alleles. Usually, such substitutions will occur at different loci. For instance, if ab is the original haplotype (gamete), also called wild type, in one population Ab may become fixed (or reach high frequency), whereas in the other population aB may become fixed. If individuals from the derived populations mate, then hybrids of type AB, which occur by recombination between aB and Ab have sometimes reduced fitness, i.e., they are incompatible (to a certain degree). Therefore, the accumulation of DMIs is considered a plausible mechanism for the evolution of so-called intrinsic postzygotic isolation between allopatric (i.e., geographically separated) populations. Since complete separation is an extreme situation, Bank et al. (2012) studied how much gene flow can inhibit the accumulation of such incompatibilities, hence parapatric speciation.
As a simple scenario, in [6] a continent-island model is assumed, i.e., a continental population, which is fixed for one haplotype (here, aB), sends migrants to an island population in which, before the onset of gene flow, Ab was most frequent or fixed. To examine how much gene flow is needed to annihilate differentiation between the two populations, the equilibrium and stability structure was derived. A stable internal equilibrium, at which all four haplotypes are maintained, is called a DMI because, only when it exists, is differentiation between the two populations maintained at both loci. Bank et al. investigated more general biological scenarios than the one outlined above, and they investigated a haploid and a diploid version of their model. Here, we deal only with the haploid case, as it admits a much more complete analysis and is also representative for an important subset of diploid models.
Our goal is to derive the equilibrium and bifurcation structure of the haploid model. We convey the main ideas and only outline most of the proofs. Detailed proofs are given in the Online Supplement S1 of [6]. Also, comprehensive biological motivation and discussion is provided in [6], as well as illuminating figures.
8.1. The haploid DMI model. Let x 1 , x 2 , x 3 , and x 4 denote the frequencies of the four haplotypes ab, aB, Ab, and AB on the island. They satisfy x i ≥ 0 for every i and 4 i=1 x i = 1. We assume that selection acts on individuals during the haploid phase of their life cycle according to   Table 8.1 is entirely general, with the fitness of the wild type normalized to 0; cf. Section 2.4. The parameters α and β measure a potential selective advantage of the island and continental haplotypes, respectively, on the island. They can be positive or negative. However, further below we will assume that α > 0 because, as will be shown, otherwise a DMI does not exist. Finally, the epistasis parameter γ measures the strength of the incompatibility among the A and B alleles. We assume γ ≥ 0, so that epistasis is negative or absent (γ = 0).
The continuous-time haplotype dynamics is given bẏ where r is the recombination rate between the two loci, D = x 1 x 4 − x 2 x 3 is the LD, η 1 = η 4 = −η 2 = −η 3 = 1, and δ i2 = 1 if i = 2 and δ i2 = 0 otherwise; cf. (7.10) and (7.55). With the fitnesses of Table 8.1, we obtaiṅ This is a dynamical system on the simplex S 4 .We always assume m ≥ 0, r ≥ 0, and γ ≥ 0. For many purposes, it is convenient to describe the dynamics in terms of the allele frequencies p = x 3 + x 4 of A, q = x 2 + x 4 of B, and the measure D of LD. Then the dynamical equations reaḋ It is straightforward to show that the condition x ∈ S 4 translates to 0 ≤ p ≤ 1, 0 ≤ q ≤ 1, and which requires r > α.
(b) If M 2 is asymptotically stable, then S A and S B are not admissible.
(c) As a function of m, boundary equilibria change stability at most once. If a change in stability occurs, then it is from unstable to stable (as m increases).
Finally, if r = 0, there is a fully polymorphic equilibrium R 0 on the edge x 1 = x 4 = 0 of S 4 . Thus, only the island and the continental haplotypes are present. Its coordinates are easily calculated.
Below we derive global asymptotic stability of boundary equilibria for various sets of parameters by applying the theory of Lyapunov functions (e.g. [60], in particular, Theorem 6.4 and Corollary 6.5). By global asymptotic stability of an equilibrium we mean that (in addition to stability) every trajectory, such that initially all alleles are present, converges to this equilibrium. By Lemma 8.1 there is at most one asymptotically stable boundary equilibrium for any given set of parameters. Hence, convergence of all trajectories to the boundary is sufficient for demonstrating global stability. Because global convergence to the boundary precludes the existence of an internal equilibrium, it yields necessary conditions for a DMI. Condition (8.8) means that w(Ab) > max{w(ab), w(AB)}, i.e., the island type has higher fitness than its one-step mutational neighbors.
Condition (8.10) follows in a similar way. Proof. We define where x 1 > 0 and x 4 > 0 is assumed. We note that Z = 1 if and only if D = 0, and Z < 1 if and only if D > 0. Theṅ the case m = γ = 0. We leave the simple first two cases to the reader. The third case is also not difficult and follows immediately from Section 3.4.1 in [19].
From here on, we assume α > 0 and γ > β and r > β − α , Our next aim is the derivation of a cubic equation from which the coordinate p of an internal equilibrium (p, q, D) can be obtained. Given p, the coordinates q and D can be computed from relatively simple explicit formulas.
By solvingṗ = 0, we find that, for given p and q, and if p = 1 − β/γ, the value of LD at equilibrium is Substituting this into (8.3b), assuming β = 0, and solvingq = 0 for q, we obtain where needs to be nonnegative to yield an admissible equilibrium. Finally, we substitute q = q 1 (p) and D = D(p, q 1 (p)) into (8.3c) and obtain that an equilibrium value p must solve the equation where If we substitute q = q 2 (p) and D = D(p, q 2 (p)) into (8.3c), we obtain instead of (8.21). By simple additional considerations, it is shown that solutions p of (8.21) or (8.23) satisfy where Because p = 1 − β/γ never gives an equilibrium of (8.3) and p = 1 − m/(γ − α) can give rise only to a single-locus polymorphism, any internal equilibrium value p must satisfy P (p) = 0. We can summarize these findings as follows. The admissibility conditions are given by simple formulas [19]. If r = 0, the equilibrium R 0 mentioned in Section 8.2 is obtained from the unique zero of P (p).
For a number of important special or limiting cases, explicit expressions or approximations can be obtained for the internal equilibria. 8.4. Bifurcations of two internal equilibria. A bifurcation of two internal equilibria can occur if and only if P (p * ) = 0, where p * ∈ (0, 1) is a critical point of P , i.e., P (p * ) = 0. There are at most two such critical points, and they are given by Solving either P (p * 1 ) = 0 or P (p * 2 ) = 0 for m, we obtain after some straightforward manipulations that the critical value m * must be a solution of the following quartic equation: [αβ(2β − γ)(α + β − γ + r) + γ(γ − β)rm] 2 [−ψ 1 + 2ψ 2 rm + 27αγ 2 (γ − β)r 2 m 2 ] = 0 , (8.28a) where The zero m = m 0 arising from the first (linear) factor in (8.28a) does not give a valid bifurcation point for internal equilibria.
The second (quadratic) factor in (8.28a) provides two potential solutions. However, because ψ 1 ≥ 0, one is negative. Therefore, the critical value we are looking for is given by m * = 1 27αγ 2 (γ − β)r −ψ 2 + ψ 2 2 + 27αγ 2 (γ − β)ψ 1 . (8.29) At this value, two equilibria with non-zero allele frequencies collide and annihilate each other. Thus, m * is the critical value at which a saddle-node bifurcation occurs. This gives an admissible bifurcation if both equilibria are internal (hence admissible) for either m < m * or m > m * . If p * 1 = p * 2 , i.e., if R = 0, a pitchfork bifurcation could occur at m * . As a function of α, β, and γ, the condition p * 1 = p * 2 can be satisfied at m * only for three different values of r, of which at most two can be positive. It can be shown that at each of these values, one of the emerging zeros of P (p) does not give rise to an admissible equilibrium (because D > 0 there). Thus, only a saddle-node bifurcation can occur.
Simple and instructive series expansions for m * may be found in the Online Supplement of [6]. 8.5. No migration. We assume m = 0. From Section 8.3, we obtain the following properties of internal equilibria (p, q, D). The LD is given by D = D(p, q) = p(1 − p) γq − α β − γ + γp , (8.30) where p = 1 − β/γ; cf. (8.18). For admissibility, (8.26) needs to be satisfied. For given p and if β = 0, the coordinate q of an internal equilibrium can assume only one of the following forms: By Theorem 8.4, for given p, at most one of q 1 = q 1 (p) or q 2 = q 2 (p) can give rise to an equilibrium. The following theorem characterizes the equilibrium and stability structure. (b) Depending on the parameters, the internal equilibrium is given by either (p, q 1 (p), D(p, q 1 (p))) or (p, q 2 (p), D(p, q 2 (p))), where p is one of the two zeros of P (p) = γ 2 r 2 p 2 −r[2β(γ−β)(γ−2α)+rγ 2 ]p+β(γ−β)(α−β−r)(α+β−γ−r) , (8.32) and q i (p) and D(p, q i (p)) are given by This theorem complements the results derived in [33] and [89] for the discretetime, haploid two-locus selection model. In [89], global convergence to a boundary equilibrium was proved if there is no internal equilibrium. If transformed to the parameters used in [89], condition (8.33) yields precisely the cases not covered by Theorem 14 in [89]. Because our model is formulated in continuous time, the internal equilibrium can be determined by solving quadratic equations. This is instrumental for our proof, as is the following index theorem.
Remark 8.7 (Hofbauer's index theorem). Theorem 2 in [48] states the following. For every dissipative semiflow on R n + such that all equilibria are regular, the sum of the indices of all saturated equilibria equals +1.
In our model, the index of a hyperbolic equilibrium is (−1) m , where m is the number of eigenvalues with positive real part. An internal equilibrium is always saturated. If it is asymptotically stable, it has index 1. Equilibria on the boundary of the simplex are saturated if and only if they are externally stable. This is the case if and only if no gamete that is missing at the equilibrium can invade. Because S 4 is attracting within R 4 + , the index of an asymptotically stable (hence, saturated) boundary equilibrium is 1.
Outline of the proof of Theorem 8.6. Statements (a) and (b) follow from four technical lemmas that provide properties of the polynomial P (p) and resulting admissibility conditions for (p, q, D).
The proof of (c) and (d) is the mathematically most interesting part. In view of the above, it remains to prove that the internal equilibrium exists if (8.33) holds and that it is unstable. This follows readily from the index theorem of Hofbauer [48] stated above. In our model, the only boundary equilibria are the four monomorphic states. (e) follows with the help of several simple Lyapunov functions (see [6]).
8.6. Weak migration. Theorem 7.7 and Remark 7.8 in conjunction with Theorem 8.6 and Lemma 8.1 yield the equilibrium configuration for weak migration.
Theorem 8.8. Let m > 0 be sufficiently small. (a) If (8.33) holds, there is one unstable and one asymptotically stable internal equilibrium (the latter is the perturbation of M 3 ). In addition, the monomorphic equilibrium M 2 is asymptotically stable. Neither S A nor S B is admissible.
(b) If (8.33) does not hold and γ = α, β = 0, and r = α − β (so that every equilibrium is hyperbolic if m = 0), the perturbation of the equilibrium M 3 is globally asymptotically stable. The equilibrium M 2 is unstable, and the equilibria S A and S B may be admissible. If S A or S B is admissible, it is unstable.
Throughout, we denote the stable internal equilibrium by I DMI . To first order in m its coordinates are The complete equilibrium and stability structure. Now we are in the position to prove our main results about the equilibrium and bifurcation structure. We continue to assume (8.17), so that a DMI can occur. Throughout, we always consider bifurcations as a function of (increasing) m. We define and note that m A , m B , and m 2 are the critical values of m above which S A , S B , and M 2 , respectively, are asymptotically stable if admissible. We set As the following theorem shows, if m < m − max , an internal equilibrium exists which, presumably, is globally asymptotically stable. Hence, a DMI will evolve from every initial condition. We note that (8.36) defines m − max for all admissible parameter combinations (8.17). We also recall the definition of m * from (8.29).
Theorem 8.9. The following three types of bifurcation patterns occur: Type 1. Here, m − max = 0 and m * > 0. • If 0 < m < m * , there exist two internal equilibria. For sufficiently small m, one is asymptotically stable (I DMI ), the other (I 0 ) is unstable. Presumably, this holds for every m < m * . The monomorphic equilibrium M 2 is asymptotically stable.
• At m = m * , the two internal equilibria collide and annihilate each other by a saddle-node bifurcation. • If m > m * , M 2 is the only equilibrium; it is asymptotically stable and, presumably, globally attracting.
Type 2. Here, 0 < m − max < m * . • If 0 < m < m − max , there is a unique internal equilibrium (I DMI ). For sufficiently small m, it is globally asymptotically stable. Presumably, this holds for every m < m − max . • At m = m − max , an unstable equilibrium (I 0 ) enters the state space by an exchange-of-stability bifurcation with a boundary equilibrium.
• If m − max < m < m * , there are two internal equilibria. Presumably, one is asymptotically stable (I DMI ), the other unstable (I 0 ). One of the boundary equilibria is asymptotically stable.
• At m = m * , the two internal equilibria merge and annihilate each other by a saddle-node bifurcation. • If m > m * , a boundary equilibrium is asymptotically stable and, presumably, globally attracting.
Type 3. Here, m − max > 0 and m * ≤ 0. • If 0 < m < m − max , a unique internal equilibrium (I DMI ) exists. For sufficiently small m, it is globally asymptotically stable. Presumably, this holds for every m < m − max . • At m = m − max , I DMI leaves the state space through a boundary equilibrium by an exchange-of-stability bifurcation.
• If m > m − max , a boundary equilibrium is asymptotically stable and, presumably, globally attracting.
The qualification "presumably" in several statements of the above theorem results from the fact that, except for sufficiently small γ or small r, we cannot exclude that between two established bifurcations in which internal equilibria are involved, an internal equilibrium undergoes two subsequent Hopf bifurcations, where the second reestablishes the stability property before the first. If γ = 0, all eigenvalues are real, and if r = 0 a global Lyapunov function exists. Therefore, these properties extend to sufficiently small γ and r [6].
Proof. Theorem 8.8 provides all equilibrium configurations for small m. Lemma 8.1 provides control over the boundary equilibria. As m increases, they can vanish but not emerge. They can also become asymptotically stable as m increases. For sufficiently large m, there is always a globally asymptotically stable boundary equilibrium. By Theorem 8.4, the number of internal equilibria is at most three. In addition, internal equilibria can emerge or vanish only either by a saddle-node bifurcation (Section 8.4) or because an equilibrium enters or leaves S 4 through one of boundary equilibria, when an exchange of stability occurs. A bifurcation involving the two internal equilibria can occur at most at one value of m, namely at m * (8.29). An exchange-of-stability bifurcation can occur only at the values m A , m B , or m 2 . If it occurs, the respective boundary equilibrium is asymptotically stable for every larger m for which it is admissible. By Hofbauer's [48] index theorem (Remark 8.7), the sum of the indices of all saturated equilibria is 1.
If Case (a) of Theorem 8.8 applies, then M 2 is asymptotically stable for every m > 0, and it is the only boundary equilibrium. Hence, its index is 1. The index of the stable internal equilibrium is also 1. Because the sum of the indices of the internal equilibria must be 0, the index of the unstable equilibrium is -1. Because at most one bifurcation involving the two internal equilibria can occur and because for large m, M 2 is globally asymptotically stable, the bifurcation must be of saddle-node type in which the equilibria collide and annihilate each other (but do not emerge). In principle, the internal equilibria could also leave S 4 through a boundary equilibrium (in this case, it must be M 2 ). However, by the index theorem, they can do so only simultaneously. This occurs if and only if m * = m 2 , which is a non-generic degenerate case. Because the sum of indices of the internal equilibria must be zero, no equilibrium can enter the state space. These considerations settle the bifurcation pattern of Type 1.
If Case (b) of Theorem 8.8 applies, then, for small m, the boundary equilibria are unstable, hence not saturated, and do not contribute to the sum of indices. Since then the indices of internal equilibria must sum up to 1, the only possible bifurcation that does not entail the stability of a boundary equilibrium would be a pitch-fork bifurcation of the internal equilibrium which, by Section 8.4, does not occur in this scenario only if the selection strength on both loci is sufficiently similar and the recombination rate is about as strong as the selection intensity.
In summary, two mechanisms can drive the evolution of parapatric DMIs. In a heterogeneous environment, a DMI can emerge as a by-product of selection against maladapted immigrants. In a homogeneous environment, selection against unfit hybrids can maintain a DMI. No DMI can be maintained if the migration rate exceeds one of the bounds given in Proposition 8.2. Therefore, it is the adaptive advantage of single substitutions rather than the strength of the incompatibility that is the most important factor for the evolution of a DMI with gene flow. In particular, neutral DMIs (α = 0) can not persist in the presence of gene flow. Interestingly, selection against immigrants is most effective if the incompatibility loci are tightly linked, whereas selection against hybrids is most effective if they are loosely linked. Therefore, opposite predictions result concerning the genetic architecture that maximizes the rate of gene flow a DMI can sustain.