A Survey on Migration-Selection Models in Population Genetics
AA Survey on Migration-Selection Models in Population Genetics Reinhard B¨urger
Department of Mathematics, University of Vienna
Contact:
Reinhard B¨urgerInstitut f¨ur MathematikUniversit¨at WienOskar-Morgenstern-Platz 1A-1090 WienAustriaE-mail: [email protected]: +43 1 4277 50631Fax: +43 1 4277 9506
Key words and phrases: differential equations, recurrence equations, stability, conver-gence, perturbation theory, evolution, geographic structure, dispersal, recombination
Running head:
Migration-selection models Acknowledgments.
This survey is based on lecture notes I compiled for a lectures series presentedduring the “Program on Nonlinear Equations in Population Biology” in June 2013 at the Center for PDE,East China Normal University, Shanghai. In turn, these lecture notes were a revised version of noteswritten for a lecture course at the University of Vienna in 2010. I am very grateful to the participantsof these courses for many helpful comments, in particular to Dr. Simon Aeschbacher, Benja Fallenstein,and Dr. Linlin Su. Comments by and useful discussions with Professor Yuan Lou helped to finalize thiswork. I am particularly grateful to Professor Wei-Ming Ni for inviting me to the Center for PDE, andto him and its members for their great hospitality. Financial support by grant P25188 of the AustrianScience Fund FWF is gratefully acknowledged. a r X i v : . [ q - b i o . P E ] S e p bstract This survey focuses on the most important aspects of the mathematical theory of pop-ulation genetic models of selection and migration between discrete niches. Such modelsare most appropriate if the dispersal distance is short compared to the scale at which theenvironment changes, or if the habitat is fragmented. The general goal of such modelsis to study the influence of population subdivision and gene flow among subpopulationson the amount and pattern of genetic variation maintained. Only deterministic mod-els are treated. Because space is discrete, they are formulated in terms of systems ofnonlinear difference or differential equations. A central topic is the exploration of theequilibrium and stability structure under various assumptions on the patterns of selectionand migration. Another important, closely related topic concerns conditions (necessaryor sufficient) for fully polymorphic (internal) equilibria. First, the theory of one-locusmodels with two or multiple alleles is laid out. Then, mostly very recent, developmentsabout multilocus models are presented. Finally, as an application, analysis and results ofan explicit two-locus model emerging from speciation theory are highlighted.
Contents
Population genetics is concerned with the study of the genetic composition of populations.This composition is shaped by selection, mutation, recombination, mating behavior andreproduction, migration, and other genetic, ecological, and evolutionary factors. There-fore, these mechanisms and their interactions and evolutionary consequences are investi-gated. Traditionally, population genetics has been applied to animal and plant breeding,to human genetics, and more recently to ecology and conservation biology. One of themain subjects is the investigation of the mechanisms that generate and maintain geneticvariability in populations, and the study of how this genetic variation, shaped by environ-mental influences, leads to evolutionary change, adaptation, and speciation. Therefore,population genetics provides the basis for understanding the evolutionary processes thathave led to the diversity of life we encounter and admire.Mathematical models and methods have a long history in population genetics, tracingback to Gregor Mendel, who used his education in mathematics and physics to draw hisconclusions. Francis Galton and the biometricians, notably Karl Pearson, developed newstatistical methods to describe the distribution of trait values in populations and to pre-dict their change between generations. Yule (1902), Hardy (1908), and Weinberg (1908)worked out simple, but important, consequences of the particulate mode of inheritanceproposed by Mendel in 1866 that contrasted and challenged the then prevailing blendingtheory of inheritance. However, it was not before 1918 that the synthesis between ge-netics and the theory of evolution through natural selection began to take shape throughFisher’s (1918) work. By the early 1930s, the foundations of modern population geneticshad been laid by the work of Ronald A. Fisher, J.B.S. Haldane, and Sewall Wright. Theyhad demonstrated that the theory of evolution by natural selection, proposed by CharlesDarwin in 1859, can be justified on the basis of genetics as governed by Mendel’s laws. Adetailed account of the history of population genetics is given in Provine (1971).In the following, we explain some basic facts and mechanisms that are needed through-out this survey. Mendel’s prime achievement was the recognition of the particulate nature4f the hereditary determinants, now called genes. Its position along the DNA is calledthe 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 theset of observable properties derived from the genotype is the phenotype . Meiosis is the process of formation of reproductive cells, or gametes (in animals, spermand eggs) from somatic cells. Under
Mendelian segregation , each gamete contains preciselyone of the two alleles of the diploid somatic cell and each gamete is equally likely to containeither one. The separation of the paired alleles from one another and their distribution tothe gametes is called segregation and occurs during meiosis. At mating, two reproductivecells fuse and form a zygote (fertilized egg), which contains the full (diploid) geneticinformation.Any heritable change in the genetic material is called a mutation . Mutations arethe ultimate source of genetic variability and form the raw material upon which selectionacts. 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. Suchmutations occur in many different ways, for instance as base substitutions, in which onenucleotide is replaced by another, as insertions or deletions of DNA, as inversions ofsequences of nucleotides, or as transpositions. For the population-genetic models treatedin this text the molecular origin of a mutant is of no relevance because they assume thatthe relevant alleles are initially present.During meiosis, different chromosomes assort independently and crossing over betweentwo homologous chromosomes may occur. Consequently, the newly formed gamete con-tains 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 allelesat different loci, recombination has the potential to combine favorable alleles of differentancestry in one gamete and to break up combinations of deleterious alleles. These prop-erties are generally considered to confer a substantial evolutionary advantage to sexualspecies relative to asexuals.The mating pattern may have a substantial influence on the evolution of gene frequen-cies. The simplest and most important mode is random mating . This means that matingstake place without regard to ancestry or the genotype under consideration. It seems tooccur frequently in nature. For example, among humans, matings within a populationappear to be random with respect to blood groups and allozyme phenotypes, but arenonrandom with respect to other traits, for example, height.5 election occurs when individuals of different genotype leave different numbers ofprogeny because they differ in their probability to survive to reproductive age ( viabil-ity ), 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 forcefor adaptation and evolution. Since selection affects the entire genome, its consequencesfor the genetic composition of a population may be complex. Selection is measured interms of fitness of individuals, i.e., by the number of progeny contributed to the nextgeneration. There are different measures of fitness, and it consists of several componentsbecause selection may act on each stage of the life cycle.Because many natural populations are geographically structured and selection variesspatially due to heterogeneity in the environment, it is important to study the conse-quences of spatial structure for the evolution of populations. Dispersal of individuals isusually modeled in one of two alternative ways, either by diffusion in space or by mi-gration between discrete niches, or demes. If the population size is sufficiently large, sothat random genetic drift can be ignored, then the first kind of model leads to partialdifferential equations (Fisher 1937, Kolmogorov et al. 1937). This is a natural choice ifgenotype frequencies change continuously along an environmental gradient, as it occursin a cline (Haldane 1948). Here we will not be concerned with this wide and fruitful areaand refer instead to Barton (1999), Nagylaki and Lou (2008), Barton and Turelli (2011),Lou et al. (2013) for recent developments and references.Instead, this survey focuses on models of selection and migration between discretedemes. Such models are most appropriate if the dispersal distance is short comparedto the scale at which the environment changes, or if the habitat is fragmented. Theyoriginated from the work of Haldane (1930) and Wright (1931). Most of the existingtheory has been devoted to study selection on a single locus in populations with discrete,nonoverlapping generations that mate randomly within demes. However, advances inthe theory of multilocus models have been made recently. The general goal is to studythe influence of population subdivision and of gene flow among subpopulations on theamount and pattern of genetic variation maintained. The models are formulated in termsof systems of nonlinear difference or differential equations. The core purpose of thissurvey is the presentation of the methods of their analysis, of the main results, and oftheir biological implications.For mathematically oriented introductions to the much broader field of populationgenetics, we refer to the books of Nagylaki (1992), B¨urger (2000), Ewens (2004), andWakeley (2008). The two latter texts treat stochastic models in detail, an important and6opical area ignored in this survey.
Darwinian evolution is based on selection and inheritance. In this section, we summarizethe essential properties of simple selection models. It will prepare the ground for thesubsequent study of the joint action of spatially varying selection and migration. Proofsand a detailed treatment may be found in Chapter I of B¨urger (2000). Our focus is onthe 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 gen-erations, such as annual plants or insects. We assume two sexes that need not be distin-guished because gene or genotype frequencies are the same in both sexes (as is always thecase in monoecious species). Individuals mate at random with respect to the locus underconsideration, i.e., in proportion to their frequency. We also suppose that the popula-tion 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 orgenotype frequencies can be described by difference or recurrence equations. These as-sumptions reflect an idealized situation which will model evolution at many loci in manypopulations or species, but which is by no means universal.
With the blending theory of inheritance variation in a population declines rapidly, andthis was one of the arguments against Darwin’s theory of evolution. With Mendelianinheritance there is no such dilution of variation, as was shown independently by the fa-mous British mathematician Hardy (1908) and, in much greater generality, by the Germanphysician Weinberg (1908, 1909).We consider a single locus with I possible alleles A i and write I = { , . . . , I } for the setof all alleles. We denote the frequency of the ordered genotype A i A j by P ij , so that thefrequency of the unordered genotype A i A j is P ij + P ji = 2 P ij . Subscripts i and j always7efer to alleles. Then the frequency of allele A i in the population is p i = I (cid:88) j =1 P ij . After one generation of random mating the zygotic proportions satisfy P (cid:48) ij = p i p j for every i and j . A mathematically trivial, but biologically very important, consequence is that (in theabsence of other forces) gene frequencies remain constant across generations, i.e., p (cid:48) i = p i for every i . (2.1)In other words, in a (sufficiently large) randomly mating population reproduction doesnot 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 amongzygotes.Evolutionary mechanisms such as selection, migration, mutation, or random geneticdrift distort Hardy-Weinberg proportions, but Mendelian inheritance restores them ifmating is random. Selection occurs when genotypes in a population differ in their fitnesses, i.e., in theirviability, mating success, or fertility and, therefore, leave different numbers of progeny.The basic mathematical models of selection were developed and investigated in the 1920sand early 1930s by Fisher (1930), Wright (1931), and Haldane (1932).We will be concerned with the evolutionary consequences of selection caused by dif-ferential viabilities, which leads to simpler models than (general) fertility selection (e.g.,Hofbauer and Sigmund 1988, Nagylaki 1992). Suppose that at an autosomal locus the al-leles A , . . . , 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 If no summation range is indicated, it is assumed to be over all admissible values; e.g., (cid:80) i = (cid:80) i ∈ I Unless stated otherwise, a prime, (cid:48) , always signifies the next generation. Thus, instead of P ij ( t ) and P ij ( t + 1), we write P ij and P (cid:48) ij (and analogously for other quantities). P ij are in Hardy-Weinberg proportions. We assume that the fitness(viability) w ij of an A i A j individual is constant, i.e., independent of time, populationsize, 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 P ∗ ij = w ij P ij ¯ w = w ij p i p j ¯ w , where we have used (2.2). Here,¯ w = (cid:88) i,j w ij P ij = (cid:88) i,j w ij p i p j = (cid:88) i w i p i (2.3)is the mean fitness of the population and w i = (cid:88) j w ij p j (2.4)is the marginal fitness of allele A i . Both are functions of p = ( p , . . . , p I ) (cid:62) . Therefore, the frequency of A i after selection is p ∗ i = (cid:88) j P ∗ ij = p i w i ¯ w . (2.5)Because of random mating, the allele frequency p (cid:48) i among zygotes of the next generationis also p ∗ i (2.1), so that allele frequencies evolve according to the selection dynamics p (cid:48) i = p i w i ¯ w , i ∈ I . (2.6)This recurrence equation preserves the relation (cid:88) i p i = 1and describes the evolution of allele frequencies at a single autosomal locus in a diploidpopulation. We view the selection dynamics (2.6) as a (discrete) dynamical system onthe simplex S I = (cid:26) p = ( p , . . . , p I ) (cid:62) ∈ R I : p i ≥ i ∈ I , (cid:88) i p i = 1 (cid:27) . (2.7)Although selection destroys Hardy-Weinberg proportions, random mating re-establishesthem. Therefore, (2.6) is sufficient to study the evolutionary dynamics. Throughout, the superscript (cid:62) denotes vector or matrix transposition. w ij is multiplied by the sameconstant. This is very useful because it allows to rescale the fitness parameters accordingto convenience (also their number is reduced by one). Therefore, we will usually considerrelative fitnesses and not absolute fitnesses.Fitnesses are said to be multiplicative if constants v i exist such that w ij = v i v j (2.8)for every i, j . Then w i = v i ¯ v , where ¯ v = (cid:80) i v i p i , and ¯ w = ¯ v . Therefore, (2.6) simplifiesto p (cid:48) i = p i v i ¯ v , i ∈ I , (2.9)which can be solved explicitly because it is equivalent to the linear system x (cid:48) i = v i x i . It iseasy 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 w ij = v i + v j (2.10)for every i, j . Then w i = v i + ¯ v , where ¯ v = (cid:80) i v i p i , and ¯ w = 2¯ v . Although this assumptionis important (it means absence of dominance; see Sect. 2.3), it does not yield an explicitsolution of the selection dynamics. Example 2.1.
Selection is very efficient. We assume (2.8). Then the solution of (2.9) is p i ( t ) = p i (0) v ti (cid:80) j p j (0) v tj . (2.11)Suppose that there are only two alleles, A and A . If A is the wild type and A is anew beneficial mutant, we may set (without loss of generality!) v = 1 and v = 1 + s .Then we obtain from (2.11): p ( t ) p ( t ) = p (0) p (0) (cid:18) v v (cid:19) t = p (0) p (0) (1 + s ) t . (2.12)Thus, A increases exponentially relative to A .For instance, if s = 0 .
5, then after 10 generations the frequency of A has increased bya factor of (1 + s ) t = 1 . ≈ . A . If s = 0 .
05 and t = 100, this factor is(1 + s ) t = 1 . ≈ .
5. Therefore, slight fitness differences may have a big long-termeffect, in particular, since 100 generations are short on an evolutionary time scale.10n important property of (2.6) is that mean fitness is nondecreasing along trajectories(solutions) , i.e., ¯ w (cid:48) = ¯ w ( p (cid:48) ) ≥ ¯ w ( p ) = ¯ w , (2.13)and equality holds if and only if p is an equilibrium. A particularly elegant proof was provided by Kingman (1961). As noted by Nagylaki(1977), (2.13) follows immediately from an inequality of Baum and Eagon (1967) by notingthat (2.6) can be written as p (cid:48) i = p i ∂ ¯ w∂p i (cid:30)(cid:88) j p j ∂ ¯ w∂p j because ∂ ¯ w/∂p i = 2 w i .The statement (2.13) is closely related to Fisher’s Fundamental Theorem of NaturalSelection , which Fisher (1930) formulated as follows: “The rate of increase in fitness of any organism at any time is equalto its genetic variance in fitness at that time.”
For recent discussion, see Ewens (2011) and B¨urger (2011).In mathematical terms, (2.13) shows that ¯ w is a Lyapunov function. This has anumber of important consequences. For instance, complex dynamical behavior such aslimit cycles or chaos can be excluded. All trajectories approach the set of points p ∈ S I that are maxima of ¯ w . This is a subset of the set of equilibria. From (2.6) it is obviousthat the equilibria are precisely the solutions of p i ( w i − ¯ w ) = 0 for every i ∈ I . (2.14)We call an equilibrium internal, or fully polymorphic, if p i > i (all alleles arepresent). The I equilibria defined by p i = 1 are called monomorphic because only oneallele is present.The following result summarizes a number of important properties of the selectiondynamics. Proofs and references to the original literature may be found in B¨urger (2000,Chap. I.9); see also Lyubich (1992, 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 if ˆ p is a critical point of the restriction of meanfitness ¯ w ( p ) to the minimal subsimplex of S I that contains the positive components of ˆ p . p is called an equilibrium, or fixed point, of the recurrence equation p (cid:48) = f ( p ) if f ( p ) = p . We usethe term equilibrium point to emphasize that we consider an equilibrium that is a single point. The termequilibrium may also refer to a (connected) manifold of equilibrium points. . If the number of equilibria is finite, then it is bounded above by I − .4. An internal equilibrium is asymptotically stable if and only if it is an isolated localmaximum of ¯ w . Moreover, it is isolated if and only if it is hyperbolic (i.e., the Jacobianhas no eigenvalues of modulus 1).5. An equilibrium point is stable if and only if it is a local, not necessarily isolated,maximum of ¯ w .6. If an asymptotically stable internal equilibrium exists, then every orbit starting inthe 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 − alleles will be absent at astable equilibrium.9. Every orbit converges to one of the equilibrium points ( even if stable manifolds ofequilibria exist ) . For the purpose of illustration, we work out the special case of two alleles. We write p and 1 − p instead of p and p . Further, we use relative fitnesses and assume w = 1 , w = 1 − hs , w = 1 − s , (2.15)where s is called the selection coefficient and h describes the degree of dominance. Weassume s > A is called dominant if h = 0, partially dominant if 0 < h < , recessive if h = 1, and partially recessive if < h < No dominance refers to h = . Absence ofdominance is equivalent to additive fitnesses (2.10). If h <
0, there is overdominance or heterozygote advantage . If h >
1, there is underdominance or heterozygote inferiority .From (2.4), the marginal fitnesses of the two alleles are w = 1 − hs + hsp and w = 1 − s + s (1 − h ) p and, from (2.3), the mean fitness is¯ w = 1 − s + 2 s (1 − h ) p − s (1 − h ) p . It is easily verified that the allele-frequency change from one generation to the next can12
Schematic selection dynamics with two alleles p
10 ,10 ≤≤<< hs
10 ,10 ≤≤<< hs
0 ,10 <<< hs
1 ,10 ><< hs Figure 2.1:
Convergence patterns for selection with two alleles. be written as ∆ p = p (cid:48) − p = p (1 − p )2 ¯ w d ¯ w d p (2.16a)= p (1 − p ) s ¯ w [1 − h − (1 − h ) p ] . (2.16b)There exists an internal equilibrium if and only if h < h > p = 1 − h − h . (2.17)If dominance is intermediate, i.e., if 0 ≤ h ≤
1, then (2.16) shows that ∆ p > < p <
1, hence p = 1 is globally asymptotically stable.If h < h >
1, we write (2.16) in the form∆ p = sp (1 − p )¯ w (1 − h )(ˆ p − p ) . (2.18)In the case of overdominance ( h < < sp (1 − p )(1 − h ) / ¯ w < < p < p is globally asymptotically stable and convergence is monotonic. If h >
1, then themonomorphic equilibria p = 0 and p = 1 each are asymptotically stable and ˆ p is unstable.The three possible convergence patterns are shown in Figure 2.1. Figure 2.2 demon-strates that the degree of (intermediate) dominance strongly affects the rate of spread ofan advantageous allele. 13igure 2.2: Selection of a dominant ( h = 0, solid line), intermediate ( h = 1 /
2, dashed), andrecessive ( h = 1, dash-dotted) allele. The initial frequency is p = 0 .
005 and the selectiveadvantage is s = 0 .
05. If the advantageous allele is recessive, its initial rate of increase isvanishingly small because the frequency p of homozygotes is extremely low when p is small.However, only homozygotes are ‘visible’ to selection. Most higher animal species have overlapping generations because birth and death occurcontinuously in time. This, however, may lead to substantial complications if one wishesto derive a continuous-time model from biological principles. By contrast, discrete-timemodels can frequently be derived straightforwardly from simple biological assumptions.If evolutionary forces are weak, a continuous-time version can usually be obtained as anapproximation to the discrete-time model.A rigorous derivation of the differential equations describing gene-frequency changeunder selection in a diploid population with overlapping generations is a formidable taskand requires a complex model involving age structure (see Nagylaki 1992, Chap. 4.10).Here, we simply state the system of differential equations and justify it in an alternativeway.In a continuous-time model, the (
Malthusian ) fitness m ij of a genotype A i A j is definedas its birth rate minus its death rate. Then the marginal fitness of allele A i is m i = (cid:88) j m ij p j , m = (cid:88) i m i p i = (cid:88) i,j m ij p i p j , and the dynamics of allele frequencies becomes˙ p i = d p i d t = p i ( m i − ¯ m ) , i ∈ I . (2.19)This is the analogue of the discrete-time selection dynamics (2.6). Its state space is againthe simplex S I . The equilibria are obtained from the condition ˙ p i = 0 for every i . Wenote that (2.19) is a so-called replicator equation (see Hofbauer and Sigmund 1998).If we set w ij = 1 + sm ij for every i, j ∈ I , (2.20)where s > w i = 1 + sm i and ¯ w = 1 + s ¯ m .Following Nagylaki (1992, p. 99), we approximate the discrete model (2.6) by thecontinuous model (2.19) under the assumption of weak selection, i.e., small s in (2.20).We rescale time according to t = (cid:98) τ /s (cid:99) , where (cid:98) (cid:99) denotes the closest smaller integer.Then s may be interpreted as generation length and, for p i ( t ) satisfying the differenceequation (2.6), we write π i ( τ ) = p i ( t ). Then we obtain formallyd π i d τ = lim s ↓ s [ π i ( τ + s ) − π i ( τ )] = lim s ↓ s [ p i ( t + 1) − p i ( t )] . From (2.6) and (2.20), we obtain p i ( t + 1) − p i ( t ) = sp i ( t )( m i − ¯ m ) / (1 + s ¯ m ). Therefore,˙ π i = π i ( m i − ¯ m ) and ∆ p i ≈ s ˙ π i = sp i ( m i − ¯ m ). We note that (2.6) is essentially the Eulerscheme for (2.19).The exact continuous-time model reduces to (2.19) only if the mathematically incon-sistent assumption is imposed that Hardy-Weinberg proportions apply for every t whichis generally not true. Under weak selection, however, deviations from Hardy-Weinbergdecay to order O ( s ) after a short period of time (Nagylaki 1992).One of the advantages of models in continuous time is that they lead to differentialequations, and usually these are easier to analyze because the formalism of calculus isavailable. An example for this is that, in continuous time, (2.13) simplifies to˙¯ m ≥ , (2.21) Throughout, we use a dot, ˙ , to indicate derivatives with respect to time. much easier to prove than (2.13):˙¯ m = 2 (cid:88) i,j m ij p j ˙ p i = 2 (cid:88) i m i ˙ p i = 2 (cid:88) i ( m i − ¯ m ) p i = 2 (cid:88) i ( m i − ¯ m ) p i . Remark 2.3.
The allele-frequency dynamics (2.19) can be written as a (generalized)gradient system (Svirezhev 1972, Shahshahani 1979):˙ p = G p grad ¯ m = G p (cid:18) ∂ ¯ m∂p , . . . , ∂ ¯ m∂p n (cid:19) (cid:62) . (2.22)Here, G p = ( g ij ) is a quadratic (covariance) matrix, where g ij = Cov( f i , f j ) = p i ( δ ij − p j ) (2.23)and f i ( A k A l ) = k = l = i , if k (cid:54) = l and k = i or l = i , . (2.24)An equivalent formulation is the following covariance form:˙ p i = Cov( f i , m ) , (2.25)where m is interpreted as the random variable m ( A k A l ) = m kl (Li 1967). It holds undermuch more general circumstances (Price 1970, Lessard 1997). We assume a population of diploid organisms with discrete, nonoverlapping generations.This population is subdivided into Γ demes (niches). Viability selection acts within eachdeme and is followed by adult migration (dispersal). After migration random matingoccurs within each deme. We assume that the genotype frequencies are the same inboth sexes (e.g., because the population is monoecious). We also assume that, in everydeme, the population is so large that gene and genotype frequencies may be treated asdeterministic, i.e., we ignore random genetic drift.
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 = { , . . . , Γ } forthe set of all demes. The presentation below is based on Chapter 6.2 of Nagylaki (1992).16e denote the frequency of allele A i in deme α by p i,α . Therefore, we have (cid:88) i p i,α = 1 (3.1)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 indeme α and the mean fitness of the population in deme α are w i,α = (cid:88) j w ij,α p j,α and ¯ w α = (cid:88) i,j w ij,α p i,α p j,α , (3.2)respectively.Next, we describe migration. Let ˜ m αβ denote the probability that an individual in deme α migrates to deme β , and let m αβ denote the probability that an (adult) individual indeme α 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 (cid:88) β ˜ m αβ = 1 and (cid:88) β m αβ = 1 for every α . (3.4)Given the backward migration matrix and the fact that random mating within eachdemes does not change the allele frequencies, the allele frequencies in the next generationare p (cid:48) i,α = (cid:88) β m αβ p ∗ i,β , (3.5a)where p ∗ i,α = p i,α w i,α ¯ w α (3.5b)describes the change due to selection alone; cf. (2.6). These recurrence equations define adynamical system on the Γ-fold Cartesian product S Γ I of the simplex S I . The investigationof 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. Inthe following, we derive their relation to the forward migration rates and discuss conditionswhen selection or migration do not change the deme proportions.17 .2 The relation between forward and backward migration rates To derive this relation, we describe the life cycle explicitly. It starts with zygotes onwhich selection acts (possibly including population regulation). After selection adultsmigrate and usually there is population regulation after migration (for instance becausethe number of nesting places is limited). By assumption, population regulation does notchange genotype frequencies. Finally, there is random mating and reproduction, whichneither changes gene frequencies (Section 2.1) nor deme proportions. The respectiveproportions of zygotes, pre-migration adults, post-migration adults, and post-regulationadults in deme α are c α , c ∗ α , c ∗∗ α , and c (cid:48) α :Zygote (cid:45) Adult (cid:45)
Adult (cid:45)
Adult (cid:45)
Zygote selection migration regulation reproduction c α , p i,α c ∗ α , p ∗ i,α c ∗∗ α , p (cid:48) i,α c (cid:48) α , p (cid:48) i,α c (cid:48) α , p (cid:48) i,α Because no individuals are lost during migration, the following must hold: c ∗∗ β = (cid:88) α c ∗ α ˜ m αβ , (3.6a) c ∗ α = (cid:88) β c ∗∗ β m βα . (3.6b)The (joint) probability that an adult is in deme α and migrates to deme β can be expressedin terms of the forward and backward migration rates as follows: c ∗ α ˜ m αβ = c ∗∗ β m βα . (3.7)Inserting (3.6a) into (3.7), we obtain the desired connection between the forward and thebackward migration rates: m βα = c ∗ α ˜ m αβ (cid:80) γ c ∗ γ ˜ m γβ . (3.8)Therefore, if ˜ M is given, an ansatz for the vector c ∗ = ( c ∗ , . . . , c ∗ Γ ) (cid:62) in terms of c =( c , . . . , c Γ ) (cid:62) is needed to compute M (as well as a hypothesis for the variation, if any, of c ). Two frequently used assumptions are the following.1) Soft selection . This assumes that the fraction of adults in every deme is fixed, i.e., c ∗ α = c α for every α ∈ G . (3.9)18his may be a good approximation if the population is regulated within each deme, e.g.,because individuals compete for resources locally (Dempster 1955).2) Hard selection . Following Dempster (1955), the fraction of adults will be propor-tional to mean fitness in the deme if the total population size is regulated. This has beencalled hard selection and is defined by c ∗ α = c α ¯ w α / ¯ w , (3.10)where ¯ w = (cid:88) α c α ¯ w α (3.11)is the mean fitness of the total population.Essentially, these two assumptions are at the extremes of a broad spectrum of possi-bilities. Soft selection will apply to plants; for animals many schemes are possible.Under soft selection, (3.8) becomes m βα = c α ˜ m αβ (cid:80) γ c γ ˜ m γβ . (3.12)As a consequence, if c is constant ( c (cid:48) = c ), M is constant if and only if ˜ M is constant.If there is no population regulation after migration, then c will generally depend on timebecause (3.6a) yields c (cid:48) = c ∗∗ = ˜ M (cid:62) c . Therefore, the assumption of constant demeproportions, c (cid:48) = c , will usually require that population control occurs after migration.A migration pattern that does not change deme proportions ( c ∗∗ α = c ∗ α ) is called con-servative . Under this assumption, (3.7) yields c ∗ α ˜ m αβ = c ∗ β m βα (3.13)and, by stochasticity of M and ˜ M , we obtain c ∗ β = (cid:88) α c ∗ α ˜ m αβ and c ∗ α = (cid:88) β c ∗ β m βα . (3.14)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 α : c ∗ α ˜ m αβ = c ∗ β ˜ m βα . (3.15)19f 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 areidentical.2) A migration scheme is called doubly stochastic if (cid:88) α ˜ m αβ = 1 for every α . (3.16)If demes are of equal size, then (3.6a) shows that c ∗∗ α = c ∗ α . Hence, with equal deme sizesa doubly stochastic migration pattern is conservative. Under soft selection, deme sizesremain constant without further population regulation. Hence, m αβ = ˜ m βα and M is alsodoubly stochastic.Doubly stochastic migration patterns arise naturally if there is a periodicity, e.g., be-cause the demes are arranged in a circular way. If we posit equal deme sizes and homoge-neous migration, i.e., ˜ m αβ = ˜ m β − α so that migration rates depend only on distance, thenthe backward migration pattern is also homogeneous because m αβ = ˜ m βα = ˜ m α − β and,hence, depends only on β − α . If migration is symmetric, ˜ m αβ = ˜ m βα , and the deme sizesare equal, then dispersion is both reciprocal and doubly stochastic. We introduce three migration patterns that play an important role in the populationgenetics and ecological literature.
Example 3.2.
Random outbreeding and site homing, or the
Deakin (1966) model.
Thismodel assumes that a proportion µ ∈ [0 ,
1] of individuals in each deme leaves their demeand is dispersed randomly across all demes. Thus, they perform outbreeding whereasa proportion 1 − µ remains at their home site. If c ∗∗ α is the proportion (of the totalpopulation) of post-migration adults in deme α , then the forward migration rates aredefined by ˜ m αβ = µc ∗∗ β if α (cid:54) = β , and ˜ m αα = 1 − µ + µc ∗∗ α . If µ = 0, migration is absent;if µ = 1, the Levene model is obtained (see below). Because this migration pattern isreciprocal, M = ˜ M holds.To prove that migration in the Deakin model satisfies (3.15), we employ (3.7) and find m βα = c ∗ α c ∗∗ β ˜ m αβ = µc ∗ α if α (cid:54) = βc ∗ β c ∗∗ β (1 − µ ) + µc ∗ β if α = β . (3.17)20rom this we deduce1 = (cid:88) α m βα = (cid:88) α (cid:54) = β µc ∗ α + c ∗ β c ∗∗ β (1 − µ ) + µc ∗ β = µ · − µ ) c ∗ β c ∗∗ β , (3.18)which immediately yields c ∗∗ β = c ∗ β for every β provided µ <
1. Therefore, we obtain˜ m βα = µc ∗ α if α (cid:54) = β and c ∗ β ˜ m βα = c ∗ β µc ∗ α = c ∗ α ˜ m αβ , i.e., reciprocity.We will always assume soft selection in the Deakin model, i.e., c ∗ α = c α . Thus, for agiven (probability) vector c = ( c , . . . , c Γ ) (cid:62) , the single parameter µ is sufficient to describethe migration pattern: m βα = ˜ m βα = (cid:40) µc α if α (cid:54) = β − µ + µc β if α = β . (3.19) Example 3.3.
The Levene (1953) model assumes soft selection and m αβ = c β . (3.20)Thus, dispersing individuals are distributed randomly across all demes in proportion tothe deme sizes. In particular, migration is independent of the deme of origin and M = ˜ M .Alternatively, the Levene model could be defined by ˜ m αβ = µ β , where µ β > (cid:80) β µ β = 1. Then (3.8) yields m αβ = c ∗ β for every α, β ∈ G . Withsoft selection, we get m αβ = c β . This is all we need if demes are regulated to constantproportions. But the proportions remain constant even without regulation, for (3.6a)gives c (cid:48) α = c ∗∗ α = µ α . This yields the usual interpretation µ α = c α (Nagylaki 1992, Sect.6.3). Example 3.4.
In the linear stepping-stone model the demes are arranged in a linearorder and individuals can reach only one of the neighboring demes. It is an extreme caseamong migration patterns exhibiting isolation by distance , i.e., patterns in which migra-tion diminishes with the distance from the parental deme. In the classical homogeneousversion, the forward migration matrix is˜ M = − m m . . . m − m m m − m m . . . m − m . (3.21)21e leave it to the reader to derive the backward migration matrix using (3.8). It is aspecial case of the following general tridiagonal form: M = n r . . . q n r q Γ − n Γ − r Γ − . . . q Γ n Γ , (3.22)where n α ≥ q α + n α + r α = 1 for every α , q α > α ≥ r α > α ≤ Γ −
1, and q = 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 eachdeme exchanges a fraction m of the population with each of its neighboring demes. Thestepping-stone model has been used as a starting point to derive the partial differentialequations for selection and dispersal in continuous space (Nagylaki 1989). Also circularand infinite variants have been investigated. Juvenile migration is of importance for many marine organisms and plants, where seedsdisperse. It can be treated in a similar way as adult migration. Also models with bothjuvenile and adult migration have been studied. Some authors investigated migrationand selection in dioecious populations, as well as selection on X-linked loci (e.g., Nagylaki1992, pp. 143, 144).Unless stated otherwise, throughout this survey we assume that the backward migra-tion matrix M is constant, as is the case for soft selection if deme proportions and theforward migration matrix are constant. Then the recurrence equations (3.5) provide aself-contained description of the migration-selection dynamics. Hence, they are sufficientto study evolution for an arbitrary number of generations. Of central interest is the identification of conditions that guarantee the maintenance ofgenetic diversity. Often it is impossible to determine the equilibrium structure in detailbecause establishing existence and, even more so, stability or location of polymorphicequilibria is unfeasible. Below we introduce an important concept that is particularlyuseful to establish maintenance of genetic variation at diallelic loci. Throughout thissection we consider a single locus with two alleles. The number of demes, Γ, can bearbitrary. 22 .1 Protected polymorphism
There is a protected polymorphism (Prout 1968) if, independently of the initial conditionsare, a polymorphic population cannot become monomorphic. Essentially, this requiresthat if an allele becomes very rare, its frequency must increase. In general, a protectedpolymorphism is neither necessary nor sufficient for the existence of a stable polymorphicequilibrium. For instance, on the one hand, if there is an internal limit cycle that attractsall solutions, then there is a protected polymorphism. On the other hand, if there aretwo internal equilibria, one asymptotically stable, the other unstable, then selection mayremove one of the alleles if sufficiently rare. A generalization of this concept to multiplealleles would correspond to the concept of permanence often used in ecological models(e.g., Hofbauer and Sigmund 1998).Because we consider only two alleles, we can simplify the notation. We write p α = p ,α for the frequency of allele A in deme α (and 1 − p α for that of A in deme α ). Let p = ( p , . . . , p Γ ) (cid:62) denote the vector of allele frequencies. Instead of using the fitnessassignments w ,α , w ,α , and w ,α , it will be convenient to scale the fitness of the threegenotypes in deme α as follows A A A A A A x α y α (4.1)( x α , y α ≥ x α = w ,α /w ,α and y α = w ,α /w ,α ,provided w ,α > w ,α = 1 − p α + x α p α and ¯ w α = x α p α + 2 p α (1 − p α ) + y α (1 − p α ) , (4.2)and the migration-selection dynamics (3.5) becomes p ∗ α = p α w ,α / ¯ w α (4.3a) p (cid:48) α = (cid:88) β m αβ p ∗ β . (4.3b)We consider this as a (discrete) dynamical system on [0 , Γ .We call allele A protected if it cannot be lost. Thus, it has to increase in frequencyif rare. In mathematical terms this means that the monomorphic equilibrium p = 0 mustbe unstable. To derive a sufficient condition for instability of p = 0, we linearize (4.3)at p = 0. If y α > α (which means that A A is nowhere lethal), a simple23alculation shows that the Jacobian of (4.3a), D = (cid:18) ∂p ∗ α ∂p β (cid:19) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) p =0 , (4.4)is a diagonal matrix with (nonzero) entries d αα = y − α . Because (4.3b) is linear, thelinearization of (4.3) is p (cid:48) = Qp , where Q = M D , (4.5)i.e., q αβ = m αβ /y β .To obtain a simple criterion for protection, we assume that the descendants of indi-viduals 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 isnonnegative. Therefore, the Theorem of Perron and Frobenius (e.g., Seneta 1981) impliesthe existence of a uniquely determined eigenvalue λ > Q such that | λ | ≤ λ holdsfor all eigenvalues of Q . In addition, there exists a strictly positive eigenvector pertainingto λ which, up to multiplicity, is uniquely determined. As a consequence, A is protected if λ > A is not protected if λ < λ = 1, then stability cannot be decided upon linearization). This maximal eigenvaluesatisfies min α (cid:88) β q αβ ≤ λ ≤ max α (cid:88) β q αβ , (4.7)with equality if and only if all the row sums are the same. Example 4.1.
Suppose that A A is at least as fit as A A in every deme and more fit inat least one deme, i.e., y α ≥ α and y β > β . Then q αβ = m αβ /y β ≤ m αβ for every β . Because M is irreducible, there is no β such that m αβ = 0 for every α .Therefore, the row sums (cid:80) β q αβ = (cid:80) β m αβ /y β in (4.7) are not all equal to one, and weobtain λ < max α (cid:88) β q αβ ≤ max α (cid:88) β m αβ = 1 . (4.8)Thus, A is not protected, and this holds independently of the choice of the x α , or w ,α .It can be shown similarly that A is protected if A A is favored over A A in at leastone deme and is nowhere less fit than A A .One obtains the condition for protection of A if, in (4.6), A is replaced by A and λ is the maximal eigenvalue of the matrix with entries m αβ /x β . Clearly, there is a protectedpolymorphism if both alleles are protected.24n the case of complete dominance the eigenvalue condition (4.6) cannot be satisfied.Consider, for instance, protection of A if A is dominant, i.e., y α = 1 for every α . Then q αβ = m αβ , (cid:80) β q αβ = (cid:80) β m αβ = 1, and λ = 1. This case is treated in Section 6.2 ofNagylaki (1992). It will be convenient to set x α = 1 − r α and y α = 1 − s α , (4.9)where r α ≤ s α ≤ α ∈ { , } . We write the backward migration matrixas M = (cid:18) − m m m − m (cid:19) , (4.10)where 0 < m α < α ∈ { , } .Now we derive the condition for protection of A . The characteristic polynomial of Q is given by ϕ ( x ) = (1 − s )(1 − s ) x − (2 − m − m − s − s + s m + s m ) x + 1 − m − m . (4.11)It is convex and satisfies ϕ (0) = 1 − m − m > , ϕ (1) = s s (1 − κ ) , (4.12a) ϕ (cid:48) (1) = (1 − s )( m − s ) + (1 − s )( m − s ) , (4.12b)where κ = m s + m s . (4.13)By Example 4.1, A is not protected if A A is less fit than A A in both demes (moregenerally, if s ≤ s ≤
0, and s + s < A will be protected if A A is fitter than A A in both demes (more generally, if s ≥ s ≥
0, and s + s > A A is fitter than A A in one deme and less fit in the other, i.e., s s < ϕ ( x ) has two real roots. We have todetermine when the larger ( λ ) satisfies λ >
1. From (4.11) and (4.12), we infer thatthis is the case if and only if ϕ (1) < ϕ (1) > ϕ (cid:48) (1) <
0. It is straightforward to25
01 1 (cid:539) =1 m m s Fig. 1s
01 1 s -1-1 (cid:525) + (cid:525) + (cid:525) (cid:525) (cid:525) (cid:525) (cid:525) (cid:525) (cid:539) =1 (cid:539) =1 (cid:539) = -1 (cid:539) = - Figure 4.1:
The region of protection of A (hatched). From Nagylaki and Lou (2008). show that ϕ (1) > ϕ (cid:48) (1) < s s <
0. Therefore, we concludethat, if s s <
0, allele A is protected if κ < κ >
1. Figure 4.1 displays the region of protectionof A for given m and m .If there is no dominance ( r α = − s α and 0 < | s α | < α = 1, 2), then furthersimplification can be achieved. From the preceding paragraph the results depicted inFigure 4.2 are obtained. The region of a protected polymorphism isΩ + = { ( s , s ) : s s < | κ | < } . (4.15)In a panmictic population, a stable polymorphism can not occur in the absence ofoverdominance. Protection of both alleles in a subdivided population requires that selec-tion 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 selectionacts in opposite direction and dominance is intermediate, i.e., r α s α < α = 1 , s s < . (4.16)26
01 1 (cid:539) =1 m m s Fig. 1s
01 1 s -1-1 (cid:525) + (cid:525) + (cid:525) (cid:525) (cid:525) (cid:525) (cid:525) (cid:525) (cid:539) =1 (cid:539) =1 (cid:539) = -1 (cid:539) = - Figure 4.2:
The regions of protection of A only (Ω ), A only (Ω ), and both A and A (Ω + )in the absence of dominance. From Nagylaki and Lou (2008). Example 4.2.
It is illuminating to study how the parameter region in which a protectedpolymorphism exists depends on the degree of dominance in the two demes. Figures4.3 and 4.4 display the regions of protected polymorphism for two qualitatively differentscenarios of dominance. In the first, the fitnesses are given by A A A A A A s hs − s − as − has as , (4.17)thus, there is deme independent degree of dominance (DIDID). In the second scenario,the fitnesses are given by A A A A A A s hs − s − as has as . (4.18)In both cases, we assume 0 < s <
1, 0 < a < a is a measure of the selection intensityin deme 2 relative to deme 1), and − ≤ h ≤ a > / (1 + 2 s ), there exists aprotected polymorphism for every h ≤ m ≤
1. Otherwise, for given a , thecritical migration rate m admitting a protected polymorphism decreases with increasing h .0 0.2 0.4 0.6 0.8 1.0 a0.10.20.30.40.5m Figure 4.3:
Influence of the degree of dominance on the region of protected polymorphism(shaded). The fitness scheme is (4.17) with s = 0 .
1, and migration is symmetric, i.e., m = m = m . The values of h are -0.95, -0.5, 0, 0.5, 0.95 and correspond to the curves from left toright (light shading to dark shading). A protected polymorphism is maintained in the shadedarea to the right of the respective curve. To the right of the white vertical line, a protectedpolymorphism exists for every m . Figure 4.4:
Influence of the degree of dominance on the region of protected polymorphism(shaded). The fitness scheme is (4.18) with s = 0 .
1, and migration is symmetric, i.e., m = m = m . The values of h are -0.75, -0.25, 0, 0.25, 0.75 and correspond to the curves from rightto left (shading from dark to light!). A protected polymorphism is maintained in the shadedarea to the right of the respective curve. A . This can be shown analyticallyby studying the principal eigenvalue.For (4.18), increasing h greatly facilitates protected polymorphism because it leads toan increase of the average fitness of heterozygotes relative to the homozygotes. The preciseargument is as follows. If we rescale fitnesses in (4.18) according to (4.1), the matrix Q in(4.5) has the entries q α = m α (1 + hs ) / (1 − s ) and q α = m α (1 + has ) / (1 − as ), which areincreasing in h . Therefore, the principal eigenvalues λ of Q increases in h (e.g., Bermanand Plemmons 1994, Chapter 1.3), and (4.6) implies that protection of A is facilitatedby increasing h . Because an analogous reasoning applies to A , the result is proved.Indeed, the above finding on the role of dominance for the fitness scheme (4.18) is aspecial case the following result of Nagylaki (personal communication). Assume fitnessesof A A , A A , and A A in deme α are 1 + s α , 1 + h α s α , and 1, respectively, where s α ≥ ≤ h α ≤
1. If the homozygote fitnesses are fixed, increasing the heterozygotefitness in each deme aids the existence of a protected polymorphism. The proof followsby essentially the same argument as above.
Example 4.3.
In the
Deakin model , the condition (4.14) for protection of allele A becomes κ = µ (cid:18) c s + c s (cid:19) < , (4.19)where s s <
0. Therefore, for given s , s , and c , there is a critical value µ such thatallele A is protected if and only if µ < µ . This implies that for two diallelic demes aprotected polymorphism is favored by a smaller migration rate. Example 4.4.
In the
Levene model , the condition for a protected polymorphism is c s + c s < c r + c r < . (4.20)We close this subsection with an example showing that already with two alleles andtwo demes the equilibrium structure can be quite complicated. Example 4.5.
In the absence of migration, the recurrence equations for the allele fre-quencies p , p in the two demes are two decoupled one-locus selection dynamics of theform (2.16). Therefore, if there is underdominance in each deme, the top convergencepattern in Figure 2.1 applies to each deme. As a consequence, in the absence of migra-tion, the complete two-deme system has nine equilibria, four of which are asymptoticallystable and the others are unstable. Under sufficiently weak migration all nine equilibria29re 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 theseequilibria are extinguished in a sequence of bifurcations (Karlin and McGregor 1972a). 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 . Theorem 4.6 (Christiansen 1974) . If there exists a permutation of demes such that det( I ( n ) − Q ( n ) ) < for some n , where ≤ n ≤ Γ , then A is protected. This theorem is sharp in the sense that if the inequality in (4.21) is reversed for every n ≤ Γ, then A is not protected.The simplest condition for protection is obtained from Theorem 4.6 by setting n = 1.Hence, A is protected if an α exists such that (Deakin 1972) m αα /y α > . (4.22)This condition ensures that, when rare, the allele A increases in frequency in deme α even if this subpopulation is the only one containing the allele. The general conditionin the above theorem ensures that A increases in the n subpopulations if they are theonly ones that contain it. Therefore, we get the following simple sufficient condition fora protected polymorphism:There exist α and β such that m αα /y α > m ββ /x β >
1. (4.23)If we apply these results to the Deakin model with an arbitrary number of demes,condition (4.23) becomesThere exist α and β such that 1 − µ (1 − c α ) y α > − µ (1 − c β ) x β >
1. (4.24)A more elaborate and less stringent condition follows from Theorem 4.6 by nice matrixalgebra: 30 orollary 4.7 (Christiansen 1974) . For the Deakin model with Γ ≥ demes, the followingis a sufficient condition for protection of A . There exists a deme α such that − y α ≥ µ (4.25) or, if (4.25) is violated in every deme α , then µ (cid:88) α c α µ + y α − > . (4.26) If (4.25) is violated for every α and the inequality in (4.26) is reversed, then A is notprotected. Corollary 4.7 can be extended to the inhomogeneous Deakin model, which allows fordifferent homing rates (Christiansen 1974; Karlin 1982, p. 182). Using Corollary 4.7, wecan generalize the finding from two diallelic demes that a lower degree of outbreeding isfavorable for protection of one or both alleles. More precisely, we show:
Corollary 4.8. If < µ ≤ µ ≤ , then allele A is protected for µ if µ satisfies theconditions in Corollary 4.7.Proof. If condition (4.25) holds for µ , it clearly holds for µ . Now suppose that 1 − y α < µ for every α (hence 1 − y α < µ ) and that µ satisfies (4.26). This is equivalent to (cid:88) α c α (1 − y α ) µ + y α − µ + y α − µ + y α − (cid:88) α c α (1 − y α ) µ + y α − > . (4.27)Because ( µ + y α − / ( µ + y α − ≥ µ /µ if and only if y α >
1, we obtain µ µ (cid:88) α c α (1 − y α ) µ + y α − ≥ (cid:88) α c α (1 − y α ) µ + y α − > , (4.28)which proves our assertion.In general, 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 polymor-phism may be destroyed by decreasing one µ α (Karlin 1982, p. 128). Example 4.9.
We apply Corollary 4.7 to the Levene model. Because (4.25) can neverbe satisfied if µ = 1, A is protected from loss if the harmonic mean of the y α is less thanone, i.e., if (Levene 1953) y ∗ = (cid:32)(cid:88) α c α y α (cid:33) − < . (4.29a)31nalogously, allele A is protected if x ∗ = (cid:32)(cid:88) α c α x α (cid:33) − < . (4.29b)Jointly, (4.29a) and (4.29b) provide a sufficient condition for a protected polymorphism.If y ∗ > x ∗ >
1, then A or A , respectively, is lost if initially rare.If A is recessive everywhere ( y α = 1 for every α ), then it is protected if (Prout 1968)¯ x = (cid:88) α c α x α > . (4.30)Therefore, a sufficient condition for a protected polymorphism is x ∗ < < ¯ x . (4.31)Whereas in the Deakin model and in its special case, the Levene model, dispersaldoes not depend on the geographic distribution of niches, in the stepping-stone model itoccurs between neighboring niches. How does this affect the maintenance of a protectedpolymorphism? Here is the answer: Example 4.10.
For the general stepping-stone model (3.22) with fitnesses given by (4.1),Karlin (1982) provided the following explicit criterion for protection of A : (cid:88) α π α y α (cid:30)(cid:88) α π α > , (4.32)where π = 1 and π α = r α − r α − · · · r q α q α − · · · q (4.33)if α ≥
2. This result is a simple consequence of Lemma 4.13.For the homogeneous stepping-stone model with equal demes sizes, i.e., M given by(3.21), (4.32) simplifies to 1Γ (cid:88) α y α > , (4.34)which is the same as condition (4.29) in the Levene model with equal deme sizes.Thus, we obtain the surprising result that the conditions for protection are the samein the Levene model and in the homogeneous stepping stone model provided all demeshave equal size. 32 .4 The continent-island model We consider a population living on an island. At first, we admit an arbitrary number ofalleles. Each generation, a proportion a of adults is removed by mortality or emigrationand a fraction b of migrants with constant allele frequencies ˆ q i is added. A simple in-terpretation is that of one-way migration from a continent to an island. The continentalpopulation is in equilibrium with allele frequencies ˆ q i >
0. Sometimes, ˆ q i is interpreted asthe average frequency over (infinitely) many islands and the model is simply called islandmodel.Assuming random mating on the island, the dynamics of allele frequencies p i on theisland becomes p (cid:48) i = (1 − m ) p i w i ¯ w + m ˆ q i , (4.35)where w i is the fitness of allele A i on the island and m = b/ (1 − a + b ). Thus, m is thefraction of zygotes with immigrant parents. Throughout we assume 0 < m < u ij = m ˆ q j and consider u ij as the mutation rate from A i to A j , a specialcase of the (multiallelic) mutation-selection model is obtained (the so-called house-of-cards model). Therefore (B¨urger 2000, pp. 102-103), (4.35) has the Lyapunov function V ( p ) = ¯ w ( p ) − m (cid:89) i p m ˆ q i i . (4.36)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 inthe population despite immigration of other alleles from the continent. Clearly, no allelecarried to the island recurrently by immigrants can be lost.We investigate the diallelic case and consider alleles A and A with frequencies p and1 − p on the island, and ˆ q = 1 on the continent. Thus, all immigrants are A A . For thefitnesses on the island, we assume w = 1 + s , w = 1 + hs , w = 1 − s , (4.37)where 0 < s ≤ − ≤ h ≤
1. Therefore, A evolves according to p (cid:48) = f ( p ) = (1 − m ) p w ¯ w . (4.38)We outline the analysis of (4.38). Since f (1) = 1 − m <
1, this confirms that A cannot be lost ( p = 1 is not an equilibrium). Because f ( p ) = p (1 − m ) w w + O ( p ) (4.39)33s p →
0, the allele A is protected if m < − w w . (4.40)Obviously, p = 0 is an equilibrium. If there is no other equilibrium, then it must beglobally asymptotically stable (as is also obvious from f (1) <
1, which implies p (cid:48) < p ).The equilibria with p (cid:54) = 0 satisfy ¯ w = (1 − m ) w , (4.41)which is quadratic in p . With the fitnesses (4.37), the solutions areˆ p ± = 14 h (cid:104) h + m (1 − h ) ± (cid:112) (1 − h ) (1 + m ) + 8 hm (1 + 1 /s ) (cid:105) , (4.42)These solutions give rise to feasible equilibria if 0 ≤ p ± ≤
1. As h →
0, (4.42) gives thecorrect limit, ˆ p − = 1 − m/s m if h = 0 . (4.43)We define h = − m − m , (4.44a) µ = 1 + h (1 − m ) , (4.44b) µ = − m − (1 − h ) (1 + m ) h . (4.44c)Then µ > µ if and only if h < h . If µ < µ , then ˆ p + is not admissible. As m increasesfrom 0 to 1, h decreases from − to −
1. Hence, if there is no dominance or the fitterallele A is (partially) dominant (0 ≤ h ≤ h > h . If A is recessive, then h < h .For fixed h and s it is straightforward, but tedious, to study the dependence of ˆ p + andˆ p − on m . The results of this analysis can be summarized as follows (for an illustration,see Figure 4.5): Theorem 4.11. (a) Let m/s < µ . Then 0 and p − are the only equilibria and p ( t ) → ˆ p − as t → ∞ . Thus, for sufficiently weak migration, a unique polymorphism is establishedindependently of the degree of dominance and the initial condition.(b) Let m/s > µ and − ≤ h < h , or m/s ≥ µ and h ≤ h ≤ , i.e., migrationis strong relative to selection. Then there exists no internal equilibrium and p ( t ) → as t → , i.e., the continental allele A becomes fixed. ) Scaled migration rate, m/s A ll e l e f r equen cy , p m/s A ll e l e f r equen cy , p Figure 4.5:
Bifurcation patterns for the one-locus continent-island model. Bold lines indicatean asymptotically stable equilibrium, dashed lines an unstable equilibrium. Figure a displaysthe case h ≤ h ≤
1, and Figure b the case − ≤ h < h . The parameters are s = 0 .
1, and h = 0 . h = − . µ ≈ .
41; in b, µ ≈ . µ ≈ . (c) Let µ < m/s ≤ µ and − ≤ h < h , i.e., migration is moderately strong relativeto selection and the fitter allele is (almost) recessive. Then there are the three equilibria0, ˆ p + , and ˆ p − , where < ˆ p + < ˆ p − , and p ( t ) → if p (0) < ˆ p + , (4.45a) p ( t ) → ˆ p − if p (0) > ˆ p + . (4.45b)For a detailed proof, see Nagylaki (1992, Chapter 6.1). The global dynamics of thediallelic continent-island model admitting evolution on the continent is derived in Nagylaki(2009a).Sometimes immigration from two continents is considered (Christiansen 1999). Someauthors call our general migration model with n demes the n -island model (e.g., Chris-tiansen 1999). The following (symmetric) island model is a standard model in investiga-tions of the consequences of random drift and mutation in finite populations: m αα = 1 − m , (4.46a) m αβ = m Γ − α (cid:54) = β . (4.46b)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.46) reduces to a specialcase of the Deakin model, (3.19), with m = µ (1 − Γ − ). In this, and only this case, theisland model is conservative and the stationary distribution of M is c = e/ Γ.35 .5 Submultiplicative fitnesses
Here, we admit arbitrary migration patterns but assume that fitnesses are submultiplica-tive, i.e., the fitnesses in (4.1) satisfy x α y α ≤ α (Karlin and Campbell 1980). We recall from Section 2.2 that with multiplica-tive fitnesses ( x α y α = 1), the diploid selection dynamics reduces to the haploid dynamics.Throughout, we denote left principal eigenvector of M by ξ ∈ S Γ ; cf. (7.24). Theorem 4.12 (Karlin and Campbell 1980, Result I) . If Γ ≥ and fitnesses are submul-tiplicative in both demes, then at most one monomorphic equilibrium can be asymptoticallystable and have a geometric rate of convergence. The proof of this theorem applies the following very useful result to the conditions forprotection of A or A (Section 4.1). Lemma 4.13 (Friedland and Karlin 1975) . If M is a stochastic n × n matrix and D isa diagonal matrix with entries d i > along the diagonal, then ρ ( M D ) ≥ n (cid:89) i =1 d ξ i i (4.48) holds, where ρ ( M D ) is the spectral radius of M D . Without restrictions on the migration matrix, the sufficient condition (cid:89) α (1 /y α ) ξ α > A is sharp. This can be verified for a circulant stepping stone model,when ξ α = 1 / Γ.The conclusion of Theorem 4.12 remains valid under some other conditions (see Karlinand Campbell 1980). For instance, if migration is doubly stochastic, then ξ α = 1 / Γ forevery α and ρ ( M D ) ≥ (cid:81) α d / Γ α . Therefore, the compound condition (cid:89) α x α y α < heorem 4.14 (Karlin and Campbell 1980, Result III) . If Γ = 2 and fitnesses aremultiplicative, then there exists a unique asymptotically stable equilibrium which is globallyattracting.
Because the proof given by Karlin and Campbell (1980) is erroneous, we present acorrected proof (utilizing their main idea).
Proof of Theorem 4.14.
Because fitnesses are multiplicative, we can treat the haploidmodel. Without loss of generality we assume that the fitnesses of alleles A and A in deme 1 are s ≥ s (0 < s <
1) and 1.(We exclude the case when one allele is favored in both demes, hence goes to fixation.)Let the frequency of A in demes 1 and 2 be p and p , respectively. Then the recursionbecomes p (cid:48) = (1 − m ) s p s p + 1 − p + m s p s p + 1 − p , (4.51a) p (cid:48) = m s p s p + 1 − p + (1 − m ) s p s p + 1 − p . (4.51b)Solving p (cid:48) = p for p , we obtain p = p [(1 − p )(1 − s ) + s m ] p (1 − p )(1 − s )(1 − s ) + m ( s + s p − s p ) . (4.52)Substituting this into p (cid:48) = p , we find that every equilibrium must satisfy p (1 − p )(1 − p + s p ) A ( p ) = 0 , (4.53)where A ( p ) = m [(1 − s )(1 − s + m s ) + m s (1 − s )] − (1 − s ) { (1 − m )(1 − s )(1 − s ) + m [1 − s (1 − m ) + s (1 − m − s )] } p + (1 − m )(1 − s ) (1 − s ) p . (4.54a)We want to show that A ( p ) has always two zeros, because then zeros in (0 ,
1) can emergeonly by bifurcations through either p = 0 or p = 1. Therefore, we calculate the discrim-inant D = (1 − m ) σ τ + 2 m (1 − m ) στ [ m (2 + σ − τ ) − στ ]+ m [ m ( σ + τ ) − m σ (2 + σ − τ ) τ + σ τ ] , (4.55a)37here we set s = 1 + σ and s = 1 − τ with σ ≥ < τ < D as a (quadratic) function in m . To show that D ( m ) > < m <
1, we compute D (0) = (1 − m ) σ τ > , D (1) = m [ τ − σ (1 − τ )] > , (4.56)and D (cid:48) (0) = 2(1 − m ) στ [ m (2 + σ − τ ) − στ ] . (4.57)We distinguish two cases.1. If m (2 + σ − τ ) ≥ στ , then D (cid:48) (0) ≥ , D is concave, then D ( m ) > ,
1] because D (1) >
0; if D is convex, then D (cid:48) ( m ) > m >
0, hence D ( m ) > , m (2 + σ − τ ) < στ , whence D (cid:48) (0) < m < στ σ − τ , we obtain D (cid:48) (1) = 2 m [ − στ (2 + σ − τ − στ ) + m ( σ (1 − τ ) + τ + στ )] < − m σ (1 + σ )(1 − τ ) τ σ − τ < . (4.58a)Therefore, D has no zero in (0 , A ( p ) = 0 has always tworeal solutions.For the rest of the proof we can follow Karlin and Campbell: If s = 1 and s < p = 0 is stable and p = 1 is unstable. As s increases from 1, a bifurcation event at p = 0 occurs before p = 1 becomes stable.At this bifurcation event, p = 0 becomes unstable and a stable polymorphic equilibriumbifurcates off p = 0. Because the constant term m [(1 − s )(1 − s + m s ) + m s (1 − s )]in A ( p ) is linear in s , no further polymorphic equilibria can appear as s increases until p = 1 becomes stable, which then remains stable for larger s .Karlin and Campbell (1980) pointed out that the dynamics (4.51) has the followingmonotonicity property: If p < q and p < q , then p (cid:48) < q (cid:48) and p (cid:48) < q (cid:48) . By considering1 − p and 1 − q , p < q and p > q implies p (cid:48) < q (cid:48) and p (cid:48) > q (cid:48) . In particular, if p < p (cid:48) and p < p (cid:48) , then p (cid:48) < p (cid:48)(cid:48) and p (cid:48) < p (cid:48)(cid:48) . Therefore, this monotonicity property impliesglobal monotone convergence to the asymptotically stable equilibrium. Remark 4.15.
Because the monotonicity property used in the above proof holds for anarbitrary number of demes and for selection on diploids, Karlin and Campbell (1980)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 demes38nd Akin (personal communication) has shown that for three diallelic demes unstableperiodic orbits can occur in the corresponding continuous-time model (cf. Section 6.2).For two diallelic demes, however, every trajectory converges to an equilibrium.Karlin and Campbell (1980) showed that for some migration patterns, Theorem 4.14remains true for an arbitrary number of demes. The Levene model is one such example. Infact, Theorem 5.16 is a much stronger result. It is an open problem whether Theorem 4.14holds for arbitrary migration patterns if there are more than two demes. An interestingand related result is the following.
Theorem 4.16 (Karlin and Campbell 1980, Result IV) . Let Γ ≥ and let M be arbi-trary but fixed. If there exists a unique, globally attracting equilibrium under multiplicativefitnesses, then there exists a unique, globally attracting equilibrium for arbitrary submul-tiplicative fitnesses. This is a particularly simple model to examine the consequences of spatially varyingselection for the maintenance of genetic variability. It can also be interpreted as a modelof frequency-dependent selection (e.g., Wiehe and Slatkin 1998). From the definition(3.20) of the Levene model and equation (3.5a), we infer p (cid:48) i,α = (cid:88) β c β p ∗ i,β , (5.1)which is independent of α . Therefore, after one round of migration allele frequencies arethe same in all demes, whence it is sufficient to study the dynamics of the vector of allelefrequencies p = ( p , . . . , p I ) (cid:62) ∈ S I . (5.2)The migration-selection dynamics (3.5) simplifies drastically and yields the recurrenceequation p (cid:48) i = p i (cid:88) α c α w i,α ¯ w α ( i ∈ I ) (5.3)for the evolution of allele frequencies in the Levene model. Therefore, there is no popula-tion structure in the Levene model although the population experiences spatially varyingselection. In particular, distance does not play any role. The dynamics (5.3) is also ob-tained if, instead of the life cycle in Section 3.2, it is assumed that adults form a commonmating pool and zygotes are distributed randomly to the demes.39 emark 5.1.
1. In the limit of weak selection, the Levene model becomes equivalent topanmixia. Indeed, with w ij,α = 1 + sr ij,α , an argument analogous to that below (2.20)shows that the resulting dynamics is of the form˙ p i = p i (cid:88) α c α ( w i,α − ¯ w α ) = p i ( z i − ¯ z ) , (5.4)where z i = (cid:80) j z ij p j = (cid:80) j (cid:80) α c α r ij,α p j is linear in p and ¯ z = (cid:80) i z i p i . Here, z ij = (cid:80) α c α r ij,α is the spatially averaged selection coefficient of A i A j (Nagylaki and Lou 2001).Therefore, the Levene model is of genuine interest only if selection is strong. For arbi-trary 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 (Lou and Nagylaki2008).3. For hard selection (3.10), the dynamics in the Levene model becomes p (cid:48) i = p i w i ¯ w ( i ∈ I ) , where w i = (cid:88) α c α w i,α , (5.5)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 protectedpolymorphism if and only if the averaged fitnesses display heterozygous advantage. Inter-estingly, it is not true that the condition for protection of an allele under hard selection isalways more stringent than under soft selection, although under a range of assumptionsthis is the case (see Nagylaki 1992, Sect. 6.3). We define ˜ w = ˜ w ( p ) = (cid:89) α ¯ w α ( p ) c α , (5.6)which is the geometric mean of the mean fitnesses in single demes. Furthermore, we define F ( p ) = ln ˜ w ( p ) = (cid:88) α c α ln ¯ w α ( p ) . (5.7)Both functions will play a crucial role in our study of the Levene model.From ∂ ¯ w α ∂p i = ∂∂p i (cid:88) i,j p i p j w ij,α = 2 (cid:88) j p j w ij,α = 2 w i,α , (5.8)40e obtain p (cid:48) i = p i (cid:88) α c α ¯ w α ∂ ¯ w α ∂p i = 12 p i (cid:88) α c α ∂ ln ¯ w α ∂p i = 12 p i ∂F ( p ) ∂p i . (5.9)Because (cid:80) i p (cid:48) i = 1, we can write (5.9) in the form p (cid:48) i = p i ∂ ˜ w ( p ) ∂p i (cid:30)(cid:88) j p j ∂ ˜ w ( p ) ∂p j = p i ∂F ( p ) ∂p i (cid:30)(cid:88) j p j ∂F ( p ) ∂p j . (5.10)A simple exercise using Lagrange multipliers shows that the equilibria of (5.10), henceof (5.3), are exactly the critical points of ˜ w , or F . Our aim is to prove the followingtheorem. Theorem 5.2 (Li 1955; Cannings 1971; Nagylaki 1992, Sect. 6.3) . (a) Geometric-meanfitness satisfies ∆ ˜ w ( p ) ≥ for every p ∈ S I , and ∆ ˜ w ( p ) = 0 if and only if p is anequilibrium of (5.3) . The same conclusion holds for F .(b) The set Λ of equilibria is globally attracting, i.e., p ( t ) → Λ as t → ∞ . If everypoint in Λ is isolated, as is generic , then p ( t ) converges to some ˆ p ∈ Λ . Generically, p ( t ) converges to a local maximum of ˜ w . This is an important theoretical result because it states that in the Levene model nocomplex dynamical behavior, such as limit cycles or chaos, can occur. One immediateconsequence of this theorem is
Corollary 5.3.
Assume the Levene model with two alleles. If there exists a protectedpolymorphism, then all trajectories converge to an internal equilibrium point.Outline of the proof of Theorem 5.2. (a) We assume c α ∈ Q for every α . Then thereexists n ∈ N such that W = ˜ w n is a homogeneous polynomial in p with nonnegativecoefficients. From ∂W∂p i = n ˜ w n − ∂ ˜ w∂p i , we infer ∂W∂p i (cid:30)(cid:88) j p j ∂W∂p j = ∂ ˜ w∂p i (cid:30)(cid:88) j p j ∂ ˜ w∂p j . From the inequality of Baum and Eagon (1967), we infer immediately W ( p (cid:48) ) > W ( p )unless p (cid:48) = p . Therefore, W and ˜ w are (strict) Lyapunov functions. By a straightforwardbut tedious approximation argument, which uses compactness of S I , it follows that ˜ w is We call a property generic if it holds for almost all parameter combinations and, if applicable, foralmost all initial conditions c α ∈ R (see Nagylaki 1992, Sect. 6.3). This proves the firstassertion 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(LaSalle 1976, p. 10). The second is obvious, and the third follows because ˜ w is nonde-creasing. (Note, however, that convergence to a maximum holds only generically becausesome trajectories may converge to saddle points. Moreover, because ˜ w may also haveminima and saddle points, the globally attracting set Λ may not be stable; only a subsetof Λ is stable.)Next, we derive a useful criterion for the existence of a unique stable equilibrium. Lemma 5.4. If ¯ w α is concave for every α , then F is concave.Proof. Because the logarithm is strictly monotone increasing and concave, a simple esti-mate shows that ln ¯ w α is concave. Hence, F is concave. Theorem 5.5 (Nagylaki and Lou 2001) . If F is concave, there exists exactly one stableequilibrium (point or manifold) and it is globally attracting. If there exists an internalequilibrium, it is the global attractor.Proof. Concavity and the fact that ∆ F ≥ S I . Because F is concave, F is constant on the line that joins them. If that line is on theboundary of S I , then these two equilibrium points are elements of the same manifold ofequilibria (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 w ij,α = v i,α + v j,α (5.11)for every i, j ∈ I and α ∈ G , whence (3.2) yields w i,α = v i,α + ¯ v α and ¯ w α = 2¯ v α , (5.12)where ¯ v α = (cid:88) i v i,α p i . (5.13)42 simple calculation shows that without dominance the dynamics (5.3) simplifies to p (cid:48) i = p i (cid:18) (cid:88) α c α v i,α ¯ v α (cid:19) ( i ∈ I ) . (5.14)Fitnesses are called multiplicative if there exist constants v i,α such that w ij,α = v i,α v j,α (5.15)for every i, j ∈ I and α ∈ G . Then (3.2) yields w i,α = v i,α ¯ v α and ¯ w α = ¯ v α . (5.16)With multiplicative fitnesses, one obtains the haploid Levene model, i.e., p (cid:48) i = p i (cid:88) α c α v i,α / ¯ v α ( i ∈ I ) . (5.17) Theorem 5.6 (Nagylaki and Lou 2001) . The function F is concave in the following cases:(a) In every deme there is no dominance.(b) In every deme fitnesses are multiplicative.(c) In every deme a globally attracting internal equilibrium exists without migration(i.e., with the dynamics p (cid:48) i,α = p i,α w i,α / ¯ w α ).Therefore, the conclusions of Theorem 5.5 apply in each of the cases.Proof. (a) Without dominance, ¯ w α = 2¯ v α is linear in every deme.(b) With multiplicative fitnesses, ln ¯ w α = 2 ln ¯ v α , which is concave because ¯ v α is linear.(c) If there is a globally attracting internal equilibrium in deme α when it is isolated,then the fact that ∆ ¯ w α ≥ w α 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 by ¯ w . Forhaploids, Strobeck (1979) proved uniqueness and asymptotic stability of an internal equi-librium.Nagylaki and Lou (2006b) provide sufficient conditions for the nonexistence of aninternal equilibrium.The following theorem shows that a globally stable internal equilibrium can be main-tained on an open set of parameters, even in the absence of overdominance. Clearly,43his requires that there are at least two demes; cf. Section 2.3. We say there is partialdominance in every deme if w ii,α < w ij,α < w jj,α or w jj,α < w ij,α < w ii,α (5.18)holds for every pair i, j ∈ I with i (cid:54) = j and for every α ∈ Γ. Theorem 5.8.
Assume an arbitrary number of alleles and Γ ≥ demes. Then thereexists a nonempty open set of fitness parameters exhibiting partial dominance such thatfor every parameter combination in this set, there is a unique, internal, asymptoticallystable equilibrium point of the dynamics (5.3) . This equilibrium is globally asymptoticallystable. This theorem arises for the special case of a single locus in Result 5.1 in B¨urger (2010).The latter is an immediate consequence of Theorem 3.1 in B¨urger (2010) and of Theorem2.2, Remark 2.3(iii), and Remark 2.4 in B¨urger (2009b); cf. Theorem 7.4 below. Theproof shows that the maintenance of a stable internal equilibrium requires a certain formof spatially averaged overdominance that, with intermediate dominance, can be achievedonly if the direction of selection and the degree of dominance vary among demes.
Throughout this subsection, we assume no dominance, i.e., (5.11). For this importantspecial case more detailed results can be proved than for general (intermediate) dom-inance. Nevertheless, internal equilibrium solutions can be determined explicitly onlyunder specific assumptions (e.g., Theorem 4.2 in Nagylaki and Lou 2006b). An interest-ing question to ask is: What determines the number of alleles that can be maintained atan equilibrium? If there is no dominance, there is a simple answer.
Theorem 5.9 (Nagylaki and Lou 2001) . Without dominance, the number of demes is ageneric upper bound on the number of alleles present at equilibrium. Any neutral demeshould not be counted in this bound.Proof.
From (5.14) we conclude that, at an internal equilibrium ˆ p , (cid:88) α c α v i,α ˆ¯ v α = 1 ( i ∈ I ) (5.19)holds. The substitution x α = c α / ˆ¯ v α , where ˆ¯ v α = (cid:88) i v i,α ˆ p i , (5.20)44inearizes (5.19): (cid:88) α v i,α x α = 1 ( i ∈ I ) . (5.21)This is a system of I inhomogeneous linear equations for the Γ unknowns x α . Therefore,a solution exists generically only if I ≤ Γ. (Note that even if there exists a solution ( x α ),it does not necessarily give rise to a solution ˆ p ∈ S I .)The statement about neutral demes follows because in a neutral deme v i,α = ¯ v α = 1for every i , hence x α = c α holds. Remark 5.10.
If in Theorem 5.9 general intermediate dominance is admitted, then thereexists an open set of parameters for which any number of alleles can be maintained at anasymptotically stable equilibrium. Theorem 7.2 provides a much more general result.
Remark 5.11.
Theorem 5.9 holds for hard selection. A slight modification of the proofshows that it holds also for multiplicative fitnesses. In fact, Strobeck (1979) proved ananalog of Theorem 5.9 for haploid species.
Example 5.12 (Nagylaki and Lou 2001) . The upper bound established in Theorem 5.9can be assumed if Γ = I . Let v i,α = u i δ iα , where u i > δ iα is the Kronecker delta.This means that allele A i has fitness u i in deme i and fitness 0 elsewhere. Hence, everyallele is the best in one deme. Then (5.14) simplifies to p (cid:48) i = p i (cid:18) u i ¯ v i (cid:19) = p i (1 + c i /p i ) = ( p i + c i ) . (5.22)The solution is p i ( t ) = c i + ( p i (0) − c i )( ) t → c i as t → ∞ .In the formulation of Theorem 5.9, the word ‘generic’ is essential. If Γ and I arearbitrary 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 globallyattracting) manifold of equilibria exists for an open set of parameter combinations. Thisholds also for the additive case, when r ij = s i + s j (Nagylaki and Lou 2001).Another interesting problem is to determine conditions when a specific allele will goto fixation. A simple and intuitive result is the following: Theorem 5.13 (Nagylaki and Lou 2006b) . Suppose there exists i ∈ I such that (cid:88) α c α v j,α v i,α < for every j (cid:54) = i . Then p i ( t ) → as t → ∞ . Phase portrait for the Levene model in Example 5.14 (Figure from Nagylaki andLou 2001).
A proof as well as additional results on the loss or fixation of alleles may be foundin Nagylaki and Lou (2006b). The following example is based on a nice application ofTheorem 5.6 and shows that migration may eliminate genetic variability.
Example 5.14 (Nagylaki and Lou 2001) . We suppose two demes of equal size ( c = c = ) and three alleles without dominance. The alleles A and A have extreme fitnessesand A is intermediate in both demes. More precisely, we assume v , = 1, v , = 0, v , = u , and v , = 0, v , = 1, v , = u , where 0 < u <
1. Without migration, A isultimately fixed in deme 1 and A is fixed in deme 2. We now establish that the Levenemodel evolves as sketched in Figure 5.1, i.e., p ( t ) → (cid:40) ( , ,
0) if 0 < u < , (0 , ,
1) if ≤ u < . (5.25)Because there is no dominance, Theorem 5.6 informs us that there can be only one stableequilibrium. Therefore, it is sufficient to establish asymptotic stability of the limitingequilibrium, which we leave as an easy exercise.Following Nagylaki (2009a), we say there is deme-independent degree of intermediatedominance (DIDID) if w ij,α = ϑ ij w ii,α + ϑ ji w jj,α (5.26a)holds for constants ϑ ij such that0 ≤ ϑ ij ≤ ϑ ji = 1 − ϑ ij (5.26b)for every α and every pair i, j . In particular, ϑ ii = .46bviously, DIDID covers complete dominance or recessiveness ( ϑ ij = 0 or ϑ ij = 1 if i (cid:54) = j ), and no dominance ( ϑ ij = ), but not multiplicativity. We also note that DIDIDincludes the biologically important case of absence of genotype-by-environment interac-tion. In general, F is not concave under DIDID. However, Nagylaki (2009a, Theorem 3.2)proved that with DIDID the evolution of p ( t ) is qualitatively the same as without dom-inance. Therefore, the conclusion of Theorem 5.5 holds, i.e., under DIDID there existsexactly one stable equilibrium (point or manifold) and it is globally attracting. If thereexists an internal equilibrium, it is the global attractor. Moreover, the number of demesis 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, thecondition for loss of an allele in Theorem 5.13 has a simple generalization to DIDID. Forproofs and further results about DIDID consult Nagylaki (2009a). Remark 5.15.
For general migration, the dynamics under DIDID may differ qualitativelyfrom that under no dominance. For instance, in a diallelic model with one way migration,two asymptotically stable equilibria may coexist if there is DIDID (Nagylaki 2009a).Moreover, in two triallelic demes, a globally asymptotically stable internal equilibriumexists for an open set of parameters (Peischl 2010). Therefore, with arbitrary migrationand DIDID, the number of alleles maintained at a stable equilibrium can be higher thanthe number of demes. This is not the case in the absence of dominance when the numberof demes is a generic upper bound (Theorem 6.1).
Here, we specialize to two alleles but admit arbitrary dominance. We assume that fitnessesare given by (4.1). Our first result yields a further class of examples when a unique stableequilibrium exists.
Theorem 5.16 (B¨urger 2009c) . For every α , let the fitnesses satisfy x α y α ≤ − y α ) and y α ≤ or x α y α ≤ − x α ) and x α ≤ . (5.27b) Then F is concave on [0 , . Hence, there exists at most one internal equilibrium. If an in-ternal equilibrium exists, it is globally asymptotically stable. If a monomorphic equilibriumis stable, then it is globally asymptotically stable.
47n the proof it is shown that ln ¯ w α is concave if and only if (5.27) is fulfilled. Then theconclusion follows from Theorem 5.5.Theorem 5.16 shows that the protection conditions (4.29) imply the existence of aglobally asymptotically stable internal equilibrium if (5.27) holds. It generalizes Result Iin Karlin (1977), who proved uniqueness of an internal equilibrium and global convergenceunder the assumption of submultiplicative fitnesses (4.47). His proof is based on a differentmethod. Submultiplicative fitnesses, hence fitnesses satisfying (5.27), include a numberof important cases: multiplicative fitnesses (hence, selection on haploids), no dominance,partial or complete dominance of the fitter allele, and overdominance. Example 5.17.
If there are two niches, no dominance, and fitnesses are given by 1 + s α ,1, and 1 − s α ( s s < | κ | < , where κ = c s + c s . (5.28)By Theorem 5.6 or 5.16 there is a unique, asymptotically stable equilibrium ˆ p . A simplecalculation yields ˆ p = (1 − κ ) / Example 5.18.
If there are two niches and multiplicative fitnesses given by 1 / (1 − s α ),1, and 1 − s α ( s s < < − κ < , (5.29)where κ is as above. The unique, asymptotically stable equilibrium is given by ˆ p = 1 − κ .A partial converse to Theorem 5.16 is the following: Theorem 5.19 (Karlin 1977, Result IA) . Suppose x α y α > { ( x α − , ( y α − } for every α , (5.30) and there exists at least one internal equilibrium. Then both monomorphic equilibria areasymptotically stable and the internal equilibrium is unique. Figure 5.2 displays the regions in which submultiplicativity, (5.27), or (5.30) hold.Table 5.1 presents numerical results which demonstrate that stronger selection, submul-tiplicativity, and a higher number of demes all facilitate protected polymorphism in theLevene model. These findings agree with intuition. For instance, if an allele has suffi-ciently low fitness in just one deme, i.e., c α /y α >
1, the other allele is protected.48igure 5.2:
In the light blue regions fitnesses are submultiplicative, in the dark blue and lightblue regions they satisfy (5.27). In the blue region(s), there exists a single stable equilibrium. Inthe purple region, fitnesses satisfy (5.30). There, both boundary equilibria are stable providedan internal equilibrium exists. In the light blue square there is overdominance, in the lightpurple square there is underdominance.
So far, we derived sufficient conditions for a unique (internal) equilibrium, but we havenot yet considered the question of how many (stable) internal equilibria can coexist. Thisturns 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) canbe written as ∆ p = p (1 − p ) (cid:88) α c α − y α + p ( x α + y α − x α p + 2 p (1 − p ) + y α (1 − p ) . (5.31)Expressing this with a common denominator, we see that the internal equilibria are thesolutions of a polynomial in p of degree 2Γ −
1. Thus, in principle, there can be up to2Γ − Theorem 5.20.
The diallelic Levene model allows for any number j ∈ { , , . . . , − } of hyperbolic internal equilibria with any feasible stability configuration. His numerical results, which admit arbitrary dominance, show that the proportion of49 weak selection moderate selection strong selectionintdom submult intdom submult intdom submult2 0.094 0.185 0.115 0.271 0.220 0.3323 0.124 0.262 0.162 0.416 0.357 0.5164 0.138 0.310 0.190 0.512 0.460 0.63610 0.158 0.468 0.270 0.788 0.787 0.917
Table 5.1:
Proportion of protected polymorphism in the Levene model with intermediate domi-nance. Fitnesses of A A , A A , A A are parameterized as 1+ s α , 1+ h α s α , 1 − s α . The data forweak selection are generated for the weak-selection limit (Remark 5.1). Then there is a protectedpolymorphism if and only if there is average overdominance, i.e., if (cid:80) α c α s α h α > | (cid:80) α c α s α | holds. For moderate or strong selection, the conditions (4.29) were evaluated. In all cases,10 parameter combinations satisfying s α ∈ ( − s, s ), h α ∈ ( − , c α ∈ ( − ,
1) were randomly(uniformly) chosen, and the values c α were normalized. For weak or strong selection, s = 1;for moderate selection, s = 0 .
2. For submultiplicative fitnesses (columns ‘submult’), the 10 parameter combinations satisfy the additional constraint (4.47), which is reformulated as a con-dition for h α . The columns ‘intdom’ contain the data for general intermediate dominance. Theproportion of submultiplicative parameter combinations among all parameter combinations withintermediate dominance is (1 − ln 2) Γ ≈ . Γ . parameter space supporting more than Γ equilibria becomes extremely small if Γ > Because for multiple alleles and arbitrary migration few general results are available, wefocus on three limiting cases that are biologically important and amenable to mathemat-ical analysis. These are weak selection and weak migration, weak migration (relative toselection), 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 tothe strong-migration limit. The latter two are based on a separation of time scales. Withthe help of perturbation theory, results about existence and stability of equilibria, butalso about global convergence, can be derived for an open subset of parameters of the fullmodel. In addition, we report results about the case of no dominance and about uniformselection, i.e., when selection is the same in every deme. We start with no dominance.
The following generalizes Theorem 5.9 for the Levene model. It holds for an arbitrary(constant) backward migration matrix. 50 heorem 6.1 (Nagylaki and Lou 2001, Theorem 2.4) . Without dominance, the numberof demes is a generic upper bound on the number of alleles that can be maintained at anyequilibrium.
This theorem also holds for hard selection. As shown by Peischl (2010), it can notbe extended to DIDID; cf. Remark 5.15. An example showing that the upper bound canbe achieved if I = Γ is obtained by setting v i,α = u i δ iα , where u i >
0. The internalequilibrium can be written as ˆ p = ( I − M ) − M , and convergence is geometric (at arate ≤ ).As we shall see below, with dominance the number of alleles that can be maintainedat equilibrium may depend on the strength and pattern of migration. Following Nagylaki and Lou (2007), we assume that both selection and migration areweak and approximate the discrete migration-selection dynamics (3.5) by a differentialequation which is easier accessible. Accordingly, let w ij,α = 1 + (cid:15)r ij,α and ˜ m αβ = δ αβ + (cid:15) ˜ µ αβ , (6.1)where r ij,α and ˜ µ αβ are fixed for every i, j ∈ I and every α, β ∈ G , and (cid:15) > w i,α = 1 + (cid:15)r i,α and ¯ w α = 1 + (cid:15) ¯ r α , (6.2a)where r i,α = (cid:88) j r ij,α p j,α and ¯ r α = (cid:88) i,j r ij,α p i,α p j,α . (6.2b)To approximate the backward migration matrix M , note that (3.10) and (6.2a) implythat, for both soft and hard selection, c ∗ α = c α + O ( (cid:15) ) (6.3)as (cid:15) →
0. Substituting (6.1) and (6.3) into (3.8) leads to m αβ = δ αβ + (cid:15)µ αβ + O ( (cid:15) ) (6.4)as (cid:15) →
0, where µ αβ = 1 c α (cid:32) c β ˜ µ βα − δ αβ (cid:88) γ c γ ˜ µ γα (cid:33) . (6.5)51ecause ˜ M is stochastic, we obtain for every α ∈ Γ,˜ µ αβ ≥ β (cid:54) = α and (cid:88) β ˜ µ αβ = 0 . (6.6)As a simple consequence, µ αβ shares the same properties.The final step in our derivation is to rescale time as in Sect. 2.4 by setting t = (cid:98) τ /(cid:15) (cid:99) and π i,α ( τ ) = p i,α ( t ). Inserting all this into the difference equations (3.5) and expandingyields π i,α ( τ + (cid:15) ) = π i,α { (cid:15) [ r i,α ( π · ,α ) − ¯ r α ( π · ,α )] } + (cid:15) (cid:88) β µ αβ π i,β + O ( (cid:15) ) (6.7)as (cid:15) →
0, where π · ,α = ( π ,α , . . . , π I,α ) (cid:62) ∈ S I . Rearranging and letting (cid:15) →
0, we arrive atd π i,α d τ = (cid:88) β µ αβ π i,β + π i,α [ r i,α ( π · ,α ) − ¯ r α ( π · ,α )] . (6.8)Absorbing (cid:15) into the migration rates and selection coefficients and returning to p ( t ), weobtain the slow-evolution approximation of (3.5),˙ p i,α = (cid:88) β µ αβ p i,β + p i,α [ r i,α ( p · ,α ) − ¯ r α ( p · ,α )] . (6.9)In contrast to the discrete-time dynamics (3.5), here the migration and selection terms aredecoupled. This is a general feature of many other slow-evolution limits (such as mutationand selection, or selection, recombination and migration). Because of the decoupling ofthe selection and migration terms, the analysis of explicit models is often facilitated.With multiple alleles, there are no general results on the dynamics of (6.9). For twoalleles, we set p α = p ,α and write (6.9) in the form˙ p α = (cid:88) β µ αβ p β + ϕ α ( p α ) . (6.10)Since µ αβ ≥ α (cid:54) = β , the system (6.10) is quasimonotone or cooperative, i.e., ∂ ˙ p α /∂p β ≥ α (cid:54) = β . As a consequence, (6.10) cannot have an exponentially stable limitcycle. However, Akin (personal communication) has proved for three diallelic demes thata Hopf bifurcation can produce unstable limit cycles. This precludes global convergence,though not generic convergence. If Γ = 2, then every trajectory converges (Hirsch 1982;Hadeler and Glas 1983; see also Hofbauer and Sigmund 1998, p. 28), as is the case in thediscrete-time model (Remark 4.15). 52 xample 6.2. Eyland (1971) provided a global analysis of (6.9) for the special case oftwo diallelic demes without dominance. As in Sect. 4.2, we assume that the fitnesses of A A , A A , and A A in deme α are 1 + s α , 1, and 1 − s α , respectively, where s α (cid:54) = 0( α = 1 , µ = µ > µ = µ >
0, and write p α for the frequencyof A in deme α . Then (6.9) becomes˙ p = µ ( p − p ) + s p (1 − p ) , (6.11a)˙ p = µ ( p − p ) + s p (1 − p ) . (6.11b)The equilibria can be calculated explicitly. At equilibrium, p = 0 if and only if p = 0,and p = 1 if and only if p = 1. In addition, there may be an internal equilibrium point.We set σ α = µ α s α , κ = σ + σ , (6.12)and B = (1 − σ σ ) / . (6.13)The internal equilibrium exists if and only if s s < | κ | <
1; cf. (4.15). If s < < s ,it is given by ˆ p = (1 + B ) − σ and ˆ p = (1 − B ) − σ . (6.14)It is straightforward to determine the local stability properties of the three possibleequilibria. Gobal asymptotic stability follows from the results cited above about quasi-monotone systems. Let p = ( p , p ) (cid:62) . Then allele A is eliminated in the region Ω inFigure 4.2, i.e., p ( t ) → (0 ,
0) as t → ∞ , whereas A is ultimately fixed in the region Ω .In Ω + , p ( t ) converges globally to the internal equilibrium point ˆ p 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 equilibriawere derived by Karlin and Campbell (1980). Some of them are treated in 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, generalconditions were derived by Nagylaki and Lou (2007). Suppose that there exists i ∈ I andconstants γ ij such that γ ij ≥ , γ ii = 0 , (cid:88) j γ ij = 1 , (6.15a)and (cid:88) j γ ij r jk,α > r ik,α (6.15b)53or every α ∈ G and every k ∈ I . Theorem 6.3 (Nagylaki and Lou 2007, Theorem 3.5 and Remark 3.7) . If the matrix ( µ αβ ) is irreducible and the conditions (6.15) are satisfied, then p i ( t ) → as t → ∞ whenever p i (0) > and p j (0) > for every j such that γ ij > . 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 (cid:88) j γ ij s j,α > s i,α , (6.16)and Theorem 6.3 applies.To highlight one of the biological implications, we follow Remark 3.12 in Nagylaki andLou (2007) and assume γ i > γ iI >
0, and γ ij = 0 for j = 2 , . . . , I −
1. Then (6.16)becomes γ i s ,α + γ iI s I,α > s i,α (6.17)for every α , and the theorem shows that p i ( t ) → t → ∞ for i = 2 , . . . , I −
1. If, inaddition to (6.17), we assume thatmin( s ,α , s I,α ) < s i,α < max( s ,α , s I,α ) (6.18)for i = 2 , . . . , I − α , then every allele A i with 1 < i < I is intermediate inevery deme, and A and A I are extreme. Thus, Theorem 6.3 implies that all intermediatealleles are eliminated. This conclusion can be interpreted as the elimination of general-ists by specialists, and it can yield the increasing phenotypic differentiation required forparapatric speciation (cf. Lou and Nagylaki 2002 for an analogous result and discussionin the context of diffusion models). If s ,α − s I,α changes sign among demes, so that everyallele is the fittest in some deme(s) and the least fit in the other(s), then both alleles maybe maintained.An other application of Theorem (6.3) is the following. If in every deme the homozy-gotes have the same fitness order and there is strict heterozygote intermediacy, then theallele with the greatest homozygous fitness is ultimately fixed (Remark 3.20 in Nagylakiand Lou 2007).
Remark 6.4.
The slow-evolution approximation of the exact juvenile-migration modelis also (6.9) (Nagylaki and Lou 2008, Nagylaki 1992, pp. 143-144).54 .3 Weak migration
If migration is sufficiently weak relative to selection, properties of the dynamics can beinferred by perturbation techniques from the well-understood case of a finite number ofisolated demes in which there is selection and random mating.To study weak migration, we assume m αβ = δ αβ + (cid:15)µ αβ , (6.19)where µ αβ is fixed for every α, β , and (cid:15) > M is stochastic,we obtain for every α ∈ Γ (cf. Section 6.2) µ αβ ≥ β (cid:54) = α and (cid:88) β µ αβ = 0 . (6.20)If there is no migration ( (cid:15) = 0), the dynamics in each deme reduces to the pure selectiondynamics p (cid:48) i,α = p i,α w i,α ¯ w α . (6.21)Because (6.21) is defined on the Cartesian product S Γ I , it may exhibit a richer equilibriumand stability structure than the panmictic selection dynamics (2.6). This was alreadyillustrated by Example 4.5, in which two asymptotically stable internal equilibria maycoexist. 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 conse-quences. Among others, it excludes complex dynamics in the full system (3.5) providedmigration is sufficiently weak. We note that hyperbolicity is a generic property for (6.21)(Appendix A in Nagylaki et al. 1999), and an internal equilibrium is hyperbolic if andonly if it is isolated (Lemma 3.2 in Nagylaki and Lou 2006a). Theorem 6.5 (Nagylaki and Lou 2007, Theorem 4.1) . Suppose that every equilibrium of (6.21) is hyperbolic, that (6.19) holds, and that (cid:15) > is sufficiently small.(a) The set of equilibria Σ ⊂ S Γ I of (6.21) contains only isolated points, as does theset of equilibria Σ (cid:15) ⊂ S Γ I of (3.5) . As (cid:15) → , each equilibrium in Σ (cid:15) converges to thecorresponding equilibrium in Σ .(b) In the neighborhood of each asymptotically stable equilibrium in Σ , there existsexactly one equilibrium in Σ (cid:15) , and it is asymptotically stable. In the neighborhood of eachunstable internal equilibrium in Σ , there exists exactly one equilibrium in Σ (cid:15) , and it isunstable. In the neighborhood of each unstable boundary equilibrium in Σ , there exists atmost one equilibrium in Σ (cid:15) , and if it exists, it is unstable.(c) Every solution p ( t ) of (3.5) converges to one of the equilibrium points in Σ (cid:15) . (cid:15) , x (cid:48) = f ( x, (cid:15) ) , (6.22)where x ∈ X ⊆ R n ( X compact and convex) and (cid:15) ∈ E ⊆ R k ( E open). We assume that ˆ x is a fixed point of (6.22) if (cid:15) = (cid:15) and that the Jacobian f (cid:48) (ˆ x, (cid:15) ) of f evaluated at (ˆ x, (cid:15) )exists and is continuous. Furthermore, we posit that ˆ x is a hyperbolic equilibrium. If wedefine the function F ( x, (cid:15) ) = f ( x, (cid:15) ) − x , (6.23)then F (ˆ x, (cid:15) ) = 0 and the Jacobian F (cid:48) (ˆ x, (cid:15) ) is continuous and nonsingular (by hyper-bolicity of ˆ x ). Therefore, the implicit function theorem shows that there exists an openneighborhood U of (cid:15) and a function φ : U → X such that φ ( (cid:15) ) = ˆ x and F ( φ ( (cid:15) ) , (cid:15) ) = 0for (cid:15) ∈ U , hence f ( φ ( (cid:15) ) , (cid:15) ) = φ ( (cid:15) ) . (6.24)Hence, for every (cid:15) ∈ U , i.e., for (cid:15) close to (cid:15) , (6.22) has a uniquely determined fixed pointˆ x (cid:15) = φ ( (cid:15) ) close to ˆ x .With the help of the Hartman-Grobman theorem, it can be shown that the stabilityproperties of the perturbed fixed points ˆ x (cid:15) are the same as those of the unperturbed,ˆ x . The reason is that if an equilibrium is hyperbolic, this property persists under smallperturbations (the Jacobian changes continuously if parameters change continuously); seealso Theorem 4.4 in Karlin and McGegor (1972b). The extension of this argument tofinitely many hyperbolic equilibria is evident.Although hyperbolic equilibria change continuously under small perturbations, limitsets of trajectories do not: they can explode. Thus, perturbations could introduce ‘new’limit sets away from the hyperbolic equilibria. What has good properties under pertur-bations is the set of chain-recurrent points introduced by Conley (1978).Let X be a compact set with metric d and let f : X → X be a continuous map. Apoint x ∈ X is called chain recurrent (with respect to f ) if, for every δ >
0, there existsa finite sequence x = x , x , . . . , x r − , x r = x (often called a δ -pseudo-orbit) such that d ( f ( x m ) , x m +1 ) < δ for m = 0 , , . . . , r −
1. The set of chain-recurrent points contains thelimit sets of all orbits and behaves well under perturbations (Akin 1993, p. 244).56 utline of the proof of Theorem 6.5. (a) and (b) follow from the implicit function theo-rem and the Hartman-Grobman theorem. That asymptotically stable boundary equilibriaremain in the state space follows from Brouwer’s fixed point theorem (for details, see Kar-lin and McGregor 1972, 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 in each deme, ∆ ¯ w α ≥ w ( p ) = (cid:88) α ¯ w α ( p · ,α ) (6.26)satisfies ∆ ¯ w ( p ) ≥ w ( p ) = 0 if and only if p ∈ F . Because, on F , ¯ w takes onlyfinitely many values, Theorem 3.14 in Akin (1993) implies that the chain-recurrent pointsof (6.21) are exactly the equilibria.Now we can follow the proof of Theorem 2.3 in Nagylaki et al. (1999) almost verbally.Because the set of chain-recurrent points consists only of hyperbolic equilibria, this isalso true for small C perturbations of the dynamics (Akin 1993, p. 244). Indeed, as animmediate consequence of the definition of chain recurrence, it follows that the chain-recurrent set of (6.21) changes in an upper semicontinuous way with (cid:15) . In particular,the chain-recurrent set for (cid:15) > δ -neighborhoods ofthe equilibria for (cid:15) = 0, with δ → (cid:15) →
0. By the implicit function theoremand the openness of hyperbolicity (Hartman-Grobman theorem), if (cid:15) > (cid:15) , the chain-recurrent set consists only of finitely many equilibria, which impliesconvergence of all trajectories.
Remark 6.6.
Boundary equilibria that are unstable in the absence of migration, candisappear under weak migration because they may leave S Γ I . A simple example is over-dominance in two diallelic demes: the unstable zero-migration equilibria ( p , , p , ) = (1 , p , , p , ) = (0 ,
1) do not survive perturbation. Indeed, if (cid:15) >
0, the perturbation of(1 ,
0) must have the form (1 − (cid:15)z , (cid:15)z ). If the fitnesses of the genotypes are as in (4.1)and the migration rates are m = (cid:15)m and m = (cid:15)m , then straightforward calculationsshow that this equilibrium is given by z = − x m − x and z = − y m − y , which is not in [0 , × [0 ,
1] if there is overdominance, i.e., if x α < y α < emark 6.7. In the absence of migration, mean fitness is monotone increasing in eachdeme (except at the equilibria); see (2.13). Therefore, by continuity, ∆ ¯ w ( p ) = (cid:80) α ∆ ¯ w α ( p · ,α ) > (cid:15) if p is bounded away from the set F of equilibria. If p is close to F , mean fitness may decrease. As an example assume two diallelic demes with overdomi-nance and (stable) equilibria ˆ p · , and ˆ p · , with ˆ p · , (cid:54) = ˆ p · , . Suppose that in some generation p · ,α = ˆ p · ,α holds for α = 1 ,
2. Since ¯ w ( p ) is maximized at ˆ p , which, with migration, is anequilibrium only if ˆ p · , = ˆ p · , , we see that ¯ w ( p (cid:48) ) < ¯ w ( p ).As another application of Theorem 6.5, we study the number of alleles that can bemaintained at an asymptotically stable equilibrium under weak migration. First we provea simple result for a single deme. For more general results, see Sect. 2 in Nagylaki andLou (2006a). 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 − . In addition, assume that w > w and there is intermediatedominance, i.e., w ii ≥ w ij ≥ w jj (6.27) for every i and every j > i . Then allele A is fixed as t → ∞ . The correspondingequilibrium is globally asymptotically stable.Proof. The assumptions imply w i ≥ w ii ≥ w Ii for every i ∈ I and w > w I or w I > w II .It follows that w = (cid:80) i w i p i ≥ (cid:80) i w Ii p i = w I and, if p (cid:54) = 0 and p I (cid:54) = 0, w > w I . (6.28)Therefore, we obtain (cid:18) p I p (cid:19) (cid:48) = p (cid:48) I p (cid:48) = p I w I p w < p I p (6.29)if p > p I >
0. Hence, p I ( t ) → t → ∞ provided p I < p >
0. Now wecan repeat this argument with p and p I − . Proceeding inductively, we obtain p ( t ) → t → ∞ . Theorem 6.9.
Suppose that migration is sufficiently weak and there is partial dominancein every deme (5.18) .(a) Generically, there is global convergence to an asymptotically stable equilibrium atwhich at most Γ alleles are present. Thus, the number of demes is a generic upper boundfor 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 segregatingat a globally asymptotically stable equilibrium. roof. Corollary 6.8 shows that, in the absence of migration, one allele is fixed in everydeme. The parameter combinations satisfying the assumptions of the theorem clearlyform an open set of all possible parameter combinations. Therefore, without migration,at most Γ alleles can be maintained at an asymptotically stable equilibrium, and this canbe achieved on an open set in the parameters space. Moreover, this equilibrium is globallyasymptotically 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 ineach 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 equilibriumcannot exceed the number of demes. If migration is much stronger than selection, we expect that rapid convergence to spatialquasi-homogeneity occurs. After this short initial phase, evolution should be approxi-mately panmictic with suitably averaged allele frequencies. Because in the absence ofselection, there exists a globally attracting manifold of equilibria, so that the dynamicsis not gradient like, the derivation of perturbation results is much more delicate than forweak migration. Since the fundamental ideas in the proofs of the most relevant resultsare essentially the same as if selection acts on many loci, we defer the analysis of strongmigration to Section 7.6, where the multilocus case is treated.
Selection is called uniform if w ij,α = w ij for every i , j , and every α . (6.30)Under spatially uniform selection, one might expect that population structure leaves nogenetic traces. However, this is not always true as shown by Example 4.5 with under-dominance and weak migration. Indeed, if we have two diallelic demes with the sameunderdominant selection in both, then under weak migration there are nine equilibria,four of which are asymptotically stable. These are the two monomorphic equilibria and59he two equilibria where each of the alleles is close to fixation in one deme and rare inthe other. Only three of the equilibria are uniform, i.e., have the same allele frequenciesin both demes. These are the two monomorphic equilibria and the ‘central’ equilibrium.In the following we state sufficient conditions under which there is no genetic indicationof population structure. We call ˆ p ∈ S Γ I a uniform selection equilibrium if every ˆ p · ,α is anequilibrium of the pure selection dynamics (6.21) and ˆ p · ,α = ˆ p · ,β for every α , β . Theorem 6.10 (Nagylaki and Lou 2007, Theorem 5.1) . If ˆ p ∈ S Γ I is a uniform selectionequilibrium, then ˆ p is an equilibrium of the (full) migration-selection dynamics (3.5) , and ˆ p 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 determinedentirely by selection and is independent of migration.Next, one may ask for the conditions when the solutions of (3.5) converge globally toa uniform selection equilibrium. One may expect global convergence under migration ifthe uniform selection equilibrium is globally asymptotically stable without migration. Sofar, only weaker results could be proved.
Theorem 6.11 (Nagylaki and Lou 2007, Sections 5.2 and 5.3) . Suppose there is a uni-form selection equilibrium that is globally asymptotically stable in the absence of migra-tion. Each of the following conditions implies global convergence of solutions to ˆ p undermigration and selection.(a) Migration is weak and, without migration, every equilibrium is hyperbolic.(b) Migration is strong, M is ergodic, and every equilibrium of the strong-migrationlimit is hyperbolic.(c) The continuous-time model (6.9) applies, µ αβ > for every pair α, β with α (cid:54) = β ,and ˆ p is internal in the absence of migration (hence, also with migration).(d) The continuous-time model (6.9) applies, µ αβ = µ βα for every pair α, β , and ˆ p isinternal. Since many phenotypic traits are determined by many genes, many loci are subject toselection. If loci are on the same chromosome, especially if they are within a short physicaldistance, they can not be treated independently. In order to understand the evolutionaryeffects of selection on multiple loci, models have to be developed and studied that take60inkage 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 specialassumptions. Whereas the former can be sometimes extended using perturbation theoryto obtain general insight, the latter are mainly useful to study specific biological questionsor to demonstrate the kind of complexity that can arise.Before developing the general model, we introduce the basic model describing the in-teraction 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.
We consider two loci, A and B , each with two alleles, A , A , and B , B . Therefore,there are four possible gametes, A B , A B , A B , A B , and 16 diploid genotypes. If, asusual, A i B j / A k B (cid:96) and A k B (cid:96) / A i B j are indistinguishable, then only 10 different unorderedgenotypes remain.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. Wedenote this recombination probability by r . The value of r usually depends on the distancebetween 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 = are called unlinked. Themaximum value of r = occurs for loci on different chromosomes, because then all fourgametes are produced with equal frequency . Thus, the recombination rate satisfies0 ≤ r ≤ .If, for instance, in the initial generation only the genotypes A B / A B and A B / A B are present, then in the next generation only these double homozygotes, as well as thetwo double heterozygotes A B / A B and A B / A B will be present. After furthergenerations of random mating, all other genotypes will occur, but not immediately attheir equilibrium frequencies. The formation of gametic types other than A B or A B 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 numberof alleles at each locus. Let the frequencies of the alleles A i at the first locus be denotedby p i and those of the alleles B j at the second locus by q j . Then p i = (cid:88) j P ij and q j = (cid:88) i P ij . (7.1)61he allele frequencies are no longer sufficient to describe the genetic composition of thepopulation because, in general, they do not evolve independently. Linkage equilibrium (LE) is defined as the state in which P ij = p i q j (7.2)holds for every i and j . Otherwise the population is said to be in linkage disequilibrium (LD). LD is equivalent to probabilistic dependence of allele frequencies between loci.Given P ij , we want to find the gametic frequencies P (cid:48) ij in the next generation afterrandom mating. The derivation of the recursion equation is based on the following basicfact of Mendelian genetics: an individual with genotype A i B j / A k B l produces gametesof parental type if no recombination occurs (with probability 1 − r ), and recombinantgametes if recombination between occurs (with probability r ). Therefore, the fraction ofgametes A i B j and A k B l is (1 − r ) each, and that of A i B l and A k B j is r each. Fromthese considerations, we see that the frequency of gametes of type A i B j in generation t + 1produced without recombination is (1 − r ) P ij , and that produced with recombination is rp i q j because of random mating. Thus, P (cid:48) ij = (1 − r ) P ij + rp i q j . (7.3)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 bythe parameter D ij = P ij − p i q j . (7.4)The D ij are called linkage disequilibria. From (7.3) and (7.4) we infer D (cid:48) ij = (1 − r ) D ij (7.5)and D ij ( t ) = (1 − r ) t D ij (0) . (7.6)Therefore, unless r = 0, linkage disequilibria decay at the geometric rate 1 − r and LE isapproached gradually without oscillation.For two alleles at each locus, it is more convenient to label the frequencies of thegametes A B , A B , A B , and A B by x , x , x , and x , respectively. A simplecalculation reveals that D = x x − x x (7.7)62 B A B A B A B Figure 7.1: The tetrahedron represents the state space S of the two-locus two-allelemodel. The vertices correspond to fixation of the labeled gamete, and frequencies aremeasured by the (orthogonal) distance from the opposite boundary face. At the centerof the simplex all gametes have frequency . The two-dimensional (red) surface is theLE manifold, D = 0. The states of maximum LD, D = ± , are the centers of the edgesconnecting A B to A B and A B to A B .satisfies D = D = − D = − D = D . (7.8)Thus, the recurrence equations (7.3) for the gamete frequencies can be rewritten as x (cid:48) i = x i − η i rD , i ∈ { , , , } , (7.9)where η = η = − η = − η .The two-locus gametic frequencies are the elements of the simplex S and may berepresented geometrically by the points in a tetrahedron. The subset where D = 0 formsa two-dimensional manifold and is called the linkage equilibrium, or Wright, manifold. Itis displayed in Figure 7.1.If r >
0, (7.6) implies that all solutions of (7.9) converge to the LE manifold alongstraight lines, because the allele frequencies x + x and x + x remain constant, wheresets such as x + x = const. represent planes in this geometric picture. The LE manifoldis invariant under the dynamics (7.9). 63 .2 Two diallelic loci under selection To introduce selection, we assume that viability selection acts on juveniles. Then recom-bination and random mating occurs. Since selection acts on diploid individuals, we assignfitnesses to two-locus genotypes. We denote the fitness of genotype ij ( i, j ∈ { , , , } )by w ij , where we assume w ij = w ji , because usually it does not matter which gameteis paternally or maternally inherited. The marginal fitness of gamete i is defined by w i = (cid:80) i =1 w ij x j , and the mean fitness of the population is ¯ w = (cid:80) i,j =1 w ij x i x j . If weassume, as is frequently the case, that there is no position effect, i.e., w = w , simpleconsiderations yield the selection-recombination dynamics (Lewontin and Kojima 1960): x (cid:48) i = x i w i ¯ w − η i w ¯ w rD , i ∈ { , , , } . (7.10)This is a much more complicated dynamical system than either the pure selection dynam-ics (2.6) or the pure recombination dynamics (7.9), and has been studied extensively (fora review, see B¨urger 2000, Sects. II.2 and VI.2). In general, mean fitness may decreaseand is no longer maximized at an equilibrium.In addition, the existence of stable limit cycles has been established for this discrete-time model (Hastings 1981b, Hofbauer and Iooss 1984) as well as for the correspondingcontinuous-time model (Akin 1979, 1982). Essentially, the demonstration of limit cyclesrequires that selection coefficients and recombination rates are of similar magnitude.There is a particularly important special case in which the dynamics is simple. Thisis the case of no epistasis , or additive fitnesses . Then there are constants u ( n ) i such that w ij = u (1) i + u (2) j for every i, j ∈ { , , , } . (7.11)In the absence of epistasis, i.e., if (7.11) holds, mean fitness ¯ w is a (strict) Lyapunovfunction (Ewens 1969). In addition, a point p is an equilibrium point of (7.10) if and onlyif it is both a selection equilibrium for each locus and it is in LE (Lyubich 1992, Nagylakiet al. 1999). In particular, the equilibria are the critical points of mean fitness.The reason for this increased complexity of two-locus (or multilocus) systems lies notso much in the increased dimensionality but arises mainly from the fact that epistaticselection generates nonrandom associations (LD) among the alleles at different loci. Re-combination breaks up these associations to a certain extent but changes gamete frequen-cies in a complex way. Thus, there are different kinds of interacting nonlinearities arisingin the dynamical equations under selection and recombination.64 .3 The general model We extend the migration-selection model of Section 3 by assuming that selection acts ona finite number of recombining loci. The treatment follows B¨urger (2009a), which wasinspired by Nagylaki (2009b). We consider a diploid population with discrete, nonover-lapping generations, in which the two sexes need not be distinguished. The population issubdivided into Γ ≥ L ≥ I n ≥ A ( n ) i n ( i n = 1 , . . . , I n ),at locus n . We use the multi-index i = ( i , . . . , i L ) as an abbreviation for the gamete A (1) i . . . A ( L ) i L . We designate the set of all loci by L = { , . . . , L } , the set of all alleles atlocus n by I n = { , . . . , I n } , and the set of all gametes by I . The number of gametes is I = | I | = (cid:81) n I n , the total number of genes (alleles at all loci) is I + · · · + I L . We useletters i, j, (cid:96) ∈ I for gametes, k, n ∈ L for loci, and α, β ∈ G for demes. Sums or productswithout ranges indicate summation over all admissible indices, e.g., (cid:80) n = (cid:80) n ∈ L .Let p i,α = p i,α ( t ) represent the frequency of gamete i among zygotes in deme α ingeneration t . We define the following column vectors: p i = ( p i, , . . . , p i, Γ ) (cid:62) ∈ R Γ , (7.12a) p · ,α = ( p ,α , . . . , p I,α ) (cid:62) ∈ S I , (7.12b) p = (cid:0) p (cid:62)· , , . . . , p (cid:62)· , Γ (cid:1) (cid:62) ∈ S Γ I . (7.12c)Here, p i , p · ,α , and p signify the frequency of gamete i in each deme, the gamete frequenciesin deme α , and all gamete frequencies, respectively. We will use analogous notation forother quantities, e.g. for D i,α .The frequency of allele A ( k ) i k among gametes in deme α is p ( k ) i k ,α = (cid:88) i | i k p i,α , (7.13)where the sum runs over all multi-indices i with the k th component fixed as i k . We write p ( k ) i k = (cid:16) p ( k ) i k , , . . . , p ( k ) i k , Γ (cid:17) (cid:62) ∈ R Γ (7.14)for the vector of frequencies of allele A ( k ) i k in each deme.65et x ij,α and w ij,α denote the frequency and fitness of genotype ij in deme α , respec-tively. We designate the marginal fitness of gamete i in deme α and the mean fitness ofthe population in deme α by w i,α = w i,α ( p · ,α ) = (cid:88) j w ij,α p j,α (7.15a)and ¯ w α = ¯ w α ( p · ,α ) = (cid:88) i,j w ij,α p i,α p j,α . (7.15b)The life cycle starts with zygotes in Hardy-Weinberg proportions. Selection acts in eachdeme on the newly born offspring. Then recombination occurs followed by adult migrationand random mating within in each deme. This life cycle extends that in Section 3.2. Todeduce the general multilocus migration-selection dynamics (Nagylaki 2009b), let x ∗ ij,α = p i,α p j,α w ij,α / ¯ w α (7.16a)be the frequency of genotype ij in deme α after selection, and p i,α = (cid:88) j,(cid:96) R i,j(cid:96) x ∗ j(cid:96),α (7.16b)its frequency after recombination. Here, R i,j(cid:96) is the probability that during gametogenesis,paternal haplotypes j and (cid:96) produce a gamete i by recombination. Finally, migrationoccurs and yields the gamete frequencies in the next generation in each deme: p (cid:48) i,α = (cid:88) β m αβ p i,β . (7.16c)The recurrence equations (7.16) describe the evolution of gamete frequencies under selec-tion on multiple recombining loci and migration. We view (7.16) as a dynamical systemon S Γ I . We leave it to the reader to check the obvious fact that the processes of migrationand recombination commute.The complications introduced by recombination are disguised by the terms R i,jl whichdepend on the recombination frequencies among all subsets of loci. To obtain an analyt-ically useful representation of (7.16b), more effort is required.Let { K , N } be a nontrivial decomposition of L , i.e., K and its complement N = L \ K areeach proper subsets of L and, therefore, contain at least one locus. (The decompositions66 K , N } and { N , K } are identified.) We designate by c K the probability of reassociationof the genes at the loci in K , inherited from one parent, with the genes at the loci in N ,inherited from the other. Let c tot = (cid:88) K c K , (7.17)where (cid:80) K runs over all (different) decompositions { K , N } of L , denote the total recombi-nation frequency. We designate the recombination frequency between loci k and n , suchthat k < n , by c kn . It is given by c kn = (cid:88) K ∈ L kn c K , (7.18)where L kn = { K : k ∈ K and n ∈ N } (B¨urger 2000, p. 55). Unless stated otherwise, weassume that all pairwise recombination rates c kn are positive. Hence, c min = min k
1. Migration and recombination rates, m αβ and c K , are fixed so that fitness differences are small compared with them. From (7.15)and (7.32), we deduce w i,α ( p · ,α ) = 1 + (cid:15)r i,α ( p · ,α ) , ¯ w α ( p · ,α ) = 1 + (cid:15) ¯ r α ( p · ,α ) , (7.33)in which r i,α ( p · ,α ) = (cid:88) j r ij,α p j,α , ¯ r α ( p · ,α ) = (cid:88) i,j r ij,α p i,α p j,α . (7.34)When selection is dominated by migration and recombination, we expect that linkagedisequilibria within demes as well as gamete- and gene-frequency differences betweendemes decay rapidly to small quantities. In particular, we expect approximately panmicticevolution of suitably averaged gamete frequencies in ‘quasi-linkage equilibrium’. We alsoshow 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 by Nagylaki et al.(1999, Theorem 3.1). For a single locus under selection and strong migration, this is thecontent of Theorem 4.5 in Nagylaki and Lou (2007). Theorem 7.2 and its proof combineand extend these results as well as some of the underlying ideas and methods.To formulate and prove this theorem, we define the vector ρ α = (cid:16) p (1)1 ,α , . . . , p (1) I ,α , . . . , p ( L )1 ,α , . . . , p ( L ) I L ,α (cid:17) T ∈ S I × · · · × S I L (7.35)of all allele frequencies at every locus in deme α , and the vector π = (cid:16) P (1)1 , . . . , P (1) I , . . . , P ( L )1 , . . . , P ( L ) I L (cid:17) T ∈ S I × · · · × S I L (7.36)of all averaged allele frequencies at every locus. We note that in the presence of selectionthe P ( k ) i k , hence π , are time dependent. Instead of p , we will use π , D , and q to analyze(7.16), and occasionally write p = ( π, D, q ).70n the LE manifold Λ ,α (7.21), which is characterized by the ρ α ( α ∈ G ), the selectioncoefficients of gamete i , allele i n at locus n , and of the entire population are r i,α ( ρ α ) = (cid:88) j r ij,α (cid:89) k p ( k ) j k ,α , (7.37a) r ( n ) i n ,α ( ρ α ) = (cid:88) i | i n r i,α ( ρ α ) (cid:89) k : k (cid:54) = n p ( k ) i k ,α , (7.37b)¯ r α ( ρ α ) = (cid:88) i r i,α ( ρ α ) (cid:89) k p ( k ) i k , (7.37c)cf. (7.34). As in (7.24), let ξ denote the principal left eigenvector of M . We introducethe average selection coefficients of genotype ij , gamete i , allele i n at locus n , and of theentire population: ω ij = (cid:88) α ξ α r ij,α , (7.38a) ω i ( π ) = (cid:88) j ω ij (cid:89) k P ( k ) j k = (cid:88) α ξ α r i,α ( π ) , (7.38b) ω ( n ) i n ( π ) = (cid:88) i | i n ω i ( π ) (cid:89) k (cid:54) = n P ( k ) i k = (cid:88) α ξ α r ( n ) i n ,α ( π ) , (7.38c)¯ ω ( π ) = (cid:88) i ω i ( π ) (cid:89) k P ( k ) i k = (cid:88) α ξ α ¯ r α ( π ) . (7.38d)For ¯ ω , we obtain the alternative representations¯ ω ( π ) = (cid:88) n (cid:88) i n ω ( n ) i n P ( n ) i n = (cid:88) i,j ω ij (cid:32)(cid:89) n P ( n ) i n (cid:33) (cid:32)(cid:89) k P ( k ) j k (cid:33) , (7.38e)and d¯ ω ( π )d P ( n ) i n = 2 ω ( n ) i n ( π ) . (7.39)For reasons that will be justified by the following theorem, we call the differentialequation ˙ P ( n ) i n = P ( n ) i n (cid:104) ω ( n ) i n ( π ) − ¯ ω ( π ) (cid:105) , (7.40a) D = 0 , q = 0 (7.40b)on S Γ I the weak-selection limit of (7.16). In view of the following theorem, it is moreconvenient to consider (7.40a) and (7.40b) on S Γ I instead of (7.40a) on S I × · · · × S I L .The differential equation (7.40a) is a Svirezhev-Shashahani gradient (Remark 2.3) with71otential function ¯ ω . In particular, ¯ ω increases strictly along nonconstant solutions of(7.40a) because ˙¯ ω = 2 (cid:88) n (cid:88) i n P ( n ) i n (cid:104) ω ( n ) i n ( π ) − ¯ ω ( π ) (cid:105) ≥ . (7.41)We will also need the assumptionAll equilibria of (7.40a) are hyperbolic. (H) Theorem 7.2 (B¨urger 2009a) . Suppose that (7.16) , (7.32) , (E) and (H) hold, the back-ward migration matrix M and all recombination rates c K are fixed, and (cid:15) > is sufficientlysmall.(a) The set of equilibria Ξ ⊂ S Γ I of (7.40) contains only isolated points, as does theset of equilibria Ξ (cid:15) ⊂ S Γ I of (7.16) . As (cid:15) → , each equilibrium in Ξ (cid:15) converges to thecorresponding equilibrium in Ξ .(b) In the neighborhood of each equilibrium in Ξ , there exists exactly one equilibriumpoint in Ξ (cid:15) . The stability of each equilibrium in Ξ (cid:15) is the same as that of the correspondingequilibrium in Ξ ; i.e., each pair is either asymptotically stable or unstable.(c) Every solution p ( t ) of (7.16) converges to one of the equilibrium points in Ξ (cid:15) . The essence of this theorem and its proof is that under weak selection (i) the exactdynamics 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.40). The latter is much easier to study because it is formally equivalentto a panmictic one-locus selection dynamics. Theorem 7.2 is a singular perturbation resultbecause in the absence of selection every point on Ψ is an equilibrium.Parts (a) and (b) of the above theorem follow essentially from Theorem 4.4 of Karlinand McGregor (1972b) which is an application of the implicit function theorem and theHartman-Grobman theorem. Part (c) is much stronger and relies, among others, on thenotion 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 hyperbolicmanifolds imply that for sufficiently small s , there exists a smooth invariant manifold Ψ (cid:15) close to Ψ , and Ψ (cid:15) is globally attracting for (7.16) at a geometric rate (see Nagylaki etal. 1999, and the references there). The manifold Ψ (cid:15) is characterized by an equation ofthe form ( D, q ) = (cid:15)ψ ( π, (cid:15) ) , (7.42)72here ψ is a smooth function of π . Thus, on Ψ (cid:15) , and more generally, for any initial values,after a long time, D ( t ) = O ( (cid:15) ) and q ( t ) = O ( (cid:15) ) . (7.43)The next step consists in deriving the recurrence equations in an O ( (cid:15) ) neighborhoodof Ψ which, in particular, contains Ψ (cid:15) . By applying (7.32) and D ( t ) = O ( (cid:15) ) to (7.16b’),straightforward calculations yield p ( n ) i n ,α = p ( n ) i n ,α + (cid:15)p ( n ) i n ,α r ( n ) i n ,α ( ρ α ) − ¯ r α ( ρ α )¯ w α ( ρ α ) + O ( (cid:15) ) (7.44)for every α ∈ G (see eq. (3.10) in Nagylaki et al. 1999). By averaging, invoking (7.16),(7.44), and (7.42) one obtains P ( n ) i n (cid:48) = µ T p ( n ) i n (cid:48) = µ T M p ( n ) i n = µ T p ( n ) i n = P ( n ) i n + (cid:15)P ( n ) i n (cid:104) ω ( n ) i n ( π ) − ¯ ω ( π ) (cid:105) + O ( (cid:15) ) . (7.45)By rescaling time as τ = (cid:15)t and letting (cid:15) →
0, the leading term in (7.45), P ( n ) i n (cid:48) = P ( n ) i n + (cid:15)P ( n ) i n (cid:104) ω ( n ) i n ( π ) − ¯ ω ( π ) (cid:105) , approximates the gradient system (7.40a). Therefore, wehave ¯ ω ( π (cid:48) ) > ¯ ω ( π ) unless π (cid:48) = π . In particular, the dynamics (7.45) on Ψ is gradientlike.The eigenvalues of the Jacobian of (7.45) have the form 1 + (cid:15)ν + O ( (cid:15) ), where ν signifies an eigenvalue of the Jacobian of (7.40a). Therefore, (H) implies that also everyequilibrium of (7.45) is hyperbolic.The rest of the proof is identical to that of Theorem 3.1 in Nagylaki et al. (1999). Theheart of its proof is the following. Since solutions of (7.16) are in phase with solutions onthe invariant manifold Ψ (cid:15) , it is sufficient to prove convergence of trajectories for initialconditions p ∈ Ψ (cid:15) . With the help of the qualitative theory of numerical approximations,it can be concluded that the chain-recurrent set of the exact dynamics (7.16) on Ψ (cid:15) is asmall perturbation of the chain-recurrent set of (7.40a). The latter dynamics, however,is gradient like. Because all equilibria are hyperbolic by assumption, the chain-recurrentset consists exactly of those equilibria. Therefore, the same applies to the chain-recurrentset of (7.16) if (cid:15) > p of (7.40) on the boundary of S Γ I satisfies ˆ p i = 0 for every i in some subset of I , and this condition is not altered bymigration. The positive components of ˆ p are perturbed within the boundary, where theequilibrium remains. 73ur next result concerns the average of the (exact) mean fitnesses over demes:¯ w ( p ) = (cid:88) α ξ α ¯ w α ( p · ,α ) . (7.46) Theorem 7.3 (B¨urger 2009a) . Suppose the assumptions of Theorem 7.2 apply. If (7.43) holds, π is bounded away from the equilibria of (7.40a) , and p is within O ( (cid:15) ) of Ψ (cid:15) , then ∆ ¯ w ( p ) > . The time to reach an O ( (cid:15) ) neighborhood of Ψ (cid:15) can be estimated (Remark 4.6 inB¨urger 2009a and Remark 4.13 in Nagylaki and Lou 2007). Unless migration is veryweak or linkage very tight, convergence occurs in an evolutionary short period. It followsthat mean fitness is increasing during the long period when most allele-frequency changeoccurs, i.e., after reaching Ψ (cid:15) and before reaching a small neighborhood of the stableequilibrium.An essential step in the proof is to show that( D (cid:48) , q (cid:48) ) − ( D, q ) = O ( (cid:15) ) (7.47)is satisfied if (7.43) holds. Therefore, linkage disequilibria and the measure q of spatialdiversity change very slowly on Ψ (cid:15) . This justifies to call states on Ψ (cid:15) spatially quasi-homogeneous and to be in quasi-linkage equilibrium.Theorem 7.2 enables the derivation of the following result about the equilibrium struc-ture of the multilocus migration-selection dynamics (7.16). It establishes that arbitrarilymany loci with arbitrarily many alleles can be maintained polymorphic by migration-selection balance: Theorem 7.4.
Let L ≥ , Γ ≥ , I n ≥ 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 forevery parameter combination in Q , there is a unique, internal, asymptotically stable equi-librium point. This equilibrium is spatially quasi-homogeneous, is in quasi-linkage equi-librium, and attracts all trajectories with internal initial condition. Furthermore, everytrajectory converges to an equilibrium point as t → ∞ .(b) Such an open set, Q (cid:48) , also exists if the set of all fitnesses is restricted to benonepistatic and to display intermediate dominance. For a more detailed formulation and a proof, see Theorem 2.2 in B¨urger (2009b). Theconstructive proof suggests that with increasing number of alleles and loci, the proportion74f parameter space that allows for a fully polymorphic equilibrium shrinks rapidly becauseat every locus some form of overdominance of spatially averaged fitnesses is required. Inaddition, the proof shows that this set Q exists in the parameter region where migrationis strong relative to selection. Remark 7.5.
The set Q (cid:48) in the above theorem does not contain parameter combinationssuch that all one-locus fitnesses are multiplicative or additive, or more generally satisfyDIDID (5.26). In each of these cases, one gamete becomes fixed as t → ∞ (Proposition2.6 and Remark 2.7 in B¨urger 2009b). Now we assume that selection and recombination are weak relative to migration, i.e., inaddition to (7.32), we posit c K = (cid:15)γ K , K ⊆ L , (7.48)where (cid:15) ≥ γ K is defined by this relation. Then every solution p ( t ) of the full dynamics (7.16) converges to a manifold close to q = 0. On this manifold,the dynamics is approximated by the differential equation˙ P i = P i [ ω i ( P ) − ¯ ω ( P )] − (cid:88) K γ K (cid:104) P i − P ( K ) i K P ( N ) i N (cid:105) , (7.49a)which we augment with q = 0 . (7.49b)This is called the strong-migration limit of (7.16).In general, it cannot be expected that the asymptotic behavior of solutions of (7.16)under strong migration is governed by (7.49) because its chain-recurrent set does not al-ways consist of finitely many hyperbolic equilibria. The differential equation (7.49a) is thecontinuous-time version of the discrete-time dynamics (7.16b’) which describes evolutionin a panmictic population subject to multilocus selection and recombination. Althoughthe recombination term is much simpler than in the corresponding difference equation,the dynamics is not necessarily less complex. Akin (1979, 1982) proved that (7.49) mayexhibit stable cycling. Therefore, under strong migration and if selection and recombina-tion are about equally weak, convergence of trajectories of (7.16) will not generally occur,and only local perturbation results can be derived (B¨urger 2009a, Proposition 4.10).The dynamics (7.49) becomes simple if there is a single locus (then the recombinationterm in (7.49a) vanishes) or if there is no epistasis, i.e., if fitnesses can be written in the75orm ω ij = (cid:88) n u ( n ) i n j n (7.50)for constants u ( n ) i n j n ≥
0. Essentially, this means that fitnesses can be assigned to eachone-locus genotype and there is no nonlinear interaction between loci. Then mean fitnessis a Lyapunov function, all equilibria are in LE, and global convergence of trajectoriesoccurs generically (Ewens 1969, Lyubich 1992, Nagylaki et al. 1999).
If recombination is weak relative to migration and selection, the limiting dynamics isformally equivalent to a multiallelic one-locus migration-selection model, as treated inSections 3 – 6. Then local perturbation results can be applied but, in general, globalconvergence of trajectories cannot be concluded and, in fact, does not always occur.
In contrast to the one-locus case treat in in Section 6.3, in the multilocus case the as-sumption of weak migration is insufficient to guarantee convergence of trajectories. Thereason is that in the absence of migration the dynamics (7.16) reduces to p (cid:48) i,α = p i,α = p i,α w i,α ¯ w α − D i,α (7.51)for every i ∈ I and every α ∈ G ; cf. (7.16b’). Therefore, we have Γ decoupled multilocusselection dynamics, one for each deme. Because already for a single deme, stable cyclinghas been established (Hastings 1981, Hofbauer and Iooss 1984), global perturbation re-sults can not be achieved without additional assumptions. Such an assumption is weakepistasis.We say there is weak epistasis if we can assign (constant) fitness components u ( n ) i n j n ,α > w ij,α = (cid:88) n u ( n ) i n j n ,α + ηs ij,α , (7.52)where the numbers s ij,α satisfy | s ij,α | ≤ η ≥ η is small enough, so that w ij,α > η = η ( (cid:15) ) , (7.53)76here η : [0 , → [0 , ∞ ) is C and satisfies η (0) = 0. Therefore, migration and epistasisneed not be ‘equally’ weak. In particular, no epistasis ( η ≡
0) is included.In the absence of epistasis and migration ( (cid:15) = η = 0), p is an equilibrium point of(7.51) if and only if for every α ∈ G , p · ,α is both a selection equilibrium for each locus andis in LE (Lyubich 1992, Nagylaki et al. 1999). In addition, the only chain-recurrent pointsof (7.51) are its equilibria (Lemmas 2.1 and 2.2 in Nagylaki et al. 1999). Therefore, onecan apply the proof of Theorem 2.3 in Nagylaki et al. (1999) to deduce the following resultwhich simultaneously generalizes Theorem 2.3 in Nagylaki et al. (1999) and Theorem 6.5: Theorem 7.6 (B¨urger 2009a, Theorem 5.4) . Suppose that in the absence of epistasisevery equilibrium of (7.51) is hyperbolic, and (cid:15) > is sufficiently small.(a) The set of equilibria Σ ⊂ S Γ I of (6.21) contains only isolated points, as does theset of equilibria Σ (cid:15) ⊂ S Γ I of (7.16) . As (cid:15) → , each equilibrium in Σ (cid:15) converges to thecorresponding equilibrium in Σ .(b) In the neighborhood of each asymptotically stable equilibrium in Σ , there existsexactly one equilibrium in Σ (cid:15) , and it is asymptotically stable. In the neighborhood of eachunstable internal equilibrium in Σ , there exists exactly one equilibrium in Σ (cid:15) , and it isunstable. In the neighborhood of each unstable boundary equilibrium in Σ , there exists atmost one equilibrium in Σ (cid:15) , and if it exists, it is unstable.(c) Every solution p ( t ) of (7.16) converges to one of the equilibrium points in Σ (cid:15) . In contrast to the case of weak selection (Theorem 7.2), unstable boundary equilibriacan leave the state space under weak migration (Karlin and McGregor 1972a). For anexplicit example in the one-locus setting, see Remark 4.2 if Nagylaki and Lou (2007).If (cid:15) = 0, then all equilibria are in Λ , i.e., there is LE within each deme, but notbetween demes. If (cid:15) > D i,α = O ( (cid:15) ) for every i and every α . Remark 7.7.
If in Theorem 7.6 the assumption of weak epistasis is not made, statements(a) and (b) remain valid, but (c) does not necessarily follow. Statement (c) follows if, for (cid:15) = 0, the chain recurrent set consists precisely of the hyperbolic equilibria.One of the consequences of the above theorem is the following, which contrasts sharplywith Theorem 7.4:
Theorem 7.8 (B¨urger 2009b, Theorem 3.9) . For an arbitrary number of loci, sufficientlyweak migration and epistasis, and partial dominance, the number of demes is the generic aximum for the number of alleles that can be maintained at any locus at any equilibrium(stable or not) of (7.16) . If all evolutionary forces are weak, i.e., if (cid:15) → c K = (cid:15)γ K for all subsets K ⊆ L , (7.54)the limiting dynamics of (7.16) on S Γ I becomes˙ p i,α = p i,α [ r i,α ( p · ,α ) − ¯ r α ( p · ,α )] − (cid:88) K γ K (cid:16) p i,α − p ( K ) i K ,α p ( N ) i N ,α (cid:17) + (cid:88) β µ αβ p i,β . (7.55)Hence, selection, recombination, and migration are decoupled. This dynamics may beviewed as the continuous-time version of (7.16). However, it may exhibit complex dynam-ical behavior already if either migration or recombination is absent. Analogs of Theorems7.1, 7.2, 7.3, 7.4, and 7.6 apply to (7.55). For the Levene model, the recurrence equations (7.16) simplify to p (cid:48) i = (cid:88) j,(cid:96),α R i,j(cid:96) p j p (cid:96) c α w j(cid:96),α / ¯ w α , i ∈ I , (7.56)where this form is sufficient to deduce the main results, i.e., it is not necessary to expressthe recombination probabilities R i,j(cid:96) in terms of the linkage disequilibria.Many of the results on the one-locus Levene model in Sections 5.1 and 5.2 can be gener-alized to the multilocus Levene model if epistasis is absent or weak. These generalizationsare based on two key results. The first is that in the absence of epistasis geometric meanfitness is again a Lyapunov function, and the internal equilibria are its stationary points(Theorems 3.2 and 3.3 in Nagylaki 2009b). The second key result is that in the absenceof epistasis, generically, every trajectory converges to an equilibrium point that is in LE(Theorem 3.1 in B¨urger 2010). This is true for the haploid and the diploid model. Then aproof analogous to that of Theorem 7.6 yields a global perturbation result for weak epis-tasis, in particular, generic global convergence to an equilibrium point in quasi-linkageequilibrium (Theorem 7.2 in B¨urger 2010).With the aid of these results, the next two theorems can be derived about the main-tenance of multilocus polymorphism in the Levene model.78 heorem 7.9 (B¨urger 2010, Result 5.1) . Assume an arbitrary number of multiallelicloci, Γ ≥ demes of given size, and let all recombination rates be positive and fixed.Then there exists a nonempty open set of fitness parameters, exhibiting partial dominancebetween every pair of alleles at every locus, such that for every parameter combination inthis set, there is a unique, internal, asymptotically stable equilibrium point of the dynamics (7.56) . This equilibrium is globally asymptotically stable. Theorem 7.9 is the multilocus extension of Theorem 5.8. From a perturbation theoremfor weak epistasis in the Levene model, analogous to Theorem 7.6 (Theorem 7.2 in B¨urger2010), it follows that such an open set exists also within the set of non-epistatic fitnesses.The following result is of very different nature and generalizes Proposition 3.18 inNagylaki (2009b):
Theorem 7.10 (B¨urger 2010, Theorem 7.4) . Assume weak epistasis and diallelic lociwith DIDID.(a) If L ≥ Γ , at most Γ − loci can be maintained polymorphic for an open set ofparameters.(b) If L ≤ Γ − and if the dominance coefficients in (5.26) are arbitrary for each locusbut fixed, there exists an open nonempty set W of fitness schemes satisfying (7.52) and ofdeme proportions ( c , . . . , c Γ − ) , such that for every parameter combination in W , thereis a unique, internal, asymptotically stable equilibrium point of (7.56) . This equilibriumis in quasi-linkage equilibrium and globally attracting. The set W is independent of thechoice of the recombination rates. Whereas Theorem 7.9 suggests that spatially varying selection has the potential tomaintain considerable multilocus polymorphism, Theorem 7.10 shows that properties ofthe genetic architecture of the trait (no dominance or DIDID) may greatly constrain thispossibility. Numerical results in B¨urger (2009c) quantify how frequent polymorphism isin the two-locus two-allele Levene model, and how this depends on the selection schemeand 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 maintainedwill certainly decrease rapidly.A few other aspects of the multilocus Levene model have also been investigated. Zhiv-otovsky et al. (1996) employed a multilocus Levene model to study the evolution of phe-notypic plasticity. Wiehe and Slatkin (1998) explored a haploid Levene model in whichLD is caused by epistasis. More recently, van Doorn and Dieckmann (2006) performed79 numerical study that admits new mutations (of variable effect) and substitutions atmany loci. They showed that genetic variation becomes increasingly concentrated on afew loci. Roze and Rousset (2008) derived recursions for the allele frequencies and forvarious types of genetic associations in a multilocus infinite-island model. Barton (2010)studied certain issues related to speciation using a generalized haploid multilocus Levenemodel that admits habitat preferences.
Local adaptation of subpopulations to their specific environment occurs only if the allelesor genotypes that perform well are maintained within the population. Similarly, differen-tiation between subpopulations occurs only if they differ in allele or gamete frequencies.Therefore, maintenance of polymorphism is indispensable for evolving or maintaining lo-cal adaptation and differentiation. The more loci and alleles are involved, the higher is thepotential for local adaptation and differentation. This is the main reason, why conditionsfor the maintenance of polymorphism played such an important role in this survey.Because for a single randomly mating population, polymorphic equilibria do not existin the absence of epistasis and of overdominance or underdominance, it is natural to con-fine attention to the investigation of the maintenance of genetic variation in subdividedpopulations if in each subpopulation epistasis is absent (or weak) and dominance is inter-mediate. In addition, intermediate dominance seems to be by far the most common formof dominance.We briefly recall the most relevant results concerning maintenance of genetic variationand put them into perspective. Particularly relevant results for a single locus are pro-vided in Sections 4.2 – 4.5, and by Theorems 5.9, 5.16, 5.19, 6.1, and 6.9. Theorems 7.4,7.8, 7.9, and 7.10 treat multiple loci. Theorems 7.4 and 7.9 establish that in structuredpopulations and in the Levene model, respectively, spatially varying selection can main-tain multiallelic polymorphism at arbitrarily many loci under conditions for which in apanmictic population no polymorphism at all can be maintained.Interestingly, for strong migration and only two demes, an arbitrary number of allelescan be maintained at each of arbitrarily many loci, whereas for weak migration, thenumber of demes is a generic upper bound for the number of alleles that can be maintained(Theorem 7.8). At first, this appears counterintuitive given the widespread opinion thatit is easier to maintain polymorphism under weak migration than under strong migration.This opinion derives from the fact that for a single diallelic locus and two demes, the80arameter region for which a protected polymorphism exists increases with decreasingmigration rate, provided there is a single parameter that measures the strength of selec-tion, as is the case in the Deakin model; see condition (4.19). Additional support forthis opinion comes from the study of weak migration in homogeneous and heterogeneousenvironments (Karlin and McGregor 1972a,b; Christiansen 1999), from the analysis of thegeneral Deakin (1966) model (Section 4.3), as well as from various numerical studies (e.g.,Spichtig and Kawecki 2004, Star et al. 2007a,b). However, Karlin (1982, p. 128) observedthat for the non-homogeneous Deakin model, in which the ‘homing probabilities’ varyamong demes, “it is possible to increase a single homing rate and reduce or even abro-gate the event of A -protection”. This is already obvious from the protection condition(4.14): If only one of the migration rates in (4.13) is reduced, it depends on the sign ofthe corresponding selection coefficient if protection is facilitated or not.Theorem 7.4 is not at variance with Theorem 7.8 or the widespread opinion expressedabove. It is complimentary, and its proof yields deeper insight into the conditions underwhich genetic variation can be maintained by migration-selection balance. With strongmigration, a stable multiallelic polymorphism requires some form of overdominance ofsuitably averaged fitnesses for the loci maintained polymorphic. The reason is that strongmigration leads to strong mixing, so that gamete and allele frequencies become similaramong demes. Because the fraction of volume of parameter space, in which averageoverdominance holds for multiple alleles, shrinks rapidly with increasing number of allelesor loci, we expect very stringent constraints on selection and dominance coefficients formaintaining polymorphism at many loci if migration is strong.The conditions for maintaining loci polymorphic under weak migration are muchweaker. With arbitrary intermediate dominance, the proof of Theorem 7.8 shows thatthose 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 in theabsence of epistasis and if dominance is intermediate in every deme (Theorem 7.9). Butthere, the maximum number of polymorphic loci depends on the pattern of dominanceand the number of demes (Nagylaki 2009b, B¨urger 2010, and Theorem 7.10). It is an openproblem if for every (ergodic) migration scheme, arbitrarily many loci can be maintainedpolymorphic in the absence of epistasis, overdominance, and underdominance.Although the fraction of parameter space, in which the conditions for maintenanceof multiallelic multilocus polymorphisms are satisfied, decreases rapidly, stable multial-lelic polymorphism involving small or moderate numbers of alleles or loci does not seem81nlikely. What is required if dispersal is weak, basically, is that there is a mosaic of direc-tional selection pressures and different genes or genotypes that are locally well adapted. Here we show how to apply some of the above perturbation results in conjunction withLyapunov functions and an index theorem to derive the complete equilibrium and stabilitystructure for a two-locus continent-island model, in which epistatic selection may bestrong. The model was developed and studied by Bank et al. (2012) to explore how muchgene 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 intotwo, for instance because at least one of the subpopulations moves to a different habi-tat, and the two subpopulations remain separated for a sufficiently long time, in each ofthem 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 highfrequency), whereas in the other population aB may become fixed. If individuals from thederived populations mate, then hybrids of type AB , which occur by recombination be-tween aB and Ab have sometimes reduced fitness, i.e., they are incompatible (to a certaindegree). Therefore, the accumulation of DMIs is considered a plausible mechanism for theevolution of so-called intrinsic postzygotic isolation between allopatric (i.e., geographicallyseparated) 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, Bank et al. (2012) assumed a continent-island model, i.e., acontinental population, which is fixed for one haplotype (here, aB ), sends migrants toan island population in which, before the onset of gene flow, Ab was most frequent orfixed. To examine how much gene flow is needed to annihilate differentiation between thetwo populations, the equilibrium and stability structure was derived. A stable internalequilibrium, at which all four haplotypes are maintained, is called a DMI because, onlywhen 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, andthey investigated a haploid and a diploid version of their model. Here, we deal only with82he haploid case, as it admits a much more complete analysis and is also representativefor 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 inthe Online Supplement S1 of Bank et al. (2012). Also comprehensive biological motivationand discussion is provided by Bank et al., as well as illuminating figures. Let x , x , x , and x denote the frequencies of the four haplotypes ab , aB , Ab , and AB on the island. They satisfy x i ≥ i and (cid:80) i =1 x i = 1. We assume that selectionacts on individuals during the haploid phase of their life cycle according to the followingscheme for Malthusian fitness values:wild type continental island recombinanthaplotype ab aB Ab AB fitness w = 0 w = β w = α w = α + β − γ frequencies x x x x Table 8.1: Haplotype fitnesses and frequencies in the DMI modelThis scheme is entirely general. The fitness of the wild type is arbitrarily normalizedto 0. The parameters α and β measure a potential selective advantage of the island andcontinental haplotypes, respectively, on the island. They can be positive or negative.However, further below we will assume that α > γ measures the strength of theincompatibility among the A and B alleles. We assume γ ≥
0, so that epistasis is negativeor absent ( γ = 0).Assuming weak evolutionary forces in continuous time, the haplotype dynamics is givenby ˙ x i = x i ( w i − ¯ w ) − η i rD − mx i + δ i m, (8.1)where r is the recombination rate between the two loci, D = x x − x x is the LD, η = η = − η = − η , and δ i = 1 if i = 2 and δ i = 0 otherwise; cf. (7.10) and (7.55).83ith the settings as in Table 8.1, we obtain˙ x = x [ − α ( x + x ) − β ( x + x ) + γx ] − rD − mx , (8.2a)˙ x = x [ − α ( x + x ) + β ( x + x ) + γx ] + rD + m (1 − x ) , (8.2b)˙ x = x [ α ( x + x ) − β ( x + x ) + γx ] + rD − mx , (8.2c)˙ x = x [ α ( x + x ) + β ( x + x ) − γ ( x + x + x )] − rD − mx . (8.2d)This is a dynamical system on the simplex S which constitutes our state space. Wealways assume m ≥ r ≥
0, and γ ≥ p = x + x of A , q = x + x of B , and the measure D of LD. Then thedynamical equations read˙ p = αp (1 − p ) − γ (1 − p )( pq + D ) + βD − mp , (8.3a)˙ q = βq (1 − q ) − γ (1 − q )( pq + D ) + αD + m (1 − q ) , (8.3b)˙ D = [ α (1 − p ) + β (1 − q )] D − γ [(1 − p )(1 − q ) − D ]( pq + D ) − rD − m [ p (1 − q ) + D ] . (8.3c)It is straightforward to show that the condition x ∈ S translates to 0 ≤ p ≤
0, 0 ≤ q ≤ − min[ pq, (1 − p )(1 − q )] ≤ D ≤ min[ p (1 − q ) , (1 − p ) q ] . (8.4) We denote the monomorphic equilibria x i = 1 by M i . If m = 0, then all monomorphicequilibria exist. However, if α > γ > β , the conditions most relevant for thisinvestigation (see (8.17) below), M and M are always unstable.It is straightforward to determine local stability of M and M . The latter can bestable only if m = 0.Next, there may exist two equilibria at which one locus is polymorphic and the other isfixed. The equilibrium S A has the coordinates ( p, q, D ) = (cid:0) − mα − γ , , (cid:1) and is admissibleif and only if m < α − γ . The equilibrium S B has coordinates (cid:0) , − mβ , (cid:1) and is admissibleif and only if m < − β . The respective local stability conditions are easy to derive.The following lemma collects the important observations from these analyses.84 emma 8.1. (a) For given m > , at most one of the boundary equilibria M , S A , or S B can be stable. M is asymptotically stable if m > max[ − β, α − γ, α − β − r ] . (8.5) S A is asymptotically stable if γ < α + β and ( α − γ )( γ − β ) α (cid:18) α + β − γr (cid:19) < m < α − γ , (8.6) which requires r > γ − β . S B is asymptotically stable if γ > α + β and − βαγ − β (cid:18) γ − β − αr (cid:19) < m < − β , (8.7) which requires r > α .(b) If M 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 changein stability occurs, then it is from unstable to stable (as m increases). Finally, if r = 0, there is a fully polymorphic equilibrium R on the edge x = x = 0of S . Thus, only the island and the continental haplotypes are present. Its coordinatesare easily calculated.In the following we derive global asymptotic stability of boundary equilibria for varioussets of parameters by applying the theory of Lyapunov functions (e.g. LaSalle 1976, inparticular, Theorem 6.4 and Corollary 6.5). By global asymptotic stability of an equilib-rium we mean that every trajectory, such that initially all alleles are present, convergesto this equilibrium. By Lemma 8.1 there is at most one asymptotically stable boundaryequilibrium for any given set of parameters. Hence, convergence of all trajectories tothe boundary is sufficient for demonstrating global stability. Because global convergenceto the boundary precludes the existence of an internal equilibrium, it yields necessaryconditions for a stable DMI. Proposition 8.2.
A DMI can exist only if each of the following conditions is satisfied: α > and γ > β , (8.8) m < α − β + r , (8.9) m < max { α − β, α, ( γ − β ) } . (8.10)85ondition (8.8) means that w ( Ab ) > max { w ( ab ) , w ( aB ) } , i.e., the island type hashigher fitness than its one-step mutational neighbors.The proof of (8.8) uses the Lyapunov functions Y = x + x x + x = 1 − qp (8.11)and X = x + x x + x = 1 − q − p (8.12)to show that convergence to S A or M occurs if (8.8) is violated.Condition (8.9) follows because, if (8.8) is satisfied,˙ x = x [ β ( x + x ) + γx − α ( x + x ) − rx ] + m ( x + x + x ) + rx x ≥ x [( m + β ) x + ( m + β − α − r ) x + ( m − α + γ ) x ] (8.13)holds if m > max {− β, α − γ, α − β + r } = α − β + r . (8.14)Therefore, global convergence to M occurs if (8.14) holds. In particular, M is globallyasymptotically stable for every m if r < β − α .Condition (8.10) follows in a similar way. Lemma 8.3.
Every trajectory eventually enters the region D ≤ and remains there.Convergence to D = 0 occurs if and only if at least one allele is eventually lost. Thus,every internal equilibrium satisfies D < . To prove this lemma, we define Z = x x x x , (8.15)where x > x > Z = 1 if and only if D = 0, and Z <
D >
0. Then˙ Z = x x x ( m + γx ) + rD ( x x x + x x x + x x x + x x x ) . (8.16)We observe that ˙ Z ≥ D ≥
0. In addition, it follows immediately that˙
Z > rD > x + x >
0. If x + x = 0 and x x >
0, then ˙ x + ˙ x < r + m >
0. Hence, all trajectories leave
D > r >
0. If rD = 0, then ˙ Z = 0 only if86 = 0 or if m = 0 and γx = 0. Thus, our result follows by investigating (i) the dynamicson x = 0 if r = 0, (ii) the dynamics on x = 0 if m = r = 0, and (iii) the case m = γ = 0.We leave the simple first two cases to the reader. The third case is also not difficult andfollows immediately from Section 3.4.1 in B¨urger and Akerman (2011).From here on, we assume α > γ > β and r > β − α , (8.17)because we have proved that internal equilibria can exist only if (8.17) is satisfied. Wenote that (8.17) holds if and only if M (island haplotype fixed) is linearly stable in theabsence of migration.Our next aim is the derivation of a cubic equation from which the coordinate p of aninternal equilibrium ( p, q, D ) can be obtained. Given p , the coordinates q and D can becomputed from relatively simple explicit formulas.By solving ˙ p = 0, we find that, for given p and q , and if p (cid:54) = 1 − β/γ , the value of LDat equilibrium is D = D ( p, q ) = p m + (1 − p )( γq − α ) β − γ + γp . (8.18)Substituting this into (8.3b), assuming β (cid:54) = 0, and solving ˙ q = 0 for q , we obtain q , ( p ) = 12 (cid:20)(cid:18) − mβ (cid:19) ± (cid:112) Q (cid:21) , (8.19)where Q = (cid:18) mβ (cid:19) − αmpβ ( γ − β ) − α ( γ − α ) β ( γ − β ) p (1 − p ) (8.20)needs to be nonnegative to yield an admissible equilibrium. Finally, we substitute q = q ( p ) and D = D ( p, q ( p )) into (8.3c) and obtain that an equilibrium value p must solvethe equation ( γ − β ) A ( p ) − (cid:112) QB ( p ) = 0 , (8.21)where A ( p ) = ( γ − β ) (cid:8) β [( α − r )( γ − α ) + γ ( α − β )] + [ β ( γ + β − α ) + r (2 β − γ )] m + βm (cid:9) + (cid:8) β [2 α (2 γ − β )( γ − α ) − βγ ( γ − β )] − (2 γ − β )(2 α − γ ) r + [2 αβ ( γ − β ) − βγ ( γ − β ) + γ (2 γ − β ) r ] m (cid:9) p − [2 αβ ( γ − α )( γ − β ) + βγ ( γ − α ) r + γ rm ] p , (8.22a) B ( p ) = ( γ − β )[ − αβ + γ ( β + r ) + βm ]+ [2 αβ ( γ − β ) − βγ ( γ − β ) + r ( β − γ )] p + γ rp . (8.22b)87f we substitute q = q ( p ) and D = D ( p, q ( p )) into (8.3c), we obtain( γ − β ) A ( p ) + (cid:112) QB ( p ) = 0 , (8.23)instead of (8.21).By simple additional considerations, it is shown that solutions p of (8.21) or (8.23)satisfy 0 = ( γ − β ) A ( p ) − QB ( p ) = 4 βγ − β ( β − γ + γp ) [ m + ( γ − α )(1 − p )] P ( p ) , (8.24)where P ( p ) = ( γ − β )( m + r + β − α )[ αβ ( α + β − γ − r ) − mr ( γ − β )]+ (cid:8) αβ ( α − β )( γ − β )( α + β − γ ) + α ( γ − β )[3 β ( γ − α ) + m (4 β − γ )] r + [ α ( γ + βγ − β ) + γ ( γ − β ) m ] r (cid:9) p − αr [ β ( γ − β )( γ − α ) + γ r ] p + αγ r p . (8.25)Because p = 1 − β/γ never gives an equilibrium of (8.3) and p = 1 − m/ ( γ − α ) can giverise only to a single-locus polymorphism, any internal equilibrium value p must satisfy P ( p ) = 0.We can summarize these findings as follows. Theorem 8.4. (a) The haploid dynamics (8.3) can have at most three internal equilibria,and the coordinate ˆ p of an internal equilibrium (ˆ p, ˆ q, ˆ D ) is a zero of the polynomial P givenby (8.25) .(b) For given ˆ p with P (ˆ p ) = 0 , only one of q (ˆ p ) or q (ˆ p ) defined in (8.19) can yieldan equilibrium value ˆ q . ˆ D is calculated from ˆ p and ˆ q by (8.18) . This procedure yields aninternal equilibrium if and only if < ˆ p < , < ˆ q < , and − min[ˆ p ˆ q, (1 − ˆ p )(1 − ˆ q )] < ˆ D < . (8.26) Remark 8.5.
In the absence of epistasis, i.e., if γ = 0, there can be at most two in-ternal equilibria. Their coordinates are obtained from a quadratic equation in p . Theadmissibility conditions are given by simple formulas (B¨urger and Akerman 2011).If r = 0, the equilibrium R mentioned in Section 8.2 is obtained from the unique zeroof P ( p ).For a number of important special or limiting cases, explicit expressions or approxi-mations can be obtained for the internal equilibria.88 .4 Bifurcations of two internal equilibria A bifurcation of two internal equilibria can occur if and only if P ( p ∗ ) = 0, where p ∗ ∈ (0 , P , i.e., P (cid:48) ( p ∗ ) = 0. There are at most two such critical points, andthey are given by p ∗ , = 13 αγ r (cid:110) α [ β ( γ − α )( γ − β ) + γ r ] ± √ R (cid:111) , (8.27a)where R = α (cid:8) − β ( γ − β ) (cid:2) αβ ( γ − α )( γ − β ) + γ (3 α + β ) − γ (3 α + β ) (cid:3) − βγ ( γ − β )( γ − α ) r + γ (3 β − βγ + γ ) r (cid:9) − mrαγ ( γ − β )[4 αβ + γ ( r − α )] . (8.27b)Solving either P ( p ∗ ) = 0 or P ( p ∗ ) = 0 for m , we obtain after some straightforward ma-nipulations that the critical value m ∗ must be a solution of the following quartic equation:[ αβ (2 β − γ )( α + β − γ + r ) + γ ( γ − β ) rm ] [ − ψ + 2 ψ rm + 27 αγ ( γ − β ) r m ] = 0 , (8.28a)where ψ = α ( γ − β )( α − β + r ) [4 αβ ( γ − α )( γ − β ) + γ r ] , (8.28b) ψ = 2 α ( γ − β ) (9 βγ − αβ − αγ ) + 3 αγ ( γ − β )(3 βγ − αβ − αγ ) r + 3 αγ ( γ − β ) r + 2 γ r . (8.28c)The zero m = m arising from the first (linear) factor in (8.28a) does not give a validbifurcation point for internal equilibria.The second (quadratic) factor in (8.28a) provides two potential solutions. However,because ψ ≥
0, one is negative. Therefore, the critical value we are looking for is givenby m ∗ = 127 αγ ( γ − β ) r (cid:18) − ψ + (cid:113) ψ + 27 αγ ( γ − β ) ψ (cid:19) . (8.29)At this value, two equilibria with non-zero allele frequencies collide and annihilate eachother. Thus, m ∗ is the critical value at which a saddle-node bifurcation occurs. Thisgives an admissible bifurcation if both equilibria are internal (hence admissible) for either m < m ∗ or m > m ∗ .If p ∗ = p ∗ , i.e., if R = 0, a pitchfork bifurcation could occur at m ∗ . As a function of α , β , and γ , the condition p ∗ = p ∗ can be satisfied at m ∗ only for three different values89f 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 > m ∗ may be found in the Online Supplementof Bank et al. (2012). We assume m = 0. From Section 8.3, we obtain the following properties of internalequilibria ( p, q, D ). The LD is given by D = D ( p, q ) = p (1 − p ) γq − αβ − γ + γp , (8.30)where p (cid:54) = 1 − β/γ ; cf. (8.18). For admissibility, (8.26) needs to be satisfied. For given p and if β (cid:54) = 0, the coordinate q of an internal equilibrium can assume only one of thefollowing forms: q , ( p ) = 12 (cid:32) ± (cid:115) − α ( γ − α ) p (1 − p ) β ( γ − β ) (cid:33) . (8.31)By Theorem 8.4, for given p , at most one of q = q ( p ) or q = q ( p ) can give rise to anequilibrium.The following theorem characterizes the equilibrium and stability structure. Theorem 8.6.
Suppose (8.17) and m = 0 .(a) The haploid dynamics (8.3) admits at most one internal equilibrium.(b) Depending on the parameters, the internal equilibrium is given by either ( p, q ( p ) , D ( p, q ( p ))) or ( p, q ( p ) , D ( p, q ( p ))) , where p is one of the two zeros of P ( p ) = γ r p − r [2 β ( γ − β )( γ − α ) + rγ ] p + β ( γ − β )( α − β − r )( α + β − γ − r ) , (8.32) and q i ( p ) and D ( p, q i ( p )) are given by (8.31) and (8.30) , respectively.(c) An internal equilibrium exists if and only if both M and M are asymptoticallystable. This is the case if and only if γ > α and β > and r > α − β . (8.33) (d) The internal equilibrium is unstable whenever it exists.(e) If (8.33) does not hold, then M is globally asymptotically stable. Remark 8.7 (Hofbauer’s index theorem) . Theorem 2 in Hofbauer (1990) states the fol-lowing. For every dissipative semiflow on R n + such that all fixed points are regular, thesum of the indices of all saturated equilibria equals +1.In our model, the index of an equilibrium is ( − m , where m is the number of negativeeigenvalues (they are always real). An internal equilibrium is always saturated. If itis asymptotically stable, it has index 1. Equilibria on the boundary of the simplex aresaturated if and only if they are externally stable. This is the case if and only if no gametethat is missing at the equilibrium can invade. Because S is attracting within R , theindex 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 technicallemmas that provide properties of the polynomial P ( p ) and resulting admissibility condi-tions for ( p, q, D ).The proof of (c) and (d) is the mathematically most interesting part. In view of theabove, it remains to prove that the internal equilibrium exists if (8.33) holds and that itis unstable. This follows readily from the index theorem of Hofbauer (1990) in Remark8.7. In our model, the only boundary equilibria are the four monomorphic states. M and M are never saturated because they are unstable within S . M and M are saturatedif and only if they are asymptotically stable within S . Then, ind( M ) = ind( M ) = 1.Hence there must exist an internal equilibrium with index -1. Such an equilibrium cannotbe stable. Because M and M are both asymptotically stable if and only if (8.33) holds,statements (c) and (d) are proved.(e) follows with the help of several simple Lyapunov functions. Theorem 7.6 and Remark 7.7 in conjunction with Theorem 8.6 and Lemma 8.1 yield theequilibrium configuration for weak migration.91 heorem 8.8. If m > is sufficiently small, the following equilibrium configurationsoccur.(a) If (8.33) holds, there exists one unstable internal equilibrium and one asymptoti-cally stable internal equilibrium (the perturbation of M ). In addition, the monomorphicequilibrium M is asymptotically stable. Neither S A nor S B is admissible.(b) Otherwise, i.e., if γ < α or β < or r < α − β , the perturbation of the equilibrium M is globally asymptotically stable (at least if γ is small). The equilibrium M 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 itscoordinates are (cid:18) − m ( α + r ) α ( α − β + r ) , m ( γ − β + r )( γ − β )( α − β + r ) , − mα − β + r (cid:19) . (8.34) Now we are in the position to prove our main results about the equilibrium and bifurcationstructure. We continue to assume (8.17), so that a DMI can occur. Throughout, we alwaysconsider bifurcations as a function of (increasing) m .We define m A = ( α − γ )( γ − β ) α (cid:18) α + β − γr (cid:19) , (8.35a) m B = − βαγ − β (cid:18) γ − β − αr (cid:19) , (8.35b) m = α − β − r (8.35c)and note that m A , m B , and m are the critical values of m above which S A , S B , and M ,respectively, are asymptotically stable provided they are admissible. If we set m − max = m A if γ − α < < γ − β < α and r > γ − β , (8.36a) m B if β < < α < γ − β and r > α , (8.36b) m if r ≤ min[ α, α − β, γ − β ], (8.36c)then, by Lemma 8.1, all boundary equilibria are repelling if and only if m < m − max . Asthe following theorem shows, if m < m − max , an asymptotically stable internal equilibriumexists which, presumably, is globally attracting. Hence, a DMI will evolve from everyinitial condition. Of course, m − max can be zero (if 0 < α < β < γ and r ≤ β − α ). We alsorecall the definition of m ∗ from (8.29). 92 heorem 8.9. The following three types of bifurcation patterns can occur:Type 1. • If < m < m ∗ , there exist two internal equilibria; one is asymptotically stable ( I DMI ) , the other ( I ) is unstable. The monomorphic equilibrium M is asymptoti-cally stable. • At m = m ∗ , the two internal equilibria collide and annihilate each other by a saddle-node bifurcation. • If m > m ∗ , M is the only equilibrium; it is asymptotically stable and, presumably,globally stable.Type 2. There exists a critical migration rate ˜ m satisfying < ˜ m < m ∗ such that: • If < m < ˜ m , there is a unique internal equilibrium ( I DMI ) . It is asymptoticallystable and, presumably, globally stable. • At m = ˜ m , an unstable equilibrium ( I ) enters the state space by an exchange-of-stability bifurcation with a boundary equilibrium. • If ˜ m < m < m ∗ , there are two internal equilibria, one asymptotically stable ( I DMI ) ,the other unstable ( I ) , and 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 asymptotically stable and, presumably, globallystable.Type 3. • If < m < m − max , a unique internal equilibrium ( I DMI ) exists. It is asymptoticallystable and, presumably, globally stable. • At m = m − max , I DMI leaves the state space through a boundary equilibrium by anexchange-of-stability bifurcation. • If m > m − max , a boundary equilibrium is asymptotically stable and, presumably,globally stable.Proof. Theorem 8.8 provides all equilibrium configurations for small m . Lemma 8.1 pro-vides control over the boundary equilibria. As m increases, they can vanish but notemerge. They can also become asymptotically stable as m increases. For sufficiently large m , there is always a globally asymptotically stable boundary equilibrium. By Theorem93 ) e) M md) mS A or S B I DMI M I DMI S A or S B b) mc)S A or S B M M I I I DMI S A or S B I DMI m M I DMI m Figure 8.1:
Bifurcation patterns in the DMI model according to Theorem 8.9. Panel a showsthe bifurcation pattern of Type 1. Panels b and c display the two cases that are of Type 2, andpanels d and e those of Type 3. S through one of boundary equilibria, when an exchangeof stability occurs. A bifurcation involving the two internal equilibria can occur at mostat one value of m , namely at m ∗ (8.29). An exchange-of-stability bifurcation can occuronly at the values m A , m B , or m . If it occurs, the respective boundary equilibrium isasymptotically stable for every larger m for which it is admissible. By the index theoremof Hofbauer (1990), Remark 8.7, the sum of the indices of all saturated equilibria is 1.If Case 1 of Theorem 8.8 applies, then M is asymptotically stable for every m >
0, andit is the only boundary equilibrium. Hence, its index is 1. The index of the stable internalequilibrium 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 thetwo internal equilibria can occur and because for large m , M is globally asymptoticallystable, the bifurcation must be of saddle-node type in which the equilibria collide andannihilate each other (but do not emerge). In principle, the internal equilibria couldalso leave S through a boundary equilibrium (in this case, it must be M ). However,by the index theorem, they can do so only simultaneously. This occurs if and onlyif m ∗ = m , which is a non-generic degenerate case. Because the sum of indices ofthe internal equilibria must be zero, no equilibrium can enter the state space. Theseconsiderations settle the bifurcation pattern of Type 1.If Case 2 of Theorem 8.8 applies, then, for small m , the boundary equilibria are un-stable, hence not saturated, and do not contribute to the sum of indices. Since then theindices of internal equilibria must sum up to 1, the only possible bifurcation that doesnot entail the stability of a boundary equilibrium would be a pitch-fork bifurcation of theinternal equilibrium which, by Section 8.4, does not occur. Indeed, because m ∗ is the onlyvalue at which a bifurcation among internal equilibria can occur, and because for large m a boundary equilibrium is globally asymptotically stable, the three equilibria emerging bya pitchfork bifurcation would have to leave the state space through boundary equilibria.This, however, cannot occur, as follows easily from the results about linear stability inSection 8.2. Thus, any further bifurcations involve a boundary equilibrium. There aretwo possibilities.(i) An equilibrium enters S at some value ˜ m (which can only be one of m A , m B , or m )through one of the unstable boundary equilibria by an exchange-of-stability bifurcation.If m > ˜ m , there is one asymptotically stable boundary equilibrium, an unstable internalequilibrium (the one that entered S ), and one asymptotically stable internal equilibrium.95ow a reasoning analogous to that applied above to Case 1 of Theorem 8.8 establishesthe bifurcation pattern of Type 2.(ii) The internal equilibrium leaves S by exchange of stability through one of theboundary equilibria at m − max . This becomes asymptotically stable then and, presumably,globally stable. At larger values of m no equilibrium can enter S through one of the(other) unstable boundary equilibria, because this would either lead to two simultaneouslystable boundary equilibria, which is impossible (Lemma 8.1), or this had to occur at thesame value at which the hitherto stable boundary equilibrium merges with M and leavesthe state space. This, too, is impossible because the sum of the indices of the new stableboundary equilibrium and the new unstable internal equilibrium would be zero. Thus, wehave established the bifurcation pattern of Type 3 and excluded all other possibilities.Our final goal is to assign the respective parameter combinations to the three types ofbifurcation patterns determined above. To this aim we define four selection scenarios thatreflect different biological situations. We say that the fitness landscape is of slope type ifat least one of the one-mutant neighbors ( ab or AB ) has fitness intermediate between thecontinental ( aB ) and the island type ( Ab ). Otherwise, there is a double-peak landscape. Selection scenario 1 : 0 < α < β < γ . It represents the parameter regime in which “se-lection against hybrids” (of continental and island haplotypes) is driving DMI evolution.The fitness landscape has two peaks with the higher at the continental type. The geneticincompatibility is strong ( γ large). Selection scenario 2 : 0 < β < α < γ . It represents an intermediate parameterregime in which the two mechanisms of “selection against hybrids” and “selection againstimmigrants” (below) interfere. This is also a double-peak landscape, but with maximumat the island type. The incompatibility is strong.
Selection scenario 3 : β < γ < min { α, α + β } . It represents part of the parame-ter regime in which “selection against immigrants” drives DMI evolution. The fitnesslandscape is of slope type and the incompatibility is weak. β may be positive or negative. Selection scenario 4 : β < < α < γ − β . It represents part of the parameter regimein which “selection against immigrants” drives DMI evolution. The fitness landscape isof slope type and there is strong local adaptation ( β < < α ).96o formulate the main result, we need the following critical recombination rates: r A = ( γ − α ) 3( γ − β ) − α γ − α , (8.37a) r B = β α + β − γβ + γ , (8.37b) r = 3 α ( γ − β ) − (cid:112) α ( γ − β )(4 βγ + 5 αγ − αβ )2 γ . (8.37c) Theorem 8.10.
1. Bifurcation patterns of Type 1 occur inSelection scenario 1;Selection scenario 2 if and only if r ≥ α − β .2. Bifurcation patterns of Type 2 occur inSelection scenario 2 if and only if r < r < α − β ;Selection scenario 3 if and only if one of the following holds:(a) r ∗ < r ≤ γ − β ,(b) r > max[ γ − β, r A ] and γ > α ,(c) γ − β < r < ∞ and γ = α > β ,(d) γ − β < r < r A and γ < α .Selection scenario 4 if and only if one of the following holds:(a) r ∗ < r ≤ α ,(b) r > max[ α, r B ] and γ > − β ,(c) α < r < ∞ and γ = − β < α + β ,(d) α < r < r B and γ < − β .3. Bifurcation patterns of Type 3 occur inSelection scenario 2 if and only if r ≤ r ;Selection scenario 3 if and only if one of the following holds:(a) r ≤ min[ γ − β, r ∗ ] ,(b) γ − β < r ≤ r A and γ > α ,(c) γ − β < r < ∞ and γ = α < β .(d) r ≥ max[ γ − β, r A ] and γ < α .Selection scenario 4 if and only if one of the following holds:(a) r ≤ min[ α, r ∗ ] ,(b) α < r ≤ r B and γ > − β ,(c) α < r < ∞ and γ = − β ≥ α + β ,(d) r ≥ max[ α, r B ] and γ < − β . γ ( α + β ) < β (2 α + β ) and γ < − β hold, then a bifurcationpattern of Type 3 occurs for every r >
0. These two conditions are satisfied whenever − β > max[2 α, γ ]. Therefore, in the strong local adaptation scenario a globally asymptot-ically stable DMI occurs whenever the selection intensity on the two loci differs by morethan a factor of two. A bistable equilibrium pattern can occur in this scenario only if theselection strength on both loci is sufficiently similar and the recombination rate is aboutas strong as the selection intensity.In summary, two mechanisms can drive the evolution of parapatric DMIs. In a hetero-geneous environment, a DMI can emerge as a by-product of selection against maladpatedimmigrants. In a homogeneous environment, selection against unfit hybrids can main-tain a DMI. No DMI can be maintained if the migration rate exceeds one of the boundsgiven in Proposition 8.2. Therefore, it is the adaptive advantage of single substitutionsrather than the strength of the incompatibility that is the most important factor for theevolution of a DMI with gene flow. In particular, neutral DMIs can not persist in thepresence of gene flow. Interestingly, selection against immigrants is most effective if theincompatibility loci are tightly linked, whereas selection against hybrids is most effectiveif they are loosely linked. Therefore, opposite predictions result concerning the geneticarchitecture that maximizes the rate of gene flow a DMI can sustain. Akin E. 1979. The Geometry of Population Genetics. Lect. Notes Biomath. 31. BerlinHeidelberg New York: Springer.Akin E. 1982. Cycling in simple genetic systems. J. Math. Biol. 13, 305-324.Akin, E. 1993. The General Topology of Dynamical Systems. Providence, R.I.: Amer.98ath. Soc.Bank, C., B¨urger, R., Hermisson, J. 2012. The limits to parapatric speciation: Dobzhansky-Muller incompatibilities in a continent-island model. Genetics 191, 845863.Barton, N.H. 1999. Clines in polygenic traits. Genetical Research 74, 223-236.Barton, N.H. 2010. What role does natural selection play in speciation? Phil. Trans. R.Soc. B 365, 1825-1840.Barton, N.H., Turelli, M. 2011. Spatial waves of advance with bistable dynamics: cyto-plasmic and genetic analogues of Allee effects. Amer. Natur. 178, No. 3, pp. E48-E75.Baum, L.E., Eagon, J.A. 1967. An inequality with applications to statistical estimationfor probability functions of Markov processes and to a model for ecology. Bull. Amer.Math. Soc. 73, 360-363.Berman, A., Plemmons, R.J. 1994. Nonnegative Matrices in the Mathematical Sciences.Philadelphia: SIAM.Bulmer, M.G. 1972. Multiple niche polymorphism. Amer. Natur. 106, 254-257.B¨urger, R. 2000. The Mathematical Theory of Selection, Recombination, and Mutation.Wiley, Chichester.B¨urger, R. 2009a. Multilocus selection in subdivided populations I. Convergence proper-ties for weak or strong migration. J. Math. Biol. 58, 939-978.B¨urger, R. 2009b. Multilocus selection in subdivided populations II. Maintenance ofpolymorphism and weak or strong migration. J. Math. Biol. 58, 979-997.B¨urger, R. 2009c. Polymorphism in the two-locus Levene model with nonepistatic direc-tional selection. Theor. Popul. Biol. 76, 214-228.B¨urger, R. 2010. Evolution and polymorphism in the multilocus Levene model with noor weak epistasis. Theor. Popul. Biol. 78, 123-138.B¨urger, R. 2011. Some mathematical models in evolutionary genetics. In: The Mathe-matics of Darwins Legacy (FACC Chalub & JF Rodrigues, eds), pp. 67-89. Birkh¨auser,Basel.B¨urger, R., Akerman, A. 2011. The effects of linkage and gene flow on local adaptation:A two-locus continent-island model. Theor. Popul. Biol. 80, 272-288.Cannings, C. 1971. Natural selection at a multiallelic autosomal locus with multipleniches. J. Genetics 60, 255-259.Christiansen, F.B. 1974. Sufficient conditions for protected polymorphism in a subdividedpopulation. Amer. Naturalist 108, 157-166.Christiansen, F.B. 1975. Hard and soft selection in a subdivided population. Amer.Natur. 109, 11-16. 99hristiansen, F.B. 1999. Population Genetics of Multiple Loci. Wiley, Chichester.Conley, C. 1978. Isolated invariant sets and the Morse index. NSF CBMS Lecture Notes38. Providence, R.I.: Amer. Math. Soc.Deakin, M.A.B. 1966. Sufficient conditions for genetic polymorphism. Amer. Natur. 100,690-692.Deakin, M.A.B. 1972. Corrigendum to genetic polymorphism in a subdivided population.Australian J. Biol. Sci. 25, 213-214.Dempster, E.R. 1955. Maintenance of genetic heterogeneity. Cold Spring Harbor Symp.Quant. Biol. 20, 25-32.Ewens, W. J. 1969. Mean fitness increases when fitnesses are additive. Nature 221, 1076.Ewens, W.J. 2004. Mathematical Population Genetics. 2nd edition. Springer, New York.Ewens, W.J. 2011. What changes has mathematics made to the Darwinian theory? In:The Mathematics of Darwins Legacy (FACC Chalub & JF Rodrigues, eds), pp. 7-26.Birkh¨auser, Basel.Eyland, E.A. 1971. Moran’s island model. Genetics 69, 399-403.Feller, W. 1968. An Introduction to Probability Theory and Its Applications, vol. I, thirdedn. Wiley, New York.Fisher, R.A. 1937. The wave of advance of advantageous genes. Ann. Eugenics 7, 355-369.Friedland, S., Karlin, S. 1975. Some inequalities for the spectral radius of nonnegativematrices and applications. Duke Math. J. 42, 459-490.Geiringer, H. 1944. On the probability theory of linkage in Mendelian heredity. Ann.Math. Stat. 15, 25-57.Hadeler, K.P., Glas, D., 1983. Quasimonotone systems and convergence to equilibrium ina population genetic model. J. Math. Anal. Appl. 95, 297-303.Haldane, J.B.S. 1930. A mathematical theory of natural and artificial selection. Part VI.Isolation. Proc. Camb. Phil. Soc. 28, 224-248.Haldane, J.B.S. 1932. The Causes of Evolution. London: Longmans, Green. (Reprintedwith a new introduction and afterword by E.G. Leigh, Jr., by Princeton UniversityPress, 1992.)Haldane, J.B.S. 1948. The theory of a cline. J. Genetics 48, 277-284.Hardy, G.H. 1908. Mendelian proportions in a mixed population. Science 28, 49-50.Hastings A. 1981a. Simultaneous stability of D = 0 and D (cid:54) = 0 for multiplicative viabili-ties at two loci: an analytical study. J. Theor. Biol. 89, 69-81.Hastings A. 1981b. Stable cycling in discrete-time genetic models. Proc. Natl. Acad. Sci.USA 78, 7224-7225. 100irsch, M. W., 1982. Systems of differential equations which are competitive or coopera-tive. I: Limit sets. SIAM J. Math. Anal. 13, 167-179.Hofbauer J., Iooss G. 1984. A Hopf bifurcation theorem of difference equations approxi-mating a differential equation. Monatsh. Math. 98, 99-113.Hofbauer, J., Sigmund, K. 1988. The Theory of Evolution and Dynamical Systems.Cambridge: University Press.Hofbauer, J., Sigmund, K. 1998. Evolutionary Games and Population Dynamics. Cam-bridge: University Press.Karlin, S. 1977. Gene frequency patterns in the Levene subdivided population model.Theor. Popul. Biol. 11, 356-385.Karlin, S. 1982. Classification of selection-migration structures and conditions for a pro-tected polymorphism. Evol. Biol. 14, 61-204.Karlin, S., Campbell, R.B. 1980. Selection-migration regimes characterized by a globallystable equilibrium. Genetics 94, 1065-1084.Karlin S., Feldman, M.W. 1978. Simultaneous stability of D = 0 and D (cid:54) = 0 for multi-plicative viabilities at two loci. Genetics 90, 813-825.Karlin, S., McGregor, J. 1972a. Application of method of small parameters to multi-nichepopulation genetics models. Theor. Popul. Biol. 3, 186-208.Karlin, S., McGregor, J. 1972b. Polymorphism for genetic and ecological systems withweak coupling. Theor. Popul. Biol. 3, 210-238.Kingman, J.F.C. 1961. An inequality in partial averages. Quart. J. Math. 12, 78-80.Kolmogoroff, A., Pretrovsky, I., Piscounoff, N. 1937. ´Etude de l’´equation de la diffusionavec croissance de la quantite de mati´ere et son application `a un probl`eme biologique.(French) Bull. Univ. Etat Moscou, Ser. Int., Sect. A, Math. et Mecan. 1, Fasc. 6, 1-25.LaSalle, J.P. 1976. The Stability of Dynamical Systems. Regional Conf. Ser. Appl. Math.25. Philadelphia: SIAM.Levene, H. 1953. Genetic equilibrium when more than one ecological niche is available.Amer. Natur. 87, 331-333.Lessard S. 1997. Fisher’s fundamental theorem of natural selection revisited. Theor. Pop.Biol. 52, 119-136.Lewontin, R.C., Kojima K.-I. 1960. The evolutionary dynamics of complex polymor-phisms. Evolution 14, 458-472.Li, C.C. 1955. The stability of an equilibrium and the average fitness of a population.Amer. Natur. 89, 281-295.Li, C.C. 1967. Fundamental theorem of natural selection. Nature 214, 505-506.101ou, Y., Nagylaki, T., 2002. A semilinear parabolic system for migration and selection inpopulation genetics. J. Differential Equations 181, 388418.Lou, Y., Nagylaki, T., Ni, W.-M. 2013. An introduction to migration-selection PDEmodels. Discrete Continuous Dynamical Systems, Series A, in press.Lyubich, Yu.I. 1971. Basic concepts and theorems of evolutionary genetics of free popu-lations. Russ. Math. Surv. 26, 51-123.Lyubich, Yu.I. 1992. Mathematical Structures in Population Genetics. Berlin HeidelbergNew York: Springer.Maynard Smith, J. 1970. Genetic polymorphism in a varied environment. Am. Nat. 104,487-490.Nagylaki, T. 1977. Selection in One- and Two-Locus Systems. Lect. Notes Biomath. 15.Berlin, Heidelberg, New York: Springer.Nagylaki, T. 1989. The diffusion model for migration and selection. Pp. 5575 in: Hastings,A. (Ed.), Some Mathematical Questions in Biology. Lecture Notes on Mathematics inthe Life Sciences, vol. 20. American Mathematical Society, Providence, RI.Nagylaki, T. 1992. Introduction to Theoretical Population Genetics. Berlin: Springer.Nagylaki, T. 1993. The evolution of multilocus systems under weak selection. Genetics134, 627-647.Nagylaki, T. 2009a. Polymorphism in multiallelic migration-selection models with domi-nance. Theor. Popul. Biol. 75,239-259.Nagylaki, T. 2009b. Evolution under the multilocus Levene model without epistasis.Theor. Popul. Biol. 76, 197-213.Nagylaki, T., Hofbauer, J., Brunovsk´y, P., 1999. Convergence of multilocus systems underweak epistasis or weak selection. J. Math. Biol. 38, 103-133.Nagylaki, T., Lou, Y. 2001. Patterns of multiallelic poylmorphism maintained by migra-tion and selection. Theor. Popul. Biol. 59, 297-33.Nagylaki, T., Lou, Y. 2006a. Multiallelic selection polymorphism. Theor. Popul. Biol.69, 217-229.Nagylaki, T., Lou, Y. 2006b. Evolution under the multiallelic Levene model. Theor.Popul. Biol. 70, 401-411.Nagylaki, T., Lou, Y. 2007. Evolution under multiallelic migration-selection models.Theor. Popul. Biol. 72, 21-40.Nagylaki, T., Lou, Y. 2008. The dynamics of migration-selection models. In: Friedman,A. (ed) Tutorials in Mathematical Biosciences IV. Lect. Notes Math. 1922, pp. 119 -172. Berlin Heidelberg New York: Springer.102ovak, S. 2011. The number of equilibria in the diallelic Levene model with multipledemes. Theor. Popul. Biol. 79, 97101.Peischl, S. 2010. Dominance and the maintenance of polymorphism in multiallelic migration-selection models with two demes. Theor. Popul. Biol. 78, 12-25.Price, G.R. 1970. Selection and covariance. Nature 227, 520-521.Prout, T. 1968. Sufficient conditions for multiple niche polymorphism. Amer. Natur. 102,493-496.Seneta, E. 1981. Non-negative Matrices, 2nd ed. London: Allen and Unwin.Spichtig, M., Kawecki, T.J. 2004. The maintenance (or not) of polygenic variation by softselection in a heterogeneous environment. Amer. Natur. 164, 70-84.Star, B., Stoffels, R.J., Spencer, H.G. 2007a. Single-locus polymorphism in a heteroge-neous two-deme model. Genetics 176, 1625-1633.Star, B., Stoffels, R.J., Spencer, H.G. 2007b. Evolution of fitnesses and allele frequenciesin a population with spatially heterogeneous selection pressures. Genetics 177, 1743-1751.Strobeck, C. 1979. Haploid selection with n alleles in mm