aa r X i v : . [ m a t h . P R ] O c t The Annals of Applied Probability (cid:13)
Institute of Mathematical Statistics, 2009
SPATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL
By N. Lanchier and C. Neuhauser
University of Minnesota
We introduce a spatially explicit model for the competition be-tween type a and type b alleles. Each vertex of the d -dimensionalinteger lattice is occupied by a diploid individual, which is in one ofthree possible states or genotypes: aa , ab or bb . We are interested inthe long-term behavior of the gene frequencies when Mendel’s lawof segregation does not hold. This results in a voter type model de-pending on four parameters; each of these parameters measures thestrength of competition between genes during meiosis. We prove thatwith or without a spatial structure, type a and type b alleles coexistat equilibrium when homozygotes are poor competitors. The inclu-sion of a spatial structure, however, reduces the parameter regionwhere coexistence occurs.
1. Introduction.
Diploid organisms carry two sets of chromosomes, whichare segregated during meiosis to produce haploid gametes. Mendel’s lawof segregation states that pairs of homologous chromosomes segregate ran-domly during meiosis and are distributed to gametes. This implies that agamete of an individual who is heterozygous at a locus will receive one orthe other allele with equal probability. This implies, for instance, that 25%of offspring of two individuals that are both heterozygous for the same pairof alleles will be homozygous for one allele, 25% will be homozygous for theother allele and 50% will be heterozygous.This law does not always hold. Some organisms have chromosomal el-ements that show a higher than expected rate of transmission. A widely-studied example of meiotic drive can be found in the house mouse (
Musmusculus ). The house mouse may have a variant form of chromosome 17,called t -haplotype, which is transmitted to more than 90% of its offspring.Such systems that show a higher than 50% representation in their offspringare called meiotic drive systems. Received October 2007; revised February 2009.
AMS 2000 subject classifications.
Primary 60K35, 82C22; secondary 92D25.
Key words and phrases.
Voter model, annihilating branching process, non-Mendeliansegregation.
This is an electronic reprint of the original article published by theInstitute of Mathematical Statistics in
The Annals of Applied Probability ,2009, Vol. 19, No. 5, 1880–1920. This reprint differs from the original inpagination and typographic detail. 1
N. LANCHIER AND C. NEUHAUSER
Simple mathematical models readily show that meiotic drive by itselfresults in fixation. Yet natural populations show a low frequency of the t -haplotype [10–15%, see Durand et al. (1997) and references therein]. Twofactors are contributing to the lower than expected frequency: male offspringwho are homozygous for the t -haplotype are sterile (which by itself precludesfixation) and most t -haplotypes carry additional recessive lethal mutants. AsBruck (1957) showed, however, these two factors do not sufficiently explainthe observed low frequencies. Meiotic drive has been observed in a number ofother systems [see articles in Lyttle and Perkins (1991)], including the fruitfly Drosophila, mosquitoes, the fungus Neurospora, humans and maize. Inorder to detect meiotic drive in natural populations, a polymorphism mustbe present.Lewontin and Dunn (1960) used computer simulations to suggest that ran-dom genetic drift due to small population size lowers Bruck’s predicted fre-quency of the t -haplotype further. Other simulations with varying assump-tions were conducted over the years [Levin, Petras and Rasmussen (1969);Petras and Topping (1983); Nunney and Baker (1993)], which are reviewedin Durand et al. (1997). Durand et al. (1997) concluded that “[N]one of thework surveyed here has fully explained the existence of the low frequency t -polymorphism seen in nature.” They then posit a number of hypothesesto explain the low observed frequency, including low levels of migration in aspatially structured population, selection against the t -haplotype and lowertransmission ratio distortion in wild versus laboratory populations. To testthese hypotheses, they built a spatially structured simulation model. Theyfound that “migration rate is the single most important factor in determin-ing the equilibrium frequency of t -haplotypes.” However, the model onlyyielded observed low frequencies for a very narrow range of parameters, andtherefore essentially failed to provide an explanation for the low frequency of t -haplotypes in the wild. Their extensive simulations of this spatially struc-tured model led this group to suggest that “the basic concept of a stable,low frequency t -polymorphism may be wrong.” Instead, there might be twostable states, one in which the t -haplotype is absent and another where it isat high frequency (as in unstructured populations). Environmental changemay be responsible for populations switching between the two states. Theyconcluded with “[I]f these [transitions] occur frequently, many populationsmoving towards a new equilibrium (either t -haplotype extinction or pan-mictic frequencies) will be observed. The allele frequency measured in suchpopulations would be low but transient.”Spatial simulations are difficult to interpret and we would like to cautionagainst drawing conclusions that are not supported by rigorous mathemati-cal analysis. Rigorous results are typically restricted to specific models, butgeneralities have emerged that allow us to speculate with some confidenceon the behavior of models. While simulations might indicate the absence of a PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL stable polymorphism for realistic parameter values and instead a presence ofmultiple stable states, we propose that a spatially explicit model on meioticdrive will either have a locally stable equilibrium, that is, a nontrivial sta-tionary distribution leading to a stable polymorphism, or fixation will occur,where the fixation states correspond to the multiple stable states of the non-spatial model. However, we suggest that the model will not alternate amongthe fixation states (except in degenerate cases). This assertion is based onthe observation that if a nonspatial model exhibits multiple, locally stableequilibria, it can often be shown that the corresponding spatially explicit,stochastic model will eventually converge to only one of these states and thatall other fixation states are metastable [see, for instance, Neuhauser (1994)].In addition, one- and two-dimensional spatial stochastic models often exhibitclustering. These clusters are not fixed in space but move around, and thusmay give the impression of local transience. The model we now introduceappears to exhibit these commonly observed behaviors as well.Our stochastic model is a continuous time Markov process in which thestate at time t is a function η t : Z d −→ { aa, ab, bb } . Each site of the d -dimensional integer lattice is occupied by an individual with η t ( x ) indicatingthe genotype ( aa , ab or bb ) of the individual present at site x at time t . Eachof the two genes at site x gives birth at varying rates to a gene of its owntype that then replaces one of the 4 d genes picked at random among the 2 d nearest neighbors. More precisely, the evolution at site x is described by thetransition rates aa → ab at rate 2 φ bb X k x − z k =1 { η ( z ) = bb } + φ ba X k x − z k =1 { η ( z ) = ab } ,ab → bb at rate φ bb X k x − z k =1 { η ( z ) = bb } + φ ba X k x − z k =1 { η ( z ) = ab } ,bb → ab at rate 2 φ aa X k x − z k =1 { η ( z ) = aa } + φ ab X k x − z k =1 { η ( z ) = ab } ,ab → aa at rate φ aa X k x − z k =1 { η ( z ) = aa } + φ ab X k x − z k =1 { η ( z ) = ab } , where k x − z k = sup i =1 , ,...,d | x i − z i | . For i, j ∈ { a, b } , the parameter φ ij isthe birth rate of alleles of type i present in sites occupied by individualsof genotype ij . To understand the first two rates, note that an individualof genotype bb produces genes of type b at rate 2 φ bb while an individual ofgenotype ab produces genes of type b at rate φ ba . If the gene (of type b ) issent to an individual of genotype aa , then it will replace one of the type a genes which always results in a transition aa → ab . But if the gene (of type b ) is sent to an individual of genotype ab , then it will replace one of the genes N. LANCHIER AND C. NEUHAUSER chosen at random which results in a transition ab → bb with probability 1 / a and b . The mean-field model.
Before investigating the consequences of the in-clusion of a spatial structure in the form of local interactions, we will look atthe mean-field model [see Durrett and Levin (1994)]. The mean-field modelcan be obtained from the spatial model by assuming that individuals arelocated on the set of sites of a complete graph with N vertices (i.e., eachsite has all other sites as its neighbors which implies that there is no spa-tial structure), by dividing the parameters φ ij by the neighborhood size N and by taking the limit as N tends to infinity. This results in a nonspatialdeterministic model in which the densities of genotypes evolve according toa system of differential equations. More precisely, by letting u aa , u ab and u bb denote the densities of individuals of genotype aa , ab and bb , respec-tively, the mean-field model is described by the following coupled system ofordinary differential equations: du aa dt = (2 φ aa u aa + φ ab u ab ) u ab / − (2 φ bb u bb + φ ba u ab ) u aa , (1) du ab dt = (2 φ aa u aa + φ ab u ab )( u bb − u ab / φ bb u bb + φ ba u ab )( u aa − u ab / ,du bb dt = (2 φ bb u bb + φ ba u ab ) u ab / − (2 φ aa u aa + φ ab u ab ) u bb . (3)Let ψ = ( φ aa − φ ba )( φ bb − φ ab ). When ψ ≤
0, the system has exactly twofixed points, the two trivial equilibria u aa = 1 and u bb = 1. By assumingthat ψ >
0, we find an additional (nontrivial) equilibrium given by φ u aa = ( φ bb − φ ab ) ,φ u bb = ( φ aa − φ ba ) ,φ u ab = 2( φ aa − φ ba )( φ bb − φ ab ) , where φ = ( φ aa − φ ba ) + ( φ bb − φ ab ). In Section 2 we will prove the followingresults:1. The trivial equilibrium u aa = 1 is locally stable if and only if φ aa > φ ba .2. By symmetry, the equilibrium u bb = 1 is locally stable if and only if φ bb >φ ab .3. The interior fixed point ( ψ >
0) is locally stable if and only if φ aa < φ ba and φ bb < φ ab .See Figure 1 for a summary of these results. PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL Fig. 1.
Phase diagram of the mean-field model. The vertical and horizontal axes representthe possible values of the rates φ aa and φ bb , respectively, ranging from zero to infinity. Inour picture, the rates φ ab and φ ba are constants whose fixed values are indicated at somearbitrary points. The spatially explicit model.
To understand the role of spatial structurecaused by local interactions, we now return to the spatially explicit modeland provide comparisons with the nonspatial model. We refer the readerto Figures 1 and 2 that show the phase diagrams of the mean-field modeland the interacting particle system, respectively. In Figure 1, the regionlabelled gene a indicates the set of parameters for which the density u aa converges to 1, the region coexistence the set of parameters for which thereis a locally stable interior fixed point, and the region founder control the setof parameters for which the interior fixed point is unstable. In Figure 2, theregion labelled gene a indicates the set of parameters for which the processstarting with infinitely many individuals of type aa converges to the “all aa ” configuration and the region coexistence the set of parameters for whichthere is a stationary distribution with a positive density of both genes.In the case when φ aa and φ bb are small, the particle system has a sta-tionary distribution with a positive density of genes of type a and b , whichis symptomatic of the existence of a locally stable interior fixed point forthe mean-field model when φ aa < φ ba and φ bb < φ ab . This indicates thatcoexistence occurs for both models. N. LANCHIER AND C. NEUHAUSER
Fig. 2.
Spatially explicit model.
On the other hand, in the case when φ ab and φ ba are small, coexistenceis not possible. This, however, translates differently for the spatial and non-spatial models. In the spatial model, we believe that, for almost all theparameters, there is fixation of one of the genes. Provided that the processstarts with infinitely many genes of both types, the “winner” depends on theparameters of the model (see Figure 2). Simulations suggest that, along thetransition curve, clustering occurs in any dimension (see the right-hand sideof Figure 3 for a picture in d = 2), and that the three regimes meet at point φ aa = φ ba and φ bb = φ ab as indicated in Figure 2. In contrast, the mean-field model exhibits founder control when φ aa > φ ba and φ bb > φ ab , that is,the interior fixed point is unstable and the “winner” depends on both theparameters and the initial densities (see the right-hand side of Figure 4).Numerical simulations also indicate that the set of parameters for whichcoexistence may occur is sensibly smaller for the spatial model. We willprove rigorously the existence of a set of parameters for which the interiorfixed point is locally stable for the mean-field model whereas coexistence isnot possible for the 1-dimensional stochastic model. This reduction of thecoexistence region of the particle system when compared to the mean-fieldmodel has already been proved theoretically for a spatially explicit versionof the Lotka–Volterra model [see Neuhauser and Pacala (1999)]. We nowstate our results in detail starting with the coexistence result. PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL Theorem 1 (Coexistence).
Assume that φ ab = φ ba > . If φ aa and φ bb are sufficiently small then there exists a translation invariant stationary dis-tribution ν such that ν ( η ( x ) = ab ) > . The same holds when φ ab < φ ba and φ ba < φ ab in the 1-dimensional case. See the left-hand side of Figure 3 for a picture. In the case when φ ab = φ ba = 1 and φ aa = φ bb = 0, the process η t reduces to an annihilating branch-ing process (ABP) in which 0 = homozygote of either type and 1 = het-erozygote, that is, the state of each site of the lattice flips from 0 to 1 andfrom 1 to 0 at rate equal to the number of neighbors in state 1. The mainingredient to prove Theorem 1 is a rescaling argument that compares theABP viewed on suitable length and time scales with a 1-dependent orientedpercolation process on the lattice G = { ( z, n ) ∈ Z : z + n is even and n ≥ } . More precisely, by defining occupied sites in the lattice G as the set of sites( z, n ) ∈ G such that there is at least one particle (state 1) in the spatial cube N ze + [ − N/ , N/ d = ( N z, , . . . ,
0) + [ − N/ , N/ d at time nT , we will prove the following theorem. Fig. 3.
Spatially explicit model with nearest neighbor interactions at time 50 on the × square with periodic boundary conditions. White (respectively, black) sites refer tohomozygotes of type a (respectively, b ), and grey sites to heterozygotes. Left: φ aa = φ bb = 4 and φ ab = φ ba = 5 . Right: φ aa = φ bb = 5 and φ ab = φ ba = 1 . In both pictures, the processstarts with a Bernoulli product measure, half of the genes being of type a . The picture on theleft shows the process at equilibrium but not the picture on the right in which clusters growin time indicating that the only equilibria are the “all aa ” and the “all bb ” configurations. N. LANCHIER AND C. NEUHAUSER
Fig. 4.
Solution curves of the mean-field model. The parameters are given by φ aa = 1 , φ ba = 3 , φ bb = 2 and φ ab = 4 (left), and φ aa = 4 , φ ba = 2 , φ bb = 3 and φ ab = 1 (right). Theorem 2 (Annihilating branching process).
For any ε > , the pa-rameters N and T can be chosen in such a way that the set of occupied sitesdominates the set of wet sites in a 1-dependent oriented site percolationprocess on G in which sites are open with probability − ε . This implies coexistence of both genes (survival of the heterozygotes)when φ aa = φ bb = 0. To conclude the proof of Theorem 1, we will apply astandard perturbation argument to extend the result to a region where φ aa and φ bb are positive but small.The last three theorems provide conditions under which there is fixationof the gene of type a and so extinction of the gene of type b . Here, “fixationof the gene of type a ” means that the process converges to the “all aa ”configuration. The reader will note that, by symmetry of the evolution rulesof the process, the conditions under which there is fixation of the gene oftype b can be deduced by reversing the roles of genes of type a and b . Thefollowing results hold for the process starting with infinitely many genes oftype a or infinitely many individuals of genotype aa . We can weaken thiscondition and assume instead that the process starts with finitely manygenes of type a or finitely many individuals of genotype aa , and concludethat fixation of the gene of type a occurs with positive probability. Theorem 3 (Fixation).
Assume that min( φ aa , φ ab ) > max( φ ba , φ bb ) andthat infinitely many genes of type a are present at time 0. Then there isfixation of the gene of type a . PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL Theorem 4 (Fixation).
Assume that φ ab is fixed, φ aa > d r d ! φ bb (4) and that infinitely many individuals of genotype aa are present at time 0.Then there is fixation of the gene of type a provided φ ba is sufficiently small. By assuming that φ aa = φ ab and φ ba = φ bb , the birth rate of a gene (ofeither type) at site x does not depend upon the type of the second genepresent at site x . This makes our process a biased voter model and allowsus to define naturally a dual process. In all the other cases, duality is nottractable. The proofs of Theorems 3 and 4 are made difficult by this lack ofa dual process. To establish Theorem 3, we will rely on a coupling argumentto compare the process with a biased voter model in which genes of type a are more competitive than genes of type b . The proof of Theorem 4 is moredifficult. The idea is to prove that, under condition (4) and for the processstarting with all sites in F N = [ − N, N ] d in state aa , the probability thatgenes of type b reach site 0 decreases exponentially fast with the length side N . This will be estimated by studying the evolution of the process along allthe paths starting outside F N and ending at site 0. Then we will use theindividual of genotype aa at site 0 as a source to produce a large cube voidof genes of type b and will conclude by comparing the process viewed onsuitable length and time scales with oriented percolation. A time rescalingin Theorem 4 suggests that φ ba being fixed if the parameters φ aa and φ bb aresufficiently large, and if condition (4) holds, then there is fixation of the geneof type a . Unfortunately, this cannot be deduced from our result since φ ba in the statement of Theorem 4 has to be smaller than some constant that(a priori) depends upon the values of the parameters φ aa and φ bb . However,numerical simulations suggest the following. Conjecture 5 (Fixation). Assume that φ ab and φ ba are fixed and thatinfinitely many genes of type a are present at time 0. Also, assume thatthe ratio φ aa /φ bb > a provided φ aa is sufficiently large.Numerical simulations also indicate that when φ aa = φ bb > φ ab = φ ba (foundercontrol for the mean-field model), clustering occurs in any dimension for thespatial model. The right-hand side of Figure 3 gives an illustration in d = 2.Theorems 1 and 3 show similarities with the mean-field model. Theorem4 implies the existence of a parameter region for which fixation occurs forthe spatial model (the “winner” depends on the choice of the parameters)while the mean-field model exhibits founder control. In both models, co-existence cannot occur for these parameters. Similarly, Conjecture 5 and N. LANCHIER AND C. NEUHAUSER
Theorem 6, formula (5), below which confirms our conjecture in d = 1 ex-tend the parameter region for which one of the genes goes extinct eventually.On the contrary, Theorem 6, formula (6), shows a major difference betweenspatial and nonspatial models, namely the set of parameters for which co-existence occurs is smaller for the spatial model (we believe that this holdsin any dimension, though our proof relies heavily on geometric argumentstrue in d = 1 only). This might provide a theoretical explanation for thelow frequency of the t -haplotype in the case of the house mouse. To modelthe evolution of the house mouse, we first choose φ ba < φ ab where a rep-resents the t -haplotype to take into account the fact that the t -haplotypeis transmitted to a majority of the offspring. Since male offspring who arehomozygous for the t -haplotype are sterile, we also fix φ aa = 0, though ourmodel does not distinguish between male and female. Finally, we let thelast parameter be φ bb = (1 / φ ab + φ ba ), assuming that ab individuals areoverall equally fit as bb individuals. This set of parameters corresponds topoints of the phase diagram which are well inside the coexistence phase forthe nonspatial model (see Figure 1) but to points which are much closer tothe boundary of the coexistence phase where the frequency of type a shouldbe expected to be lower for the spatial model (see Figure 2). In particular,the spatial component might be responsible for the low frequency of the t -haplotype in the example of the house mouse. Theorem 6 (Fixation).
Assume that d = 1 and that infinitely manyindividuals of genotype aa are present at time 0. Then whenever one thefollowing two conditions holds: φ aa > φ bb + φ ba + p φ bb φ ba , (5) φ aa > max (cid:18) φ ba , φ ba − φ ab (cid:19) and φ bb > small , (6) there is fixation of the gene of type a . The rest of this paper is devoted to proofs. In Section 2, we investigatein detail the existence and the stability of the fixed points of the mean-fieldequations (1)–(3). In Sections 3–6, we will prove the results regarding thespatially explicit model.
2. The mean-field model.
This section is devoted to the study of themean-field model introduced in Section 1. First of all, by setting the right-hand sides of (1) and (3) equal to 0, we obtain2 u aa (2 φ bb u bb + φ ba u ab ) = u ab (2 φ aa u aa + φ ab u ab ) , u bb (2 φ aa u aa + φ ab u ab ) = u ab (2 φ bb u bb + φ ba u ab ) , PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL so that 4 u aa u bb = u ab . Note that the previous equation holds when 2 φ aa u aa + φ ab u ab = 0 as well since in this case u aa = u ab = 0. Substituting this into (1)and setting the right-hand side equal to 0, leads to du aa dt = ( φ aa − φ ba ) u aa u ab + (1 / φ ab u ab − φ bb u aa u bb = ( φ aa − φ ba ) u aa u ab − φ bb − φ ab ) u aa u bb = 0 , from which it follows that u aa = 0 or ( φ aa − φ ba ) u ab = 2( φ bb − φ ab ) u bb at equilibrium. Similarly, using (3), we find u bb = 0 or ( φ bb − φ ab ) u ab = 2( φ aa − φ ba ) u aa at equilibrium. The condition u aa = 0 gives rise to the trivial equilibrium u bb = 1, and the condition u bb = 0 to the trivial equilibrium u aa = 1. If u aa u bb = 0, then2( φ aa − φ ba ) u aa = ( φ aa − φ ba )( φ bb − φ ab ) u ab = 2( φ bb − φ ab ) u bb . Setting φ = ( φ aa − φ ba ) + ( φ bb − φ ab ), we get φ u aa = ( φ aa − φ ba ) u aa + 2( φ aa − φ ba )( φ bb − φ ab ) u aa + ( φ bb − φ ab ) u aa = ( φ bb − φ ab ) u bb + ( φ bb − φ ab ) u ab + ( φ bb − φ ab ) u aa = ( φ bb − φ ab ) . By using the symmetry of the model, we deduce that at equilibrium φ u aa = ( φ bb − φ ab ) ,φ u bb = ( φ aa − φ ba ) ,φ u ab = 2( φ aa − φ ba )( φ bb − φ ab ) . Therefore, the condition for the existence of a nontrivial equilibrium is givenby ψ := ( φ aa − φ ba ) × ( φ bb − φ ab ) > . The stability of the equilibria can be analyzed using standard linearizationtechniques. Namely, denoting by J eq the Jacobian matrix at eq , the equilib-rium is locally stable if and only if all the eigenvalues of J eq associated witheigenvectors v = ( v , v , v ) such that v + v + v = 0 are negative. We firstlook at the stability of u aa = 1. The Jacobian matrix is given by J a = φ aa − φ ba − φ bb − φ aa + φ ba φ aa + 2 φ bb − φ aa . N. LANCHIER AND C. NEUHAUSER
The eigenvalues of the matrix J a are 0, − φ aa + φ ba and − φ aa . The eigenspacesassociated with the eigenvalues 0 and − φ aa + φ ba are respectively, given bySpan((1 , , , − , u aa = 1 is locally stable if and only if φ aa > φ ba . Bysymmetry, the equilibrium u bb = 1 is locally stable if and only if φ bb > φ ab .Now, by letting ψ aa = ( φ aa − φ ba ) ,ψ bb = ( φ bb − φ ab ) ,ψ ab = ( φ aa − φ ba )( φ bb − φ ab ) , we find that the Jacobian matrix at the interior fixed point is given by J ab = 1 φ − ψ aa φ ab ψ ab ( φ bb + φ ab ) − ψ bb φ bb ψ aa ( φ aa + φ ab ) − ψ ab ( φ aa + φ bb + φ ab + φ ba ) 2 ψ bb ( φ bb + φ ba ) − ψ aa φ aa ψ ab ( φ aa + φ ba ) − ψ bb φ ba ! , where φ = ( φ aa − φ ba ) + ( φ bb − φ ab ) as above. First of all, it is easy to checkthat 0 is an eigenvalue, and that it is associated to the eigenvector(( φ bb − φ ab ) , φ aa − φ ba )( φ bb − φ ab ) , ( φ aa − φ ba ) )whose coordinates add up to 0 only if ψ ≤
0. In particular, since the condition ψ > { ( u aa , u ab , u bb ) : u aa = p , u ab = 2 p (1 − p ) ,u bb = (1 − p ) for some p ∈ [0 , } . This curve is clearly visible in Figure 4. A population of diploid individualsin which the genotype frequencies lie on Γ is said to be in Hardy–Weinbergequilibrium. When Mendel’s law of segregation holds, points in this curvecorrespond to the limiting genotype frequencies of a nonspatial population(mean-field model), with p denoting the initial frequency of genes of type a .Note however that the inclusion of local interactions translates into spatialcorrelations that lead to an increase of the density of homozygous indi-viduals, and a decrease of the density of heterozygous individuals. For thegeneral mean-field model (when Mendel’s law does not hold), points in Γare no longer equilibrium points, but the curve is globally invariant underthe time evolution. Observing in addition that Γ contains the interior fixedpoint (take p = ( φ bb − φ ab ) φ − ), we find that(( φ bb − φ ab ) , ( φ aa − φ ba ) − ( φ bb − φ ab ) , − ( φ aa − φ ba ))is an eigenvector, and that it is associated to the eigenvalue φ := ( φ aa − φ ba )( φ bb − φ ab )(( φ aa − φ ba ) + ( φ bb − φ ab )) . This shows that, in the direction of the curve Γ, the interior fixed point is:
PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL
1. locally stable if φ aa < φ ba and φ bb < φ ab ,2. unstable if φ aa > φ ba and φ bb > φ ab .This proves the instability when φ aa > φ ba and φ bb > φ ab . See Figure 4 for apicture of the solution curves. Finally, using that 0 and φ are eigenvalues,it is straightforward to conclude that (as the third root of the characteristicpolynomial) the third eigenvalue is given by φ := − φ aa − φ ba )( φ bb − φ ab )( φ aa + φ bb )+ ( φ bb − φ ab ) φ ba + ( φ aa − φ ba ) φ ab )which is always negative whenever ψ = ( φ aa − φ ba )( φ bb − φ ab ) >
0, the con-dition required for the existence of an interior fixed point. This proves thestability when φ aa < φ ba and φ bb < φ ab .
3. Proof of Theorem 1.
The aim of this section is to prove that bothalleles may coexist provided φ aa and φ bb are sufficiently small. We will startby assuming that φ ab = φ ba > φ ab < φ ba and φ ba < φ ab .First, we observe that when φ ab = φ ba > φ aa = φ bb = 0, the local dy-namics of heterozygous sites can be defined regardless of the type of nearbyhomozygous sites. More precisely, by setting 0 = homozygote of either typeand 1 = heterozygote, or, equivalently, by letting ¯ η t denote the stochasticprocess defined for any x ∈ Z d and any time t ≥ η t ( x ) = (cid:26) , if η t ( x ) = ab ,0 , otherwise , (7)we obtain a new Markov process with state space { , } Z d whose state ateach lattice site flips from 0 to 1 and from 1 to 0 at rate φ ab times thenumber of neighbors in state 1. This process is known as the annihilatingbranching process (ABP) and has been investigated by Sudbury (1990) andBramson, Ding and Durrett (1991). Their results show that starting withany measure, which puts no mass on the empty configuration, the ABPconverges in distribution to the Bernoulli product measure with density 1 / a and b coexist when φ aa = φ bb = 0. To proveTheorem 1, we will employ a rescaling argument [see Bramson and Durrett(1988)] where we tile the space–time region Z d × Z + into large boxes, definea “good event” in each box and show that if a good event occurs in onesuch box, it will with probability close to 1 occur in neighboring boxes ata later time. Comparison with an oriented percolation process on Z × Z + then completes the proof. The “good event” is the presence of heterozygousindividuals in a given box if one thinks of the process η t , or the presence of N. LANCHIER AND C. NEUHAUSER η t . To show that heterozygositycan spread, we define a repositioning algorithm that shows that heterozygousindividuals spread out linearly in the direction of any of the coordinate axes(Lemma 3.1). The reason for invoking the rescaling argument is that oncethe process is embedded in a block structure, it is straightforward to apply aperturbation argument and deduce that coexistence occurs as well when φ aa and φ bb are sufficiently small. The main ingredient to compare the ABP witha percolation process is that particles can be moved along a given straightline ∆ while staying, with probability close to 1, within a reasonably shortdistance of ∆ (see Lemma 3.1 below).For i ∈ { , , . . . , d } , we let e i denote the i th unit vector, π i denote theorthogonal projection on the i th axis, and H i denote the hyperplane orthog-onal to e i , that is, π i ( x ) = h x, e i i = x i and H i = Span( e , e , . . . , e i − , e i +1 , . . . , e d ) , where h · , · i is the usual inner product on R d . We also set H i,k = [ − k . , k . ] d ∩ H i and E i,k = k [ j = − k ( H i,k + je i ) , where k is a large integer whose value will vary (increase) from lemma tolemma. For ω ∈ {− , +1 } and j = 1 , , . . . , k , we introduce the hitting times τ ω,i,j,k = inf { t ≥ η t ( x ) = 1 for some x ∈ H i,k + ωje i } and τ ω,i,k = sup j =1 , ,...,k τ ω,i,j,k . Roughly speaking, the next lemma tells usthat, with probability close to 1 when k is large, a particle at site 0 canbe moved of a distance k along the i th axis while staying within a distance k . of this axis. This can be done in less than Bk units of time for someconstant B that depends on the parameters of the process but not on k . Lemma 3.1.
Let ¯ η (0) = 1 and i ∈ { , , . . . , d } . There exist B < ∞ , C < ∞ and α > independent of k such that P ( τ ω,i,k > Bk ) ≤ C exp( − α k . ) and P (¯ η t ( x ) = 0 for all x ∈ E i,k and some t ≤ τ ω,i,k ) ≤ C exp( − α k . ) for k sufficiently large. Proof.
Note that, by symmetry, it suffices to prove Lemma 3.1 for ω = +1. The idea is to construct a jump process X t that keeps track ofa particle moving towards the hypercube H i,k + ke i while staying withina distance k . of the i th axis. The process X t starts at X = 0 and theevolution rules are formulated according to whether the process moves inthe direction of the vector e i or in the hyperplane X t + H i as follows: PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL
1. Whenever ¯ η t ( X t + e i ) = 1, the process X t jumps instantaneously to X t + e i . In particular, site X t + e i is empty, that is, in state 0, at any time. Ifthe particle at X t is killed by a particle located at X t − e i (due to the birthof a gene originated from site X t − e i ), then X t jumps instantaneously to X t − e i .2. If the particle at site X t gives birth to a particle which is sent to x ∈ X t + H i , then X t jumps to site x . If the particle at site X t is killed bya particle located at x ∈ X t + H i (due to the birth of a gene originatedfrom site x ∈ X t + H i ), then X t jumps to x .Rules 1 and 2 are illustrated in Figure 5 (note that the definition of X t implies that at time t site X t is always occupied). The picture shows aconfiguration of the process ¯ η t together with the position of X t . Open circlesrefer to empty sites (state 0) and full circles to occupied sites (state 1). Ifa particle gives birth through a j -arrow, then the process X t jumps to site x j , j = 1 , , ,
4. Births through 1- and 2-arrows illustrate rule 1 (jumps inthe direction of e i ), while births through 3- and 4-arrows illustrate rule 2(jumps in the direction of H i ). Lemma 3.1 is a straightforward consequenceof Lemmas 3.2 and 3.3 in which we investigate the projections of the process X t on the i th axis and the hyperplane H i , respectively.In the next lemma, we prove that rule 1 above implies that, with proba-bility close to 1 for k large, the projection π i ( X t ) reaches k in less than Bk units of time and before reaching − k . Fig. 5.
Evolution rules of the process X t . N. LANCHIER AND C. NEUHAUSER
Lemma 3.2.
Let σ ω,i = inf { t ≥ π i ( X t ) ≥ k } . Then there exists B < ∞ such that, for k sufficiently large, P ( σ ω,i > Bk ) ≤ C exp( − α k ) and P ( π i ( X t ) ≤ − k for some t ≤ σ ω,i ) ≤ C exp( − α k ) for suitable constants C < ∞ and α > . Proof.
We observe that π i ( X t ) increases by one whenever a particleis sent to site X t + e i (type 1 arrow in Figure 5), which occurs at rate atleast φ ab = φ ba since X t is occupied. Rule 1 above also implies that theprocess π i ( X t ) can decrease by one only if X t − e i is occupied, the decreaseresulting from the birth at X t of a gene that originated from site X t − e i (type 2 arrow in Figure 5). This occurs (when X t − e i is occupied) at rate φ ab / φ ba / φ ab . In other words, if we introduce the gap process G t = inf { j ≥ η t ( X t − ( j + 1) e i ) = 1 } , then π i ( X t ) jumps to the interval [ π i ( X t ) + 1 , ∞ ) at rate at least φ ab re-gardless of the value of the process G t and jumps to π i ( X t ) − φ ab only when G t = 0. In addition, a straightforward calculation shows that G t makes transitions as follows:0 → I at rate at least 2 φ ab ,I → dφ ab , where I = { , , . . . } and where “0 → I ” (respectively, “ I → m (respectively, from m to 0) for some m ∈ I .When site X t − e i is empty, the first rate is the rate at which particles at X t − e i and X t kill each other. When X t − e i is occupied, it is the rate atwhich particles at X t − e i and X t kill the particle at site X t − e i . The upperbound for the second rate is reached when all the neighbors of sites X t − e i and X t + e i are occupied. The evolution rules imply that P ( G t + h = 0) ≤ P ( G t = 0) × [1 − φ ab h ] + P ( G t = 0) × dφ ab h + O ( h ) . By subtracting P ( G t = 0) from both sides, dividing by h and taking thelimit as h → ddt P ( G t = 0) ≤ − φ ab × P ( G t = 0) + 4 dφ ab × P ( G t = 0) ≤ dφ ab − (4 d + 2) φ ab P ( G t = 0)from which it follows that P ( G t = 0) ≤ d d + 1 + (cid:20) P ( G = 0) − d d + 1 (cid:21) × exp( − (4 d + 2) φ ab t ) . PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL Taking the lim sup as t → ∞ , we getlim sup t →∞ P ( G t = 0) ≤ d d + 1 . Finally, by using that π i ( X t ) increases by one at rate at least φ ab , anddecreases by one at rate at most φ ab when G t = 0 and at rate 0 when G t = 0,we obtainlim inf t →∞ π i ( X t ) t ≥ φ ab − φ ab × lim sup t →∞ P ( G t = 0) ≥ φ ab d + 1 . (8)In particular, by setting B = 4 dφ − ab , we can conclude that P ( σ ω,i > Bk ) ≤ P ( π i ( X t ) ≤ k for all t ≤ Bk ) ≤ P ( π i ( X Bk ) ≤ k )= P ( π i ( X Bk ) ≤ ( φ ab / d ) × Bk ) ≤ C exp( − α k )for suitable C < ∞ and α > k sufficiently large. The limit in (8) alsoimplies that P ( π i ( X t ) ≤ − k for some t ≤ σ ω,i ) ≤ P ( π i ( X t ) ≤ − k for some t ≥ ≤ C exp( − α k )for suitable C < ∞ and α > k sufficiently large. The lemma follows. (cid:3) The next lemma shows that rule 2 describing the evolution of X t impliesthat the norm of the projection of X t on the hyperplane H i is bounded by k . until time Bk . Lemma 3.3.
Let B be given by Lemma 3.2 and j ∈ { , , . . . , i − , i +1 , . . . , d } . Then, for k sufficiently large, P ( π j ( X t ) / ∈ [ − k . , k . ] for some t ≤ Bk ) ≤ C exp( − α k . ) for suitable C < ∞ and α > . Proof.
Let j = i . We observe that, when site x = X t + e j is empty, theprocess X t jumps to x whenever the particle at X t gives birth to a particlewhich is sent to x (type 3 arrow in Figure 5), which occurs at rate φ ab . Whensite x is occupied, the process X t jumps to x whenever the particle at x killsthe particle at X t (type 4 arrow in Figure 5), which also occurs at rate φ ab .The same holds by replacing X t + e j by X t − e j . This implies that π j ( X t )is a continuous-time symmetric random walk starting at π j ( X ) = 0 run atrate φ ab . This is the main ingredient to prove the lemma. Let M k denotethe number of times π j ( X t ) jumps by time Bk , and Y n the discrete-timeversion of π j ( X t ) so that π j ( X Bk ) = Y M k . The expectation of M k is given N. LANCHIER AND C. NEUHAUSER by m k = 2 Bkφ ab . By decomposing the event to be estimated according towhether M k > m k or M k ≤ m k , we obtain P ( π j ( X t ) / ∈ [ − k . , k . ] for some t ≤ Bk ) ≤ P ( M k > m k ) + P ( Y n / ∈ [ − k . , k . ] for some n ≤ m k ; M k ≤ m k ) . Large deviation estimates imply that the first term can be bounded as fol-lows: P ( M k > m k ) ≤ C exp( − α k )for suitable C < ∞ and α >
0. The second term can be estimated by firstusing the reflection principle and then Chebyshev’s inequality. Given θ > P ( Y n / ∈ [ − k . , k . ] for some n ≤ m k ; M k ≤ m k ) ≤ P ( Y m k / ∈ [ − k . , k . ]) ≤ − θk . ) E exp( θY m k ) ≤ − θk . ) m k Y n =1 E exp( θ ( Y n − Y n − )) ≤ − θk . + 2 m k log φ ( θ )) , where φ ( θ ) is the moment generating function of Y . Since E Y = 0 andVar Y < ∞ , we have that log φ ( θ ) ≤ C θ for some C < ∞ and for θ smallenough. By taking θ = k − . in the last expression, we can conclude that P ( Y n / ∈ [ − k . , k . ] for some n ≤ m k ; M k ≤ m k ) ≤ − k . + 4 C φ ab Bk − . ) ≤ − k . / k sufficiently large. This completes the proof. (cid:3) To deduce Lemma 3.1, we first observe that Lemmas 3.2 and 3.3 implythat P ( τ ω,i,k > Bk ) ≤ P ( σ ω,i > Bk )+ P ( π j ( X t ) / ∈ [ − k . , k . ] for some t ≤ Bk and some j = i ) ≤ C exp( − α k ) + ( d − × C exp( − α k . ) . By Lemmas 3.2 and 3.3, we also have P ( η t ( x ) = 0 for all x ∈ E i,k and some t ≤ τ ω,i,k ) ≤ P ( π i ( X t ) ≤ − k for some t ≤ σ ω,i )+ P ( π j ( X t ) / ∈ [ − k . , k . ] for some t ≤ σ ω,i and some j = i ) ≤ P ( π i ( X t ) ≤ − k for some t ≤ σ ω,i ) + P ( σ ω,i > Bk )+ P ( π j ( X t ) / ∈ [ − k . , k . ] for some t ≤ Bk and some j = i ) ≤ C exp( − α k ) + ( d − C exp( − α k . ) PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL which concludes the proof of Lemma 3.1. (cid:3) Note that rule 1 describing the evolution of X t allows to give the process X t a drift in the direction of the vector e i whereas rule 2 implies that theprojection of X t on the orthogonal hyperplane H i evolves according to arandom walk without drift. However, rule 2 can be defined in such a waythat the tagged particle has an additional drift towards the i th axis. Thisapproach would probably strengthen Lemma 3.1 as follows: a particle at site0 can be moved of a distance k along the i th axis while staying within adistance independent of k of the axis. This result is obviously more difficultto prove than our Lemma 3.1 but would probably shorten the proofs ofLemmas 3.4 and 3.5 below.The next step is to show that, for given ε >
0, the particle system, whenviewed on suitable length and time scales, dominates (in a sense to be madeprecise) the set of wet sites of an oriented percolation process on the lattice G = { ( z, n ) ∈ Z : z + n is even and n ≥ } in which sites are open with probability 1 − ε [see Durrett (1984), for moredetails about oriented percolation]. For κ = 1 , , . . . , let J κ = [ − N/κ, N/κ ) d where N is a large integer to be fixed later and let T >
0. A site ( z, n ) ∈ G is said to be occupied if there is at least one particle in the translated cube N ze + J at time nT , that is,¯ η nT ( x ) = 1 for some site x ∈ N ze + J . The objective is to show that if site ( z, n ) ∈ G is occupied then ( z − , n + 1)and ( z + 1 , n + 1) are occupied with probability at least 1 − ε . Since theevolution rules of the particle system are translation invariant in space andtime, it suffices to prove the result when z = n = 0. To prepare for therescaling argument, we need the following lemma, which will be appliedrepeatedly and in which we prove that, with probability close to 1 when N is large, a particle in J κ can be moved to the smaller cube J κ in less than dBN units of time and before stepping out of J κ . Lemma 3.4.
Let κ be a positive integer and consider the stopping times υ in = inf { t ≥ η t ( x ) = 1 for some x ∈ J κ } and υ out = inf { t ≥ η t ( x ) = 0 for all x ∈ J κ } . If there is at least one particle in J κ at time 0, then P ( υ in ≤ dBN ; υ in < υ out ) ≥ − C exp( − α N . ) , where B is given by Lemma 3.1 and for suitable C < ∞ and α > . N. LANCHIER AND C. NEUHAUSER
Proof.
The idea is to apply Lemma 3.1 d times in each direction i ∈{ , , . . . , d } to move a particle from anywhere in J κ to J κ . We refer thereader to Figure 6 for an illustration. To make the argument precise, we let X ∈ J κ with ¯ η ( X ) = 1 and, for j = 0 , , . . . , d , set D j = [ − jN . , jN . ] d + π j +1 ( X ) e j +1 + π j +2 ( X ) e j +2 + · · · + π d ( X ) e d ,υ j = inf { t ≥ η t ( x ) = 1 for some x ∈ D j } and ω j = − sign( π j ( X )). Note that the sequence starts at D = X . By ob-serving that D d ⊂ J κ for all N sufficiently large and ∅ = { x + E i,N/ κ } ∩ H i ⊂ D i for all x ∈ D i − and applying Lemma 3.1 with i = 1 , , . . . , d , k = N/ κ and ω = ω i , we obtain P ( υ in > dBN ) ≤ P ( υ d > dBN/ κ ) ≤ d X i =1 P ( υ i − υ i − > BN/ κ ) ≤ d X i =1 P ( τ ω i ,i,k > BN/ κ ) ≤ C exp( − α N . ) , where C = C × d and α = α / (2 κ ) . . Since x + E i,N/ κ ⊂ J κ for all x ∈ D i − and the evolution rules of the process are translation invariant in spaceand time, we can also deduce that P ( υ in > υ out ) ≤ P ( υ d > υ out ) ≤ d X i =1 P ( υ i > υ out | υ i − < υ out ) ≤ d X i =1 sup x ∈ D i − P (¯ η t ( z ) = 0 for all z ∈ x + E i,N/ κ Fig. 6.
Picture of the block argument.
PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL and some t < υ i | ¯ η ( x ) = 1) ≤ d X i =1 P (¯ η t ( z ) = 0 for all z ∈ E i,N/ κ and some t < τ ω i ,i,k | ¯ η (0) = 1) ≤ C exp( − α N . ) , where C < ∞ and α > (cid:3) The next lemma provides the main estimate for the rescaling argument.
Lemma 3.5.
Let T = (3 d + 1) BN and assume that ¯ η ( X ) = 1 for some X ∈ J . Then P (¯ η T ( x ) = 0 for all x ∈ J + N e ) ≤ C exp( − α N . ) for suitable C < ∞ and α > . Proof.
For any integer i ≥ , we introduce the stopping times T j = (cid:26) inf { t > T j − : ¯ η t ( x ) = 1 for some x ∈ J + N e } , for j odd,inf { t > T j − : ¯ η t ( x ) = 0 for all x ∈ J + N e } , for j even , with the convention T = 0. In words, the stopping times T j ’s tell us alter-nately when a particle appears in J + N e and when J + N e becomesempty. Let T ⋆ be the first exit time T ⋆ = inf { t ≥ T : ¯ η t ( x ) = 0 for all x ∈ J + N e } . Then, to establish the lemma, it suffices to prove that T < T and T ⋆ > T with probability arbitrarily close to 1 when N is large since this indicatesthat a particle appears in J + N e by time T and then stays inside J + N e until time T .To estimate the event that T < T , we start by applying successivelyLemma 3.4 to bring a particle from J to J , then from J to J , andfinally from J to J in less than 3 dBN units of time. More precisely, wehave P (¯ η t ( x ) = 0 for all x ∈ J and all t ≤ dBN ) ≤ × C exp( − α N . ) . Once there is a particle in J , we apply Lemma 3.1 with i = 1 and k = N to conclude that P ( T > T ) = P ( T > (3 d + 1) BN ) ≤ × C exp( − α N . ) + C exp( − α N . )(9) ≤ C exp( − α N . ) N. LANCHIER AND C. NEUHAUSER for N large enough so that J + [ − N . , N . ] d ⊂ J and suitable C < ∞ and α > T ⋆ > T , which is the event that box J + N e contains at least one particle from T to T . Due to nearest neigh-bor interactions, there is at least one particle within distance 1 from theboundary of J + N e , and so contained in J + N e , at time T j − when j is odd. In particular, by applying Lemma 3.4 with κ = 4, we get P ( T j > T ⋆ | T j − < T ⋆ ) ≤ C exp( − α N . ) for j odd . (10)Now, the definition of the T j ’s implies that, between time T j − and time T j for j even, there is at least one particle in the box J + N e . It follows that P ( T j > T ⋆ | T j − < T ⋆ ) = 0 for j even . (11)Roughly speaking, (10) and (11) tell us that the number of times a particleappears in J + N e and J + N e becomes empty before J + N e becomesempty can be made arbitrarily large by choosing N large. We still need toprove that the number of steps up to time T is smaller than a constantuniform in N by estimating K = inf { j ≥ T j > T } which is about twice the number of times J + N e becomes empty bytime T . Since particles die at rate at most 2 dφ ab and the number of deathsrequired to empty J + N e starting with at least one particle in J + N e is larger than N/
32, we obtain E [ T j − T j − ] ≥ N × (64 dφ ab ) − for j even . Using the fact that T is equal to a constant times N , large deviation es-timates imply the existence of some m = m ( d, φ ab ) independent of N suchthat P ( K > m ) ≤ C exp( − α N )(12)for suitable C < ∞ and α > P (¯ η T ( x ) = 0 for all x ∈ J + N e ) ≤ P ( T > T ) + P ( T ⋆ < T ) ≤ P ( T > T ) + K X j =2 P ( T j > T ⋆ | T j − < T ⋆ ) ≤ P ( T > T ) + P ( K > m ) + m X j =2 P ( T j > T ⋆ | T j − < T ⋆ ) ≤ C exp( − α N . ) + C exp( − α N ) + m × C exp( − α N . ) . PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL This completes the proof of the lemma. (cid:3)
Lemma 3.5 implies that, for any ε >
0, the parameters N and T canbe chosen in such a way that, provided that site (0 ,
0) is occupied, (1 , − ε . Now, since the range of theinteractions is finite, the events described in Lemmas 3.2–3.5 have a finiterange of dependence. In particular, with probability at least 1 − ε for N sufficiently large, the configuration of the process at time T in the box N e +[ − N, N ] d only depends upon the graphical representation in the space–timeregion B , = ( N e , T ) + { [ − N, N ] d × [0 , T ] } . This implies the existence of an event G , that has probability at least1 − ε , which is measurable with respect to the graphical representation ofthe process restricted to B , , and such that, on the event that site (0 ,
0) isoccupied and G , occurs, site (1 ,
1) is occupied. Finally, since the evolutionrules of the process are translation invariant in space and time, we concludethat there exists a collection of events { G z,n : ( z, n ) ∈ G} such that, for any ε >
0, the parameters N and T can be chosen in such a way that:1. The event G z,n is measurable with respect to the graphical representationof the process restricted to the space–time region B z,n = ( N ze , nT ) + { [ − N, N ] d × [0 , T ] } .2. P ( G z,n ) ≥ − ε .3. If ( z, n ) is occupied and G z,n occurs, then ( z − , n + 1) and ( z + 1 , n + 1)are occupied.This implies that the set of occupied sites dominates the set of wet sites ina 1-dependent oriented site percolation process on G in which sites are openwith probability p = 1 − ε [see Durrett (1995) for a rigorous proof]. Thiscompletes the proof of Theorem 2.To deduce Theorem 1, the last step is to apply a perturbation argumentto prove that the previous property holds when φ aa and φ bb are positive butsmall. The probability that a gene originated from a homozygote is sent tothe space–time box B z,n can be bounded by X x ∈ [ − N, N ] d Z T φ exp( − φs ) ds = (4 N + 1) d × (1 − exp( − φT )) , (13)where φ = 2( φ aa + φ bb ). The parameters N and T being fixed so that, when φ aa = φ bb = 0, the set of occupied sites dominates the set of wet sites in a1-dependent percolation process with parameter 1 − ε , the rates φ aa and φ bb can be chosen so small that the probability in (13) is smaller than ε .With this choice, the set of occupied sites now dominates the set of wet sites N. LANCHIER AND C. NEUHAUSER in a percolation process with parameter 1 − ε . To conclude, we take ε > ab ” configuration. Then the limit ν is a stationary distribution.Moreover, the process being coupled with a supercritical oriented percolationprocess, the density of heterozygotes under the measure ν is strictly positive.This completes the proof of Theorem 1 when φ ab = φ ba .To conclude this section, we now prove Theorem 1 when φ ab = φ ba and d = 1. Since the assumption φ ab = φ ba has only been used in the proof ofLemma 3.1, it suffices to extend Lemma 3.1 (see Lemma 3.7 below) to the1-dimensional case when φ ab < φ ba and φ ba < φ ab . The difficulty arisingfrom the fact that φ ab = φ ba is that births of heterozygous individuals ontosites in state aa occur at rate φ ab whereas births of heterozygous individualsonto sites in state bb occur at rate φ ba . Without loss of generality, we canassume from now on that φ ab ≥ φ ba . To facilitate the writing of the proof,we will use as above the process ¯ η t defined in (7). Note however that, dueto the general assumption φ ab = φ ba , the process ¯ η t is no longer a Markovprocess. It should be thought of as an auxiliary process that allows us tokeep track of heterozygous sites. Let X + t = sup { x ∈ Z : ¯ η t ( x ) = 1 } and X − t = inf { x ∈ Z : ¯ η t ( x ) = 1 } denote the right and left edges of the process at time t . Then the followinglemma holds. Lemma 3.6.
Assume that φ ab < φ ba and start the process with a singleparticle at site 0. Then there exists a constant C > depending only onthe rates φ ab and φ ba such that lim inf t →∞ t − X + t > C and lim inf t →∞ t − X − t < − C . Proof.
By symmetry of the evolution rules, it suffices to prove thelemma for X + t only. The proof is similar to the proof of Lemma 3.2. Weintroduce the gap process G t = inf { j ≥ η t ( X + t − j −
1) = 1 } , Table 1
Evolution rules of G t . In our table, I k = { k, k + 1 , . . . } , ◦ = state 0, and • = state 1.( → I indicates a jump from state 0 to some state in I , etc.) The rows refer to thecases when G t = 0 , G t = 1 and G t ≥ , respectively Configurations Transition Rate Drift ( X + t ) ◦••◦◦◦ · · · •••◦◦◦ · · · → I ≥ φ ab + φ ba ≥ ( φ ba − φ ab ) / ◦•◦•◦◦ · · · ••◦•◦◦ · · · → ≤ φ ab ≥ φ ba •◦◦•◦◦ · · · •◦◦◦•◦ · · · I → ≤ φ ab ≥ φ ba PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL and observe that X + t jumps to the left at rate ( φ ab + φ ba ) / G t = 0, andjumps to the right at rate at least φ ba regardless of the value of G t . Moreover,a straightforward calculation shows that the process G t makes transitionsas indicated in Table 1, where • ’s refer to 1’s and ◦ ’s to 0’s, and where I k denotes the set { k, k + 1 , . . . } . The transition from 0 to I is obtained bykilling one of the particles at X + t and X + t − X + t − X + t − X + t + 1. Since φ ab ≥ φ ba ,upper bounds are reached when each empty site is in state bb . These rulesimply that P ( G t + h = 0) ≤ P ( G t = 0) × [1 − ( φ ab + φ ba ) h ] + P ( G t = 0) × φ ab h + O ( h ) . The same argument as in Lemma 3.2 leads tolim inf t →∞ t − X + t ≥ φ ba − φ ab + φ ba × lim t →∞ P ( G t = 0) ≥ φ ba − ( φ ab + φ ba ) × φ ab φ ab + φ ba )= 2 φ ba + 5 φ ab φ ba − φ ab φ ab + φ ba ) . The roots of Q ( X ) = 2 X + 5 φ ab X − φ ab are − φ ab and φ ab /
2. This impliesthat the last term of the previous inequality is positive if and only if φ ab < φ ba and proves the lemma. (cid:3) We can now use Lemma 3.6 to demonstrate that with high probability 1’sspread out linearly in time and can be found within this growing cone. Thesecond part of Theorem 1 then follows as before from a rescaling argumentand we omit the details.
Lemma 3.7.
Assume that φ ab < φ ba and ¯ η (0) = 1 . For ω ∈ {− , +1 } ,we let τ ω = inf { t ≥ η t ( ωk ) = 1 } . Then there exist constants
B < ∞ , C < ∞ and α > such that P ( τ ω > Bk ) ≤ C exp( − α k ) and P (¯ η t ( x ) = 0 for all x ∈ [ − k, k ] and some t ≤ τ ω ) ≤ C exp( − α k ) . Proof.
By symmetry, it suffices to prove the result for ω = +1. Let Y + t = sup { x ≤ k : ¯ η t ( x ) = 1 } and Y − t = inf { x > k : ¯ η t ( x ) = 1 } . N. LANCHIER AND C. NEUHAUSER
Since Y + t and Y − t do not interfere by time τ ω , the process Y + t − Y +0 is equalin distribution to the right edge X + t of the particle system starting witha single particle at site 0 by time τ ω . In particular, using large deviationestimates, the fact that Y +0 ≥ P ( τ ω > Bk ) ≤ P ( X + t < k for all t ≤ Bk ) ≤ P ( X + Bk < k ) ≤ C exp( − α k )for B = 2 /C and suitable C < ∞ and α >
0. Lemma 3.6 also impliesthat P (¯ η t ( x ) = 0 for all x ∈ [ − k, k ] and some t ≤ τ ω ) ≤ P ( X + t < − k for some t ≥ ≤ C exp( − α k ) . This completes the proof of Lemma 3.7. (cid:3)
4. Proof of Theorem 3.
To prove Theorem 3, it is more convenientto consider our particle system as a set of genes evolving on Z d × { , } rather than a set of genotypes evolving on the lattice Z d . More precisely,we introduce the Markov process ξ t : Z d × { , } −→ { a, b } whose evolutionat site ( x, g ) is given by the following transition rates, where we have set h = g + 1 mod 2: a → b at rate 2 φ bb X k x − z k =1 { ξ ( z, g ) = ξ ( z, h ) = b } + φ ba X k x − z k =1 { ξ ( z, g ) = a, ξ ( z, h ) = b } + φ ba X k x − z k =1 { ξ ( z, g ) = b, ξ ( z, h ) = a } ,b → a at rate 2 φ aa X k x − z k =1 { ξ ( z, g ) = ξ ( z, j ) = a } + φ ab X k x − z k =1 { ξ ( z, g ) = a, ξ ( z, h ) = b } + φ ab X k x − z k =1 { ξ ( z, g ) = b, ξ ( z, h ) = a } . The process η t is a genotype-based model whereas the process ξ t is a gene-based model, however both describe the same genetic system in the sensethat they can be defined on the same probability space in such a way that,for any site x ∈ Z d and any time t > η t ( x ) = ξ t ( x, ξ t ( x,
1) or η t ( x ) = ξ t ( x, ξ t ( x, PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL provided the property holds at time 0. In view of (14), it suffices to proveour results for either one of the two processes. We will prove Theorem 3 forthe gene-based model.The proof of Theorem 3 relies on duality techniques supplemented with astandard coupling argument. Let φ a = min( φ aa , φ ab ) and φ b = max( φ ba , φ bb ).The first step is to construct the gene-based model and the biased votermodel with parameters φ a and φ b on the same probability space in such away that ξ t has more genes of type a than the biased voter model (denotedby ζ t later). To do this, we will construct both processes from the samecollections of independent Poisson processes, which is referred to as Harris’graphical representation [Harris (1972)]. We will conclude by showing that ζ t converges to the “all a ” configuration provided φ a > φ b .To construct the gene-based model ξ t we put, for all g, h ∈ { , } andall sites x , z ∈ Z d such that k x − z k = 1, a specific arrow from site ( x, g )to site ( z, h ) at the arrival times of independent Poisson processes withsuitable rates. Specifically, for i, j ∈ { a, b } , we put, at the arrival times ofa Poisson process with parameter φ ij , a type ij arrow from site ( x, g ) tosite ( z, h ) to indicate that if site ( x, g ) is occupied by a gene of type i andsite ( x, g + 1 mod 2) is occupied by a gene of type j then site ( z, h ) becomesoccupied by a gene of type i . This graphical representation allows us toconstruct the gene-based model starting from any initial configuration ξ .Figure 7 shows a realization of the graphical representation of the gene-basedprocess.To construct the biased voter model ζ t , we put, for all g, h ∈ { , } andall sites x, z ∈ Z d such that k x − z k = 1, a type i arrow from site ( x, g ) tosite ( z, h ) at the arrival times of independent Poisson processes with rate φ i , i ∈ { a, b } . This indicates that if site ( x, g ) is occupied by a gene of type i then[regardless of the gene present at site ( x, g + 1 mod 2)] site ( z, h ) becomesoccupied by a gene of type i .To construct both processes on the same probability space, we use Poissonprocesses with rates depending on the signs of φ aa − φ ab and φ bb − φ ba . Ta-ble 2 gives an explicit description of the graphical representation in the fourpossible cases. The type i arrows, i ∈ { a, b } , in the construction of the gene-based model have the same interpretation as the ones in the construction ofthe biased voter model. In particular, a type i arrow in the construction ofthe gene-based model can be seen as the superposition of a type ia arrowand a type ib arrow. In view of the well-known properties of the Poisson pro-cess, each of these graphical representations produces the desired flip ratesto construct the processes we are interested in. Moreover, it follows directlyfrom Table 2 that ζ t ( x, i ) = a = ⇒ ξ t ( x, i ) = a for all x ∈ Z d and i ∈ { , } (15) N. LANCHIER AND C. NEUHAUSER
Fig. 7.
Construction of the gene-based model from the Harris’ graphical representation.The thick lines refer to sites occupied by a gene of type a . at any time t > φ ij , i, j ∈ { a, b } , are no longerindependent.Assume that φ a > φ b and that infinitely many genes of type a are presentat time 0. The last step is to prove that the biased voter model ζ t : Z d ×{ , } −→ { a, b } with parameters φ a and φ b converges almost surely to the“all a ” configuration. The proof relies on duality techniques. We constructthe process by drawing a type a arrow (respectively, an unlabeled arrow)from site ( x, g ) to site ( z, h ) at the arrival times of a Poisson process withrate φ a − φ b (respectively, φ b ). Unlabeled arrows indicate that the gene atsite ( z, h ) is replaced by the gene at site ( x, g ) regardless of the type ofboth genes. The dual process is constructed by going backwards in time andjumps at the tips of unlabeled arrows while branches at the tips of type a arrows as indicated in Figure 8. More precisely, if ˆ ζ s denotes the dual processat dual time s then:1. if there is an unlabeled arrow from ( x, g ) to ( z, h ) at time t − s for some( z, h ) ∈ ˆ ζ s − , then ˆ ζ s can be deduced from ˆ ζ s − by removing site ( z, h ) andadding site ( x, g ) while PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL Table 2
Couplings of the gene-based model and the biased voter model
Rate Gene-based Biased voter φ aa ≤ φ ab and φ bb ≤ φ ba φ aa a −→ a −→ φ ab − φ aa ab −→ φ bb b −→ b −→ φ ba − φ bb ba −→ b −→ φ aa ≤ φ ab and φ bb ≥ φ ba φ aa a −→ a −→ φ ab − φ aa ab −→ φ ba b −→ b −→ φ bb − φ ba bb −→ b −→ φ aa ≥ φ ab and φ bb ≤ φ ba φ ab a −→ a −→ φ aa − φ ab aa −→ φ bb b −→ b −→ φ ba − φ bb ba −→ b −→ φ aa ≥ φ ab and φ bb ≥ φ ba φ ab a −→ a −→ φ aa − φ ab aa −→ φ ba b −→ b −→ φ bb − φ ba bb −→ b −→
2. if there is a type a arrow from ( x, g ) to ( z, h ) at time t − s for some( z, h ) ∈ ˆ ζ s − , then ˆ ζ s can be deduced from ˆ ζ s − by adding site ( x, g ).The convergence to the “all a ” configuration follows from the analysisof the dual process. A stronger version of this result is proved in Bramsonand Griffeath (1980, 1981) for the biased voter model on Z d . Their proof,however, easily extends to the process on Z d × { , } . This, together with(14) and (15), concludes the proof of Theorem 3.
5. Proof of Theorem 4.
Similar to the proof of Theorem 1, the proof ofTheorem 4 relies on a rescaling argument. First, we set F N = [ − N, N ] d for N ≥ P F N the law of the process starting with all sites in F N occupied by individuals of genotype aa , and all sites outside F N occupiedby individuals of genotype bb . The first step is to prove the following lemma. N. LANCHIER AND C. NEUHAUSER
Fig. 8.
Picture of the dual coalescing branching random walks.
Lemma 5.1.
Assume (4) and φ ba = 0 . Then P F N ( η t (0) = aa at some time t ≥ ≤ C exp( − α N ) for suitable constants C < ∞ and α > . The proof of Lemma 5.1 is made difficult by the lack of a dual process.To estimate the probability that a gene of type b reaches site 0, we willconstruct a path γ ⊂ Z d starting outside the cube F N and ending in theneighborhood of site 0 along which genes of type b spread. The constructionwill be done by going backwards in time. By investigating the evolution ofthe process along the invasion path γ by going forward in time, we will provethat the probability that such a path exists is smaller than a constant times c − K for some c > d , where K denotes the length of γ . We will conclude theproof of Lemma 5.1 by summing over all the possible paths starting at site x for some x ∈ Z d \ F N and ending in the neighborhood of site 0.Lemma 5.1 implies that, with probability arbitrarily close to 1 when N is large, the individual at site 0 is of type aa at any time t ≥
0. The secondstep is to use this individual as a source to fill the cube F N with genes oftype a . This occurs by some time T < ∞ with probability arbitrarily closeto 1 for N large. By applying the same arguments as in the proof of Lemma PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL F N instead of the individual at site 0, we willdeduce the following lemma. Lemma 5.2.
Assume (4) and φ ba = 0 . For any ε > , P F N ( η T ( x ) = aa for some x ∈ F N ) ≤ ε for a suitable integer N ≥ and a (deterministic) time T < ∞ . Lemma 5.2 implies the existence of a partition of Z d into cubes of sidelength 2 N + 1 such that the set of cubes filled with a ’s at time 2 nT , n ≥ n of a 1-dependent oriented perco-lation process with parameter 1 − ε . This can be proved by applying thesame arguments as in the proof of Theorem 1 (for more details, see ourconstruction after the proof of Lemma 3.5). To prove that genes of type a actually outcompete genes of type b , namely that there is an in-all-directionsexpanding region which is void of genes of type b , we apply a result fromDurrett (1992) that shows that sites that are not occupied do not perco-late when ε is close enough to 0. Since genes of either type cannot appearspontaneously, once a region is void of one type, this type can only reappearin the region through invasion from the outside. To complete the proof ofTheorem 4, the last step is to apply a perturbation argument similar to theone described in the proof of Theorem 1 to extend the result to the casewhen φ ba is small. We now prove Lemmas 5.1 and 5.2 in detail.5.1. Construction of the invasion path γ . The first step in proving Lemma5.1 is to demonstrate the existence of a path γ ⊂ Z d along which genes oftype b spread in the case when the stopping time T N = inf { t ≥ η t (0) = aa } is finite. We first let the process starting with genes of type aa in F N evolveuntil time T N and then do the construction inductively by going backwardsin time.We set S = T N . The definition of S implies that there is a gene of type b that is sent to site 0 at time S . This gene must have originated from asite Γ with k Γ k = 1. Since φ ba = 0, we deduce that site Γ is necessarily instate bb at time S . Then we let S = sup { s ≤ S : η s (Γ ) = aa } . Again, there is a site Γ with k Γ − Γ k = 1 which is in state bb at time S .Let S = sup { s ≤ S : η s (Γ ) = aa } . N. LANCHIER AND C. NEUHAUSER
This procedure allows us to construct a sequence (Γ i , S i ) such that time S i is almost surely finite and site Γ i is in state bb at time S i . The constructionstops at i = K where K = min { i ≥ η s (Γ i ) = aa for all s ∈ [0 , S i ] } . Note that this implies that Γ K / ∈ F N and that K ≥ N . Our constructionimplies that 0 < S K < S K − < · · · < S = T N . Since it will be more convenient to work forward in time, we construct theinvasion path γ by reversing the direction of time, that is we set γ ( i ) = Γ K − i and s i = S K − i for i = 0 , , . . . , K , and γ = ( γ (0) , γ (1) , . . . , γ ( K )).5.2. Evolution of the genotypes along γ . The next step is to define aprocess ω t starting at ω = 0 such that γ ( ω t ) follows the invasion of genesof type b in the cube F N . The process ω t is stochastically smaller thana process whose probability to reach K decreases exponentially in K (seeLemma 5.3 below). This allows us to conclude that the probability that thepath γ exists decreases exponentially fast with its length, which is the keyto prove Lemma 5.1.The process ω t starts at ω = 0 and has values in { , , . . . , K } . The tran-sition rates depend on whether site γ ( ω t + 1) is in state aa or not.1. If γ ( ω t + 1) is in state aa (at time t ), then ω t → ω t + 1 , when the state of γ ( ω t + 1) jumps to ab, sup { i ≤ ω t − η t ( γ ( i )) = aa } , when the state of γ ( ω t ) jumps to aa .2. If γ ( ω t + 1) is not in state aa (at time t ), then ω t → sup { i ≤ ω t − η t ( γ ( i )) = aa } when the state of γ ( ω t ) jumps to aa .Note that ω t is well defined until time T N and sup { i ≤ ω t − η t ( γ ( i )) = aa } ≥ γ exists. Our construction also implies that ω t jumps from i to i + 1 attime s i (see the left-hand side of Figure 9) so that ω t reaches state K at time s K = T N , and that site γ ( ω t ) is in state ab or bb up to time s K . To keep trackof the genotype, we define a process ¯ ω t with state space ( Z / ∩ ( − , K ] byletting ¯ ω t = ω t − { η t ( γ ( ω t )) = ab, η t ( γ ( ω t + 1)) = aa } . See the right-hand side of Figure 9. We now describe the process ¯ ω t in details.The transition rates of the process depend on the state of sites γ ( ω t ) and γ ( ω t + 1), and we distinguish four different cases that we call configurations1–4 (see Figure 10 for an illustration): PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL Fig. 9.
Illustration of ω t (on the left-hand side) and ¯ ω t (on the right-hand side) indashed lines from a realization of the particle system. In both pictures, the horizontal axesrepresents both the state space of γ (on the bottom) and the state space of the processes ω t and ¯ ω t (on the top). Vertical thick black lines refer to sites in state bb , thick grey lines tosites in state ab and dotted lines to sites in state aa .
1. Config. 1—site γ ( ω t ) is in state bb and site γ ( ω t + 1) is in state aa .In this case, ¯ ω t can either jump to ¯ ω t + 1 / ω t + 1 when the state ofsite γ ( ω t + 1) jumps to ab , or jump to ¯ ω t − / Fig. 10.
Pictures of the configurations 1–4. N. LANCHIER AND C. NEUHAUSER γ ( ω t ) jumps to ab which implies that¯ ω t → (cid:26) ≤ ¯ ω t + 1 , at rate at most 4 dφ bb ,¯ ω t − / , at rate at least 2 φ aa .
2. Config. 2—site γ ( ω t ) is in state bb and site γ ( ω t + 1) is not in state aa .In this case, the process ¯ ω t is frozen.3. Config. 3—site γ ( ω t ) is in state ab and site γ ( ω t + 1) is in state aa .In this case, ¯ ω t can either jump to ¯ ω t + 1 / γ ( ω t )jumps to bb or the state of site γ ( ω t + 1) jumps to ab , or jump to the leftwhen the state of site γ ( ω t ) jumps to aa which implies that¯ ω t → (cid:26) ¯ ω t + 1 / , at rate at most 6 dφ bb , ≤ ¯ ω t − / , at rate at least φ aa .
4. Config. 4—site γ ( ω t ) is in state ab and site γ ( ω t + 1) is not in state aa .In this case, the process ¯ ω t can only jump to the left.The process ¯ ω t is well defined and ≥ − / T N condition onthe event that the invasion path γ exists. However, when ω t = 0 and theconfiguration is of type 3, the state of γ (0) may still jump to aa at the rateindicated above, thus contradicting the existence of γ . This motivates, inorder to bound the probability that γ exists, the introduction of the Markovprocess X t with state space Z / X t → (cid:26) ⌊ X t ⌋ + 1 , at rate 6 dφ bb ,X t − / , at rate φ aa ,where ⌊·⌋ denotes the integer part. The choice of the transition rates of X t and the observation above imply that the probability that the invasion path γ exists is bounded by the probability that, starting at 0, the process X t reaches state K before state −
1. To estimate the latter, we consider thediscrete-time version of X t , which is the Markov process Y n with transitionprobabilities Y n +1 = (cid:26) ⌊ Y n ⌋ + 1 , with probability r , Y n − / , with probability l, (16)where the probabilities r (right) and l (left) are given by r = 6 dφ bb φ aa + 6 dφ bb and l = φ aa φ aa + 6 dφ bb . Then, denoting by τ i the first time the process Y n hits state i , we have P ( γ = ( γ (0) , γ (1) , . . . , γ ( K )) is an invasion path)(17) ≤ P ( τ K < τ − | Y = 0) , which decreases exponentially in K according to the following lemma. PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL Lemma 5.3.
Assume (4). Then there exists a constant c > d such that P ( τ K < τ − | Y = 0) = ( c − · ( c K +1 − − ≤ c − K . Proof.
Let P i denote the law of the process starting at Y = i . Thenthe transition probabilities of the process indicated in (16) imply that P i ( τ i − < τ i +1 ) = l ∞ X j =0 ( l · r ) j = l − lr and P i ( τ i − > τ i +1 ) = r − lr for any i ∈ Z . Let p i = P i ( τ K < τ − ) and w ( i ) = ( l /r ) i . By decomposingaccording to whether the process first jumps to i − i + 1, and using thefact that ( l + r ) w ( i ) = l w ( i −
1) + rw ( i + 1) , we deduce that w ( Y n ) is a martingale. The martingale stopping theoremthen implies that w ( i ) = p i w ( K ) + (1 − p i ) w ( − p i = ( w ( i ) − w ( − · ( w ( K ) − w ( − − = ( c i +1 − · ( c K +1 − − with c = l /r . This implies that P ( τ K < τ − ) = p = ( c − · ( c K +1 − − .Finally, since l /r > d if and only if φ aa > d φ bb ( φ aa + 6 dφ bb ) , if and only if φ aa > d r d ! φ bb which is condition (4) in Theorem 4; the lemma follows. (cid:3) Proof of Lemma 5.1.
First of all, combining Lemma 5.3 with (17),we obtain P ( γ = ( γ (0) , γ (1) , . . . , γ ( K )) is an invasion path) ≤ c − K with c > d given in Lemma 5.3. Now, let x ∈ ∂F L +1 = F L +1 \ F L for some L ≥ x and ending in the neighbor-hood of 0 by Θ x . By the triangle inequality, for any γ = ( γ (0) , γ (1) , . . . , γ ( K )) ∈ Θ x , L ≤ k γ ( K ) − γ (0) k ≤ K − X i =0 k γ ( i + 1) − γ ( i ) k = K, N. LANCHIER AND C. NEUHAUSER that is, each path in Θ x has length at least L . Furthermore, since eachsite has 2 d nearest neighbors, the number of possible paths of length K isbounded by (2 d ) K . It follows that P (there exists γ ∈ Θ x such that γ is an invasion path) ≤ ∞ X K = L P (there exists γ ∈ Θ x such that γ is an invasion path of length K ) ≤ ∞ X K = L (2 d/c ) K = cc − d × (2 d/c ) L . Finally, since the birth of a gene of type b at site 0 implies the existence ofan invasion path starting at some site x / ∈ F N for the process starting withall sites in F N occupied by individuals of genotype aa , by summing over allsites x / ∈ F N , we obtain P F N ( T N < ∞ ) ≤ X x/ ∈ F N P (there exists γ ∈ Θ x such that γ is an invasion path) ≤ ∞ X L = N X x ∈ ∂F L +1 P (there exists γ ∈ Θ x such that γ is an invasion path) ≤ ∞ X L = N ((2 L + 3) d − (2 L + 1) d ) × cc − d × (2 d/c ) L ≤ C (2 d/c ) N ≤ C exp( − α N )for suitable C < ∞ and α >
0. This completes the proof of Lemma 5.1.5.4.
Proof of Lemma 5.2.
Lemma 5.2 is a straightforward consequenceof Lemma 5.1. Let ε > T ⋆ N the first time the cube F N isvoid of genes of type b , that is T ⋆ N = inf { t ≥ η t ( x ) = aa for all x ∈ F N } . By applying Lemma 5.1 to each of the sites belonging to ∂F N +1 at time T ⋆ N , we obtain P ( η t ( x ) = aa for some x ∈ F N at some time t ≥ T ⋆ N | T ⋆ N < ∞ ) ≤ P ( η t ( x ) = aa for some x ∈ ∂F N +1 at some time t ≥ T ⋆ N | T ⋆ N < ∞ ) ≤ ((2 N + 3) d − (2 N + 1) d ) × C exp( − α N ) ≤ C exp( − α N ) PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL for suitable C < ∞ and α >
0. Pick N sufficiently large so that C exp( − α N ) + C exp( − α N ) ≤ ε. Now that N is fixed, we will define a time T satisfying Lemma 5.2 as afunction of N . First of all, we observe that if site 0 is in state aa at sometime t then there is a positive probability (uniform in t ) that the cube F N is void of genes of type b at time t + 1. This implies that P ( T ⋆ N < ∞| T N = ∞ ) = 1 . Since the sequence of events { T ⋆ N ≤ n } , n ≥
0, is nondecreasing, we canapply Beppo–Levi theorem to obtain the existence of a large enough deter-ministic time
T < ∞ such that P ( T ⋆ N < T | T N = ∞ ) = 1 − ε. In conclusion, P F N ( η T ( x ) = aa for some x ∈ F N ) ≤ P F N ( T N < ∞ ) + P ( η T ( x ) = aa for some x ∈ F N ; T N = ∞ ) ≤ P F N ( T N < ∞ ) + P ( T N = ∞ ; T ⋆ N ≥ T )+ P ( η T ( x ) = aa for some x ∈ F N ; T ⋆ N < T ) ≤ P F N ( T N < ∞ ) + P ( T N = ∞ ; T ⋆ N ≥ T )+ P ( η t ( x ) = aa for some x ∈ F N at some time t ≥ T ⋆ N ; T ⋆ N < ∞ ) ≤ ε + C exp( − α N ) + C exp( − α N ) ≤ ε. This completes the proof of Lemma 5.2 and Theorem 4.
6. Proof of Theorem 6.
Theorem 6 predicts convergence to the “all aa ” configuration for a suitable choice of the parameters and for the 1-dimensional process starting from any initial configuration with infinitelymany individuals of genotype aa . In Sections 6.1 and 6.2, we will prove forthe process starting with only aa on the half line ( −∞ ,
0] that lim t →∞ r t = ∞ where r t = sup { x ∈ Z : η t ( z ) = aa for all z ≤ x } . This implies that (i) starting with a single individual of genotype aa at site0 lim t →∞ Z − t = −∞ and lim t →∞ Z + t = + ∞ , almost surely on the event { Z − t ≤ ≤ Z + t for all t ≥ } where Z − t = inf { x ≤ η t ( z ) = aa for all x ≤ z ≤ } ,Z + t = sup { x ≥ η t ( z ) = aa for all x ≥ z ≥ } , N. LANCHIER AND C. NEUHAUSER and that (ii) P ( Z − t ≤ ≤ Z + t for all t ≥ > C >
0. From (i) and (ii), itis straightforward to deduce that, with probability 1, the process startingwith infinitely many individuals of genotype aa converges to the “all aa ”configuration. Coming back to our initial objective, we assume from now onthat the process starts with only individuals of genotype aa in the interval( −∞ , t →∞ r t = ∞ , we first set K t = s t − r t − s t = inf { x ≥ r t + 1 : η t ( x ) = bb } and observe that the process K t evolves as indicated in Table 3. The proofis carried out in two steps according to whether φ aa > φ ba or φ aa ≤ φ ba .6.1. φ aa > φ ba . In this case, the process r t has a positive drift when K t = 0 (see Table 3) and the proof consists in estimating the fraction oftime K t = 0 when φ aa > φ bb + φ ba + p φ bb φ ba . First of all, we observe that, since φ aa > max( φ ba , φ bb ), the process K t jumpsfrom I to 0 at rate at most 2 φ aa , which implies that P ( K t + h = 0) ≤ P ( K t = 0) × [1 − φ aa + φ bb ) h ]+ P ( K t = 0) × φ aa h + O ( h ) . By first taking the limit as h → t → ∞ , it followsthat ddt P ( K t = 0) ≤ φ aa − (4 φ aa + 2 φ bb ) P ( K t = 0)and then that lim sup t →∞ P ( K t = 0) ≤ φ aa φ aa + φ bb . Table 3
Evolution rules of K t . In our table, I k = { k, k + 1 , . . . } , ◦ = state aa , • = state ab and • = state bb Configuration Transition Rate Drift ( r t ) ◦◦◦• · · · → I ≥ φ aa + φ bb ) ≥ − φ bb ◦◦•• · · · → φ aa + φ bb ≥ φ aa − φ ba ◦◦•• · · · → I ≥ φ ab + φ ba ≥ φ aa − φ ba ◦◦•• · · · I → φ ba / ≥ φ aa + φ ab / − φ ba ◦◦•◦ · · · I → ≤ φ aa ≥ φ aa − φ ba ◦◦•• · · · I → ≤ φ aa + φ bb + φ ab / φ ba / ≥ φ aa + φ ab / − φ ba ◦◦•◦ · · · I → ≥ φ aa − φ ba PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL Since the drift is bounded from below by φ aa − φ ba when K t = 0 (see Table3), we obtainlim inf t →∞ r t t ≥ − φ bb lim sup t →∞ P ( K t = 0) + ( φ aa − φ ba ) lim inf t →∞ P ( K t = 0) ≥ − φ bb φ aa φ aa + φ bb + ( φ aa − φ ba ) φ aa + φ bb φ aa + φ bb = [ − φ bb + ( φ aa − φ ba )(1 + φ bb /φ aa )] φ aa φ aa + φ bb = φ aa − ( φ bb + φ ba ) φ aa − φ bb φ ba φ aa + φ bb . This leads to the requirement that φ aa − ( φ bb + φ ba ) φ aa − φ bb φ ba ≥
0, whichis equivalent to φ aa ≥ ( φ bb + φ ba ) + q φ bb φ ba + ( φ bb + φ ba ) . Since, by subadditivity of the function square root, q φ bb φ ba + ( φ bb + φ ba ) ≤ p φ bb φ ba + ( φ bb + φ ba ) , this shows that condition (5) implies that lim t →∞ r t = ∞ .6.2. φ aa ≤ φ ba . The process r t may have a negative drift when K t = 1but has a positive drift when K t ≥ K t = 1 and K t ≥
2. We start by proving the result when φ aa > max (cid:18) φ ba , φ ba − φ ab (cid:19) and φ bb = 0 . In this case, the drift when K t = 0 is nonnegative (the worst case is drift =0) so we shall think of state 0 as a “neutral” state, state 1 as a “bad” stateand states in I as “good” states. We will construct a Markov process X t with state space { , , } such thatlim sup t P ( K t = 1)lim inf t P ( K t ≥ ≤ lim t P ( X t = 1)lim t P ( X t = 2) . (18)Noticing that:1. K t jumps from 1 to I at rate at least φ ab + φ ba ,2. K t jumps from I to { , } at rate at most φ aa + φ ab / φ ba ,and having in mind that 0 is neutral, 1 is a bad state and 2 is a good state,we impose that:1. X t cannot jump from 0 to 2 but jumps from 0 to 1 at some positive rate, N. LANCHIER AND C. NEUHAUSER X t jumps from 1 to 2 at rate φ ab + φ ba ,3. X t jumps from 2 to { , } at rate φ aa + φ ab / φ ba .The process X t describes a “worst case” scenario in which the time spent inbad states is rather large and the time spent in good states is rather smallcompared to K t so that (18) holds (note that the time spent in state 0 isunimportant in this inequality). Moreover, regardless of the specific valuesof the other transition rates of the process X t , the conditions above implythat the stationary distribution π = ( π , π , π ) of X t satisfies0 · π + ( φ ab + φ ba ) · π − ( φ aa + φ ab / φ ba ) · π = 0 . In particular, we can conclude thatlim sup t P ( K t = 1)lim inf t P ( K t ≥ ≤ lim sup t P ( X t = 1)lim inf t P ( X t = 2)= lim t P ( X t = 1)lim t P ( X t = 2) = π π = φ aa + φ ab / φ ba φ ab + φ ba which we rewrite as( φ ab + φ ba ) × lim sup t →∞ P ( K t = 1)(19) ≤ ( φ aa + φ ab / φ ba ) × lim inf t →∞ P ( K t ≥ . To conclude the proof, we now distinguish two cases.1. If φ aa ≤ (1 / φ ab then (19) implies thatlim sup t →∞ P ( K t = 1) ≤ lim inf t →∞ P ( K t ≥ . Moreover, the drift when K t ≥ φ aa − φ ba and the rateat which the gap process K t jumps from 0 is strictly positive, so thereexists C > t →∞ r t t ≥ C ( φ aa − φ ba ) lim sup t →∞ P ( K t = 1)+ C (4 φ aa − φ ba ) lim inf t →∞ P ( K t ≥ ≥ C (5 φ aa − φ ba ) lim sup t →∞ P ( K t = 1) > φ aa > (2 / φ ba .2. If φ aa > (1 / φ ab , then, by using the fact that φ aa ≤ φ ba together with(19), we obtain lim sup t →∞ P ( K t = 1) ≤ × lim inf t →∞ P ( K t ≥ . PATIALLY EXPLICIT NON-MENDELIAN DIPLOID MODEL Moreover, the drift when K t ≥ φ aa + (1 / φ ab − φ ba andthe rate at which the gap process K t jumps from 0 is strictly positive, sothere exists C > t →∞ r t t ≥ C ( φ aa − φ ba ) lim sup t →∞ P ( K t = 1)+ C ( φ aa + (1 / φ ab − φ ba ) lim inf t →∞ P ( K t ≥ ≥ C (3 / φ aa + (1 / φ ab − φ ba ) lim sup t →∞ P ( K t = 1) > φ aa > φ ba − (1 / φ ab .This establishes the result in the case when φ bb = 0. In addition, our proofimplies that, properly rescaled in space and time, the process can be embed-ded in a supercritical oriented percolation process (we omit the details ofthe proof) which allows us to apply a perturbation argument similar to theone in the proof of Theorem 1 and extend the result to a set of parameterswhere φ bb is small but positive. This completes the proof of Theorem 6. Acknowledgment.
The authors thank an anonymous referee for his/hernumerous suggestions to improve some of our results, and three other refereesfor additional interesting comments.REFERENCES
Bramson, M. , Ding, W. D. and
Durrett, R. (1991). Annihilating branching processes.
Stochastic Process. Appl. Bramson, M. and
Durrett, R. (1988). A simple proof of the stability criterion of Grayand Griffeath.
Probab. Theory Related Fields Bramson, M. and
Griffeath, D. (1980). On the Williams–Bjerknes tumour growthmodel. II.
Math. Proc. Cambridge Philos. Soc. Bramson, M. and
Griffeath, D. (1981). On the Williams–Bjerknes tumour growthmodel. I.
Ann. Probab. Bruck, D. (1957). Male segregation ratio advantage as a factor in maintaining lethalalleles in wild populations of house mice.
Genetics Durand, D. , Ardlie, K. , Buttel, L. , Levin, S. A. and
Silver, L. M. (1997). Impactof migration and fitness on the stability of lethal t -haplotype polymorphism in Musmusculus: A computer study. Genetics
Durrett, R. (1984). Oriented percolation in two dimensions.
Ann. Probab. Durrett, R. (1992). Multicolor particle systems with large threshold and range.
J. The-oret. Probab. Durrett, R. (1995). Ten lectures on particle systems. In
Lectures on Probability The-ory (Saint-Flour, 1993) . Lecture Notes in Mathematics
Durrett, R. and
Levin, S. (1994). The importance of being discrete (and spatial).
Theor.Popul. Biol. N. LANCHIER AND C. NEUHAUSER
Harris, T. E. (1972). Nearest-neighbor Markov interaction processes on multidimensionallattices.
Adv. in Math. Lewontin, R. C. and
Dunn, L. C. (1960). The evolutionary dynamics of a polymorphismin the house mouse.
Genetics Levin, B. , Petras, M. and
Rasmussen, D. (1969). The effect of migration on the main-tenance of a lethal polymorphism in the house mouse.
Am. Nat.
Lyttle, T. W. and
Perkins, D. (1991). Symposium: The genetics and evolutionarybiology of meiotic drive.
Am. Nat.
Neuhauser, C. (1994). A long range sexual reproduction process.
Stochastic Process.Appl. Neuhauser, C. and
Pacala, S. W. (1999). An explicitly spatial version of the Lotka–Volterra model with interspecific competition.
Ann. Appl. Probab. Nunney, L. and
Baker, E. (1993). The role of deme size, reproductive patterns anddispersal in the dynamics of t -lethal haplotypes. Evolution Petras, M. and
Topping, J. C. (1983). The maintenance of polymorphism at two lociin house mouse (Mus musculus) populations.
Can. J. Genet. Cytol. Sudbury, A. (1990). The branching annihilating process: An interacting particle system.
Ann. Probab. School of Mathematical andStatistical SciencesArizona State UniversityPO Box 871804Tempe, Arizona 85287-1804USAE-mail: [email protected]