Selection, recombination, and the ancestral initiation graph
aa r X i v : . [ q - b i o . P E ] J a n SELECTION, RECOMBINATION, AND THE ANCESTRAL INITIATIONGRAPH
FREDERIC ALBERTI, CAROLIN HERRMANN, AND ELLEN BAAKE
Abstract.
Recently, the selection-recombination equation with a single selected site and anarbitrary number of neutral sites was solved by Alberti and Baake (2020+) via the ancestralselection-recombination graph. Here, we introduce a more accessible approach, namely theancestral initiation graph, in the special case where the fitness of an individual is determinedby the leading site. The construction is based on a discretisation of the differential equation.Furthermore, we revisit the dynamics of the linkage equilibria between two neutral locihitchhiking along with a selected one, and the dependence on the position of the selectedsite. keywords: selection-recombination differential equation; ancestral selection-recombinationgraph; ancestral initiation graph; linkage disequilibria; hitchhiking; population genetics.1.
Introduction
The recombination equation is a large nonlinear dynamical system that describes the evo-lution of the distribution of genetic types within an infinite population under the influenceof recombination , that is, the reshuffling of genetic information that occurs in the process ofmeiosis during the formation of germ cells (or gametes) in sexually reproducing populations.Since its introduction by Jennings (1917), Robbins (1918), and Geiringer (1944), the recom-bination equation has posed a major challenge to mathematical population geneticists. Itwas finally solved by Baake and Baake (2016) by considering a backward-time partitioningprocess that describes the random ancestry of a single individual.The logical next step was to attack the selection -recombination equation, which describesthe evolution under the additional influence of natural selection. This was previously consid-ered unsolvable; in fact, the monograph by Akin (1979) starts with the words ‘The selection-recombination equation is a large, nonlinear system of differential equations that cannot besolved explicitly’. An approximate solution was given by Stephan et al. (2006) for the caseof a single selected locus linked to two neutral ones, with two alleles (or letters) at each ofthe three loci. The solution displays an interesting behaviour, which depends on whether theselected locus lies outside or between the neutral ones. The approximation relies on an ap-proximation of the incomplete Beta function and, while sufficiently precise, is quite involvedand does not convey any hope of generalisation.Recently, Alberti and Baake (2020+) have obtained an exact recursive solution, again for asingle selected locus but an arbitrary number of neutral loci, and an arbitrary position of theselected locus within the sequence. This solution involves intricate probabilistic constructions,
PSfrag replacements 1 2 3 4 5 61 2 3 4 5
Figure 1.
Top: Six loci within a genetic sequence. Bottom: the links be-tween the loci, numbered from 1 to 5.based on the ancestral recombination graph by Griffiths and Marjoram (1996, 1997), as wellas a generalisation of the notion of product measure. On a more abstract level, the authorsproved a formal duality between the solution of the selection-recombination equation and aso-called initiation process .The purpose of the present article is to complement this work in a number of ways. First, weassume throughout that the selected locus is the first locus in the sequence, which streamlinesthe notation considerably; we will hint at how this generalises to an arbitrary position of theselected locus as we go along. Secondly, we introduce a novel ancestral initiation graph , whichsimplifies the genealogical arguments based on the ASRG, while still retaining the underlyingintuition; it can be understood as a bridge between the ASRG and the abstract initiationprocess. We will see how both the ASRG and the ancestral initiation graph arise from adiscretisation of the selection-recombination equation. Lastly, we apply the results to theevolution of linkage disequilibria in the context of genetic hitchhiking, thus complementingthe discussions by Stephan et al. (2006) and Pfaffelhuber et al. (2008).This paper is organised as follows. We first recall the selection-recombination equation,along with the necessary concepts around it (Section 2). Then, in Section 3, we introduce theancestral selection-recombination graph, which was the central tool used by Alberti and Baake(2020+), and the ancestral initiation graph. Subsequently, in Section 4, we use the ancestralinitiation graph to give two proofs of the recursive solution in the special case consideredhere; namely, a probabilistic and a purely analytic one. Finally, in Section 5, we discuss theapplication to linkage disequilbria. The generalisation to an arbitrary position of the selectedlocus will be made via a series of remarks.2.
The selection-recombination equation and its solution
The selection-recombination equation is a system of ordinary differential equations describ-ing the evolution of the distribution of the genetic types in an infinitely large, haploid popula-tion. Equivalently, one may consider a diploid population in the absence of dominance and inHardy-Weinberg equilibrium, and work at the level of gametes. The finite set S := { , . . . , n } represents the set of genetic loci (or sites; we will use ‘locus’ and ‘sites’ synonymously) ofinterest, and L := { , . . . , n − } the set of links between them (see Fig. 1).We assume that there are two possible alleles at each locus, denoted by 0 and 1. Thus, wethink of ( genetic ) types as binary sequences of length n , i.e. elements x = ( x , . . . , x n ) of the ELECTION, RECOMBINATION, AND THE ANCESTRAL INITIATION GRAPH 3 type space X = { , } n . Since we will later also consider the evolution of subsets of loci, we de-fine, for any nonempty U ⊆ S , the corresponding marginal type x U := ( x i ) i ∈ U ∈ { , } U =: X U ,where X U is the marginal type space .We identify a population with its type distribution , a probability measure (or vector) ν = (cid:0) ν ( x ) (cid:1) x ∈ X ∈ P ( X ), where P ( X ) denotes the set of all probability measures on X ; so, ν ( { x } ) > x in the population and ν ( E ) := P x ∈ E ν ( { x } ) for E ⊆ X ; we often abbreviate ν ( { x } ) as ν ( x ). Clearly, k ν k := k ν k = P x ∈ X ν ( x ) = 1.We define the marginal distribution ν U of ν with respect to U ⊆ S via ν U ( E ) := ν ( E × X S \ U ) for all E ⊂ U. In particular, for x ∈ X U , ν U ( x ) = ν ( x, ∗ ), where ‘ ∗ ’ is a shorthand for the summation overall elements of X S \ U .We describe the evolution of the type distribution by a time-dependent family ω = (cid:0) ω t (cid:1) t > of distributions on X (so ω t ( x ) is the proportion of type x at time t ) that satisfies the selection-recombination equation (SRE)(1) ˙ ω t = Ψ sel ( ω t ) + Ψ reco ( ω t ) =: Ψ( ω t );here, the operators Ψ sel and Ψ reco describe the (independent) action of selection and recombi-nation and will be detailed below. As to selection, we assume that the fitness of an individualis determined by its allele at a single, fixed, site i ◦ ∈ S ; an individual of type x is fit if x i ◦ = 0 and unfit if x i ◦ = 1. Unless stated otherwise, we assume throughout that i ◦ = 1. Fitindividuals reproduce at rate 1 + s with s >
0, while unfit ones reproduce at rate 1, so s is theselective advantage. When an individual reproduces, the single offspring inherits the parent’stype and replaces a randomly chosen individual in the population. We write f ( ν ) := ν i ◦ (0)for the proportion of fit individuals in a population ν and b ( ν ) ( d ( ν )) for the type distribution within the subpopulation of fit (unfit) individuals . That is, for x ∈ X , b ( ν )( x ) ( d ( ν )( x )) is theprobability of sampling type x from the population ν , conditional on x i ◦ = 1. More formally, b ( ν ) and d ( ν ) can be defined via(2) f ( ν ) b ( ν )( x ) = (1 − x i ◦ ) ν ( x ) =: F ( ν )and (cid:0) − f ( ν ) (cid:1) d ( ν )( x ) := x i ◦ ν ( x ) = ν − f ( ν ) b ( ν )( x ) = (id − F )( ν ) , respectively (note that the map ν F ( ν ) is linear). The selection part of the right-hand sidein (1) therefore reads(3) Ψ sel ( ν ) = sf ( ν ) (cid:0) b ( ν ) − ν (cid:1) = s (cid:0) F ( ν ) − f ( ν ) ν (cid:1) ;note that it only reflects selective reproduction, because the ‘neutral’ reproduction (at rate 1per individual) cancels out.Recombination is modelled as follows. At rate ̺ ℓ > ℓ ∈ L , each individual is replacedby a new (offspring) sequence that is formed via a (single) crossover at link ℓ , that is, betweensites ℓ and ℓ + 1 (see Fig. 1), of two parental sequences (say, y and z ), chosen randomly from FREDERIC ALBERTI, CAROLIN HERRMANN, AND ELLEN BAAKE the population. The new individual copies the letters of the sequence y at the sites in [1 : ℓ ]and those of the sequence z at the sites in [ ℓ + 1 : n ], where for two integers a and b with a b ,we denote by [ a : b ] the discrete interval { a, a + 1 , . . . , b } . The offspring sequence thus reads( y [1: ℓ ] , z [ ℓ +1 ,n ] ). Viewed differently, the offspring sequence reads x whenever its (randomlychosen) parents are of the form y = ( x [1: ℓ ] , ∗ ) and z = ( ∗ , x [ ℓ +1 ,n ] ); this ancestral perspective(from the offspring to the parents) will become essential in what follows. Under the assump-tion of random mating, the offspring’s type is x with probability ν ( x [1 ,ℓ ] , ∗ ) ν ( ∗ , x [ ℓ +1 ,n ] ), where ν is the current type distribution. Put differently, the distribution of the offspring’s type isgiven by the product measure R ℓ ( ν ) := ν [1: ℓ ] ⊗ ν [ ℓ +1: n ] , where the operator R ℓ : P ( X ) → P ( X ) thus defined is called a recombinator (Baake and Baake2016, 2020+; see also Baake and Baake 2003). The recombination part in Eq. (1) thereforereads Ψ reco ( ν ) = n − X ℓ =1 ̺ ℓ (cid:0) R ℓ ( ω t ) − ω t (cid:1) . Note that in (Alberti and Baake 2020+), recombination rates and recombinators are indexedby sites rather than links; this is because stochastic dual processes are associated with eachsite, describing their individual ancestries. In the case i ◦ = 1, this results in an index shiftby one.In (Alberti and Baake 2020+), the solution ω of the SRE (1) was constructed iterativelyvia the solutions ω ( k ) = (cid:0) ω ( k ) t (cid:1) t > of the SRE truncated at k , for 0 k < n :(4) ˙ ω ( k ) t = Ψ ( k ) ( ω ( k ) t ) with Ψ ( k ) := Ψ sel + Ψ ( k )rec and Ψ ( k )rec := k X ℓ =1 ̺ ℓ (cid:0) R ℓ − id (cid:1) , where the initial condition is ω ( k )0 = ω for all k . For k = 0, Eq. (4) reduces to the pureselection equation ˙ ω t = Ψ sel ( ω t ); for 0 < k < n −
1, sites k + 1 , . . . , n are ‘glued together’, andfor k = n −
1, we recover Eq. (1).The following is Theorem 5.4 of Alberti and Baake (2020+) in the special case i ◦ = 1. Theorem 2.1.
For all k n − , the solutions ω ( k ) of Eq. (4) satisfy the recursion ω ( k ) t = e − ̺ k t ω ( k − t + ω ( k − k ] ,t ⊗ Z t ̺ k e − ̺ k τ ω ( k − k +1: n ] ,τ d τ, starting with the solution ω (0) t = e st F ( ω ) + (id − F )( ω )e st f ( ω ) + (1 − f ( ω )) of the pure selection equation. (cid:3) Remark 2.2.
If we assume i ◦ = n instead of i ◦ = 1, we have by symmetry that ω ( k ) t = e − ̺ ( k ) t ω ( k − t + ω ( k − n − k +1: n ] ,t ⊗ Z t ̺ ( k ) e − ̺ ( k ) τ ω ( k − n − k ] ,τ d τ, ELECTION, RECOMBINATION, AND THE ANCESTRAL INITIATION GRAPH 5 where ̺ ( k ) is the recombination rate associated with the link that connects site n − k to n − k + 1; thus effectively reversing the order of the sites. More generally, when droppingthe assumption that i ◦ = 1, we introduce, in a suitable way, new labels i , . . . , i n − for thesites (see Alberti and Baake (2020+) for details) and let ̺ ( k ) be the recombination rate of thelink that joins i k − and i k . We partition, for every 0 < k < n , the set of sites into two parts { , . . . , i k − } and { i k , . . . , n } , where C ( k ) denotes the part that contains i ◦ and D ( k ) denotesthe one that does not. The general recursion formula then reads ω ( k ) t = e − ̺ ( k ) t ω ( k − t + ω ( k − C ( k ) ,t ⊗ Z t ̺ ( k ) e − ̺ ( k ) τ ω ( k − D ( k ) ,τ d τ. ♦ That ω (0) t indeed satisfies the pure selection equation can be verified by a straightforwardcomputation. Our goal is to prove the recursion for 0 < k < n , using a stochastic rep-resentation of ω , which is related to Eq. (1) via discretisation; this approach differs fromthat taken by Alberti and Baake (2020+), who relied on an underlying approximation bystochastic models for finite populations and the corresponding duals.3. The stochastic perspective
The ancestral selection-recombination graph.
We will now provide a probabilisticinterpretation of Eq. (1) via a variant of the ancestral selection-recombination graph (ASRG) in the deterministic limit (Alberti and Baake 2020+, Sec. 7); the ASRG is a combinationof the ancestral selection graph (ASG) introduced by Krone and Neuhauser (1997) and theancestral recombination graph (ARG) (Hudson 1983; Griffiths and Marjoram 1996, 1997).The fundamental idea is to trace the (random) ancestry of an individual back in time. Here,we provide a different perspective, proceeding in two steps. First, we express the (discrete)Euler polygon (with step size h ) of the SRE (1) in terms of a random graphical construction,which turns into the ASRG in the limit h → Lemma 3.1.
Let ψ = ( ψ t ) t > be a continuous semigroup acting on a subset (Ξ , k · k ) of anormed space and let T = ( T t ) t > be a family of maps on Ξ such that ψ h = T h + O ( h ) as h → , in the sense of the uniform norm; this means that there exists a constant C > such that h sup ξ ∈ Ξ k ψ h ( ξ ) − T h ( ξ ) k C for sufficiently small h . Furthermore, assume that the family of maps ( T h − id ) /h is uniformlyLipschitz continuous. Then, (5) ψ hk = T kh + O ( h ) , as a function from ( R > , | · | ) to (Ξ , k · k ) FREDERIC ALBERTI, CAROLIN HERRMANN, AND ELLEN BAAKE uniformly in k ∈ N , where T kh denotes the k -fold composition of T h . In particular, lim h → T ⌊ th ⌋ h = ψ t , uniformly in t . Here and in what follows, asymptotic statements involving operators such as Eq. (5) arealways understood in the sense of the uniform norm.
Proof.
Let
C > L be the uniformLipschitz constant of the family of maps ( T h − id) /h . Then, for all k ∈ N > , ξ ∈ Ξ andsufficiently small h ,∆ k +1 ,h := k ψ ( k +1) h − T k +1 h k ∞ = k ψ ( k +1) h − T h ( ψ kh ) + T h ( ψ kh ) − T k +1 h k ∞ k ψ h ( ψ kh ) − T h ( ψ kh ) k ∞ + k T h ( ψ kh ) − T h ( T kh ) k ∞ Ch + Lh ∆ k,h , where ∆ ,h := 0. Thus, ∆ k,h Ch k − X j =0 ( Lh ) j Ch − Lh , and the first claim follows immediately, while the second claim follows from the continuity of ψ . (cid:3) Remark 3.2.
Assume Ξ = R and consider the general autonomous ordinary differentialequation(6) ˙ ξ = g ( ξ )with g ∈ C ( R ) with bounded first and second derivatives. Assume that its unique solution,for any initial value, exists globally. Then, let ψ be the associated flow semigroup, that is,the family of maps ( ψ t ) t > that map any ξ ∈ Ξ to the solution of Eq. (6) with initial value ξ ,evaluated at time t . Setting T h ( ξ ) = ξ + hg ( ξ ) , Lemma 3.1 reduces to the classical result that the solution of an initial value problem for thisequation can be locally uniformly approximated by its (discrete) Euler polygon; see Atkinson(1978) for background. ♦ .In our setting, the role of (Ξ , k · k ) will be played by P ( X ) (interpreted as the standardsimplex in R S ) equipped with k · k (although the precise choice will not be important), andthe semigroup ψ will be the flow associated with the SRE (1). Accordingly, we set for all h > ν ∈ P ( X )(7) T h ( ν ) := ν + h Ψ( ν ) . ELECTION, RECOMBINATION, AND THE ANCESTRAL INITIATION GRAPH 7
Note that ν Ψ( ν ) is Lipschitz continuous by virtue of being continuously differentiable onthe compact manifold P ( X ). By Taylor’s theorem, we have(8) ψ h = T h + O ( h ) . Therefore, Lemma 3.1 applies and we have(9) lim h → T ⌊ th ⌋ h ( ω ) = ω t , uniformly in t as well as in the initial condition ω .To reveal the probabilistic content of the map T h and its iterates, we write(10) T h ( ν ) = (cid:16) − hs − h n − X ℓ =1 ̺ ℓ (cid:17) ν + hs (cid:0) (1 − f ( ν )) ν + f ( ν ) b ( ν ) (cid:1) + h n − X ℓ =1 ̺ ℓ R ℓ ( ν )and assume throughout that h ( s + P n − ℓ =1 ̺ ℓ ) <
1. Then, T h ( ν ) is a convex combination of theprobability measures ν , R ( ν ) , . . . , R n − ( ν ) and b ( ν ). This means that we can sample from T h ( ν ) via the following multistage experiment.Let G h be a random graphical element that is either • a branching with probability hs , • an ℓ - splitting PSfrag replacements ℓ with probability h̺ ℓ for ℓ ∈ L , • or a straight line segment otherwise.The horizontal length of every graphical element is h , and we call the rightmost point its root and each left endpoint a leaf . After choosing a realisation of G h , we sample the type of eachleaf independently from ν . The types are subsequently propagated from left to right alongeach line, according to the following rules, which are summarised in Fig. 2.(i) When two edges meet in a branching, the type of the descendant line is the type of thebottom (or incoming ) line in case that type is fit . Otherwise, the descendant inheritsthe type of the top (or continuing ) line.(ii) When two edges meet in an ℓ -splitting, the type to the right of the square (that is,on the descendant line) is obtained by combining the ℓ leading letters of the type ofthe continuing line with the n − ℓ trailing letters of the type of the incoming line.Comparing with (10), note that in (i), the type of the descendant line indeed has thedistribution (cid:0) − f ( ν ) (cid:1) ν + f ( ν ) b ( ν ); FREDERIC ALBERTI, CAROLIN HERRMANN, AND ELLEN BAAKE
PSfrag replacements xyℓ CI D ( x , . . . , x ℓ , y ℓ +1 , . . . , y n ) x + (1 − y i ◦ )( y − x )PSfrag replacements xy ℓCI D ( x , . . . , x ℓ , y ℓ +1 , . . . , y n ) Figure 2.
The propagation of types along a realisation of G h ; in case of an ℓ -splitting, the leading ℓ loci of x are combined with the trailing n − ℓ loci of y . In a branching, the type x on the continuing line ( C ) competes against theindividual of type y on the incoming line ( I ). If y is fit, i.e if y i ◦ = 0, then thetype at the descendant line ( D ) is x + (1 − y i ◦ )( y − x ) = x + ( y − x ) = y . If y is unfit, then x + (1 − y i ◦ )( y − x ) = x and the type on the continuing lineprevails. PSfrag replacements000100001111 2 00010011 00112 h h Figure 3.
A realisation of an h -ASRG of length 2 h . To its leaves, we assignedthe types 0001 , i ◦ = 1 is marked in light brownand the top line, which does not contribute to the type of the root (i.e. is‘non-ancestral’), is dashed.with probability f ( ν ), the type on the incoming line is fit, and thus has distribution b ( ν ).With probability 1 − f ( ν ), however, it is unfit and the type on the descendant line agreeswith the type of the continuing line, which was sampled from ν . In (ii), it is clear that thetype of the descendant line has distribution R ℓ ( ν ). We summarise this in the following Lemma 3.3.
Let ν ∈ P ( X ) . The type of the root of a random graphical element G h with thetypes of its leaves sampled independently from ν has distribution T h ( ν ) of (10) . (cid:3) It is now straightforward to iterate this construction. To sample from T h ( ν ), we can drawa realisation of G h and sample the types at its leaves according to T h ( ν ). Equivalently, wecan attach to each leaf a new independent copy of G h and take instead the type delivered bythat copy, whose leaves we sample again from ν ; see Fig. 3. Definition 3.4. An h -ASRG of length kh , for k ∈ N , is a random graph constructedrecursively as follows. For k = 0, an h -ASRG of length 0 consists only of its single root, ELECTION, RECOMBINATION, AND THE ANCESTRAL INITIATION GRAPH 9
PSfrag replacements t Figure 4.
A realisation of an ASRG of length t , along with an assignment oftypes to its leaves, which are then propagated through the graph, yielding thetype 0110 of the root. The site i ◦ = 1 is marked in light brown. The lettersthat are contributed by any given line are underlined, and the lines that donot contribute any letters are dashed.which coincides with its single leaf. For k ∈ N > , an h -ASRG of length kh is constructed byattaching independent realisations of G h to each of the leaves of an h -ASRG of length ( k − h .An inductive application of Lemma 3.3 immediately shows that a sample from T kh ( ν ) canbe constructed via an h -ASRG of length k . Lemma 3.5.
Let ν ∈ P ( X ) . The type delivered by an h - ASRG of length kh with the typesat its leaves sampled independently from ν has distribution T kh ( ν ) . (cid:3) By Eq. (9), choosing k = ⌊ t/h ⌋ and letting h →
0, the distribution of the type deliveredby an h -ASRG of length kh with leaves sampled from ω converges to ω t ; at the same time(see also the sketch of proof of Theorem 3.7), these h -ASRGs converge to a limiting object,which we simply call the ASRG. Definition 3.6.
The ASRG (of length t ) is a random graph of length t that is grown fromright to left, starting with a single line. Each line is hit by branching events at rate s and,for ℓ ∈ L , it is hit by ℓ -splittings at rate ̺ ℓ . Branchings and splittings occur (conditionally)independently on each line. For an illustration, see Fig. 4.After constructing a realisation of the ASRG, we proceed as in the discrete setting. Weassign types to its leaves, sampled independently from ω , and propagate them through thegraph, from left to right, according to the rules (i) and (ii). Theorem 3.7.
The root of an
ASRG of length t with the types of its leaves sampled from ω has distribution ω t . Sketch of proof.
We have to show that the h -ASRG converges to the ASRG in the limit h → h -)ASRGand its state space, which would be needed for a rigorous proof, are beyond the scope ofthis work. We instead content ourselves with the following heuristics. Fix some γ > h ) h> of random point sets constructed by letting each t ∈ h Z be in Λ h with probability γh (assuming, of course, that h is small enough to make γh < h → Λ h = Λ , where Λ is a Poisson point process with rate γ ; for background on Poisson processes, we referthe reader to Kingman (1993). When considering any given line segment in the h -ASRG,the set of branching ( ℓ -splitting) events on that line segment can be understood as such a setΛ h with γ = s ( γ = ̺ ℓ ), restricted to that line segment. In the limit h →
0, Eq. (11) thenimplies that branching and ℓ -splitting events are laid down on each line segment according toconditionally independent Poisson point processes with rates s and ̺ ℓ . (cid:3) The ancestral initiation graph.
The ASRG can be used to prove Theorem 2.1, whichwas carried out by Alberti and Baake (2020+); however, the constructions involved are fairlysophisticated. We therefore introduce the ancestral initiation graph (AIG) here, which issimilar in spirit to the ASRG, but more convenient technically. Put simply, it describes forhow long the individuals ancestral to different sites were subjected to selection and thus buildsa bridge to the initiation process considered by Alberti and Baake (2020+).We start by discretising Eq. (1) in a slightly different way. Letting T sel ,h ( ν ) := ν + h Ψ sel ( ν ) ,T h ( ν ) can be rewritten as T h ( ν ) = (cid:16) T sel ,h ( ν ) − h n − X ℓ =1 ̺ ℓ ν (cid:17) + h n − X ℓ =1 ̺ ℓ R ℓ ( ν ) . Since h id = hT sel ,h + O ( h ) and, consequently,(12) h R ℓ ( ν ) = hν [1: ℓ ] ⊗ ν [ ℓ +1: n ] = hT sel ,h ( ν ) [1: ℓ ] ⊗ ν [ ℓ +1: n ] + O ( h )uniformly in ν , it follows that T h = e T h + O ( h ) , where(13) e T h ( ν ) := (cid:16) − h n − X ℓ =1 ̺ ℓ (cid:17) T sel ,h ( ν ) + h n − X ℓ =1 ̺ ℓ T sel ,h ( ν ) [1: ℓ ] ⊗ ν [ ℓ +1: n ] . Our motivation to rewrite T h in this way is that we want all summands to have the sameprobability f (cid:0) T sel ,h ( ν ) (cid:1) of sampling a fit individual; we will shortly see why this is important.Thus, Eq. (8) becomes(14) ψ h = e T h + O ( h ) . ELECTION, RECOMBINATION, AND THE ANCESTRAL INITIATION GRAPH 11
It is easy to check that once more, the maps ( e T h − id) /h are uniformly Lipschitz continuouswhence Lemma 3.1 applies and we havelim h → e T ⌊ th ⌋ h ( ω ) = ω t , uniformly in ω and t . As in Section 3.1, we can represent e T h graphically. For this purpose,we define another random graphical element e G h for sufficiently small h . In contrast to G h , e G h has labelled leaves, but no branchings. The law is as follows. • With probability 1 − h P n − ℓ =1 ̺ ℓ , e G h is a straight line of length h , and the label of theleaf is also h . • With probability h̺ ℓ for ℓ ∈ L , e G h is an ℓ -splitting, again of length h . The upperleaf (which contributes the sites 1 , . . . ℓ ) carries the label h ; the lower leaf (whichcontributes sites ℓ + 1 , . . . , n ) carries the label 0.Having chosen a realisation of e G h , we sample independent types for each leaf. But in contrastto Section 3.1, they are not identically distributed; rather, for a leaf with label kh , k ∈{ , } , we sample from T k sel ,h ( ν ), where T ,h := id and ν is a given probability measure,representing the initial type distribution. The types are then propagated as before. Thus,the type of the root has distribution T sel ,h ( ν ) with probability 1 − h P n − ℓ =1 ̺ ℓ and distribution T sel ,h ( ν ) [1: ℓ ] ⊗ ν [ ℓ +1: n ] with probability h̺ ℓ for ℓ ∈ L ; altogether, it thus has the desireddistribution e T h ( ν ).We now want to iterate this construction as well. To construct a sample from e T h ( ν ), wemust replace ν with e T h ( ν ). In other words, we need to sample the types of the leaves from T k sel ,h (cid:0) e T h ( ν ) (cid:1) rather than from T k sel ,h ( ν ) for k ∈ { , } (depending on the labels). We thereforeneed to understand how to construct a sample from T sel ,h (cid:0) e T h ( ν ) (cid:1) . As we will see in the proofof Proposition 3.11, all that is required for further iterations, too, is T k sel ,h (cid:0) e T h ( ν ) (cid:1) , k ∈ N , asdelivered by the following lemma. Lemma 3.8.
For all ν ∈ P ( X ) and all m ∈ N , T m sel ,h (cid:0) e T h ( ν ) (cid:1) = (cid:16) − h n − X ℓ =1 ̺ ℓ (cid:17) T m +1sel ,h ( ν ) + h n − X ℓ =1 ̺ ℓ T m +1sel ,h ( ν ) [1: ℓ ] ⊗ ν [ ℓ +1: n ] . Proof.
This result is an easy consequence of the following two elementary observations. First,for arbitrary ν and ν ′ ∈ P ( X ), we have(15) T sel ,h ( ν [1: ℓ ] ⊗ ν ′ [ ℓ +1: n ] ) = T sel ,h ( ν ) [1: ℓ ] ⊗ ν ′ [ ℓ +1: n ] . To see this, note that for all x ∈ X , F (cid:0) ν [1: ℓ ] ⊗ ν ′ [ ℓ +1: n ] (cid:1) ( x ) = (1 − x ) ν ( x , . . . , x ℓ , ∗ ) ν ′ ( ∗ , x ℓ +1 , . . . , x n )= ( F ν )( x , . . . , x ℓ , ∗ ) ν ′ ( ∗ , x ℓ +1 , . . . , x n )= (cid:0) ( F ν ) [1: ℓ ] ⊗ ν ′ [ ℓ +1] (cid:1) ( x ) and use the bilinearity of ⊗ . Second, for all ν and ν ′ with f ( ν ) = f ( ν ′ ) and all α ∈ [0 , T sel ,h (cid:0) αν + (1 − α ) ν ′ (cid:1) = αT sel ,h ( ν ) + (1 − α ) T sel ,h ( ν ′ ) . Indeed, using the last expression for Ψ sel in Eq. (3), we have T sel ,h ( ν ) := ν + hs (cid:0) F ( ν ) − f ( ν ) ν (cid:1) . Now, Eq. (16) follows from f (cid:0) αν + (1 − α ) ν ′ (cid:1) = f ( ν ) = f ( ν ′ ), together with the linearity of F .For m = 1, the claim now follows by combining (13), (16), and (15). For m >
1, we getinductively T m sel ,h (cid:0) e T ( ν ) (cid:1) = T sel ,h (cid:2) T m − ,h (cid:0) e T ( ν ) (cid:1)(cid:3) = T sel ,h h(cid:16) − h n − X ℓ =1 ̺ ℓ (cid:17) T m sel ,h ( ν ) + h n − X ℓ =1 ̺ ℓ T m sel ,h ( ν ) [1: ℓ ] ⊗ ν [ ℓ +1: n ] i = (cid:16) − n − X ℓ =1 h̺ ℓ (cid:17) T m +1sel ,h ( ν ) + h n − X ℓ =1 ̺ ℓ T m +1sel ,h ( ν ) [1: ℓ ] ⊗ ν [ ℓ +1: n ] , where we have used the induction hypothesis in the second step and (16) together with (15)in the last. (cid:3) Remark 3.9.
The idea behind this lemma is that the restriction of the nonlinear operator T sel ,h to the level sets of f is linear. This also motivated our definition of e T h in Eq. (13):because we assumed that i ◦ = 1, we chose to replace ν by T sel ,h ( ν ) in the first factor ofthe measure product, so that f agrees on both summands. More generally, if we drop theassumption i ◦ = 1, we would define instead e T h ( ν ) := (cid:16) − h n − X ℓ =1 ̺ ( ℓ ) (cid:17) T sel ,h ( ν ) + h n − X ℓ =1 ̺ ( ℓ ) T sel ,h ( ν ) C ( ℓ ) ⊗ ν D ( ℓ ) with C ( ℓ ) and D ( ℓ ) as in Remark 2.2. One would then show that T m sel ,h (cid:0) e T h ( ν ) (cid:1) = (cid:16) − h n − X ℓ =1 ̺ ( ℓ ) (cid:17) T m +1sel ,h ( ν ) + h n − X ℓ =1 ̺ ( ℓ ) T m +1sel ,h ( ν ) C ( ℓ ) ⊗ ν D ( ℓ ) , in analogy to Lemma 3.8. ♦ The following iterative construction, which is analogous to Definition 3.4, will turn out toallow to sample from e T h ( ν ); see Fig. 5 for an illustration. Definition 3.10.
For k ∈ N , an h -AIG of length kh is a random graph constructed recur-sively as follows. For k = 0, it consists only of its root, which coincides with its single leaf,with label 0. For k ∈ N > , an h -AIG of length kh is constructed recursively from an h -AIGof length ( k − h by attaching to each leaf (with label mh , where m ∈ N ) • a straight line of length h with a leaf labelled ( m + 1) h , with probability 1 − h P n − ℓ =1 ̺ ℓ ; ELECTION, RECOMBINATION, AND THE ANCESTRAL INITIATION GRAPH 13
PSfrag replacements 3 hh h h Figure 5.
A realisation of the h -AIG of length 3 h ; its leaves (and thoseof the intermediate h -ASRGs of lengths h and 2 h ) are visualised by circleswith inscribed labels. The squares correspond to splittings, and the inscribednumbers indicate the links where the splittings occur. • an ℓ -splitting, again of length h , with probability h̺ ℓ for ℓ ∈ L . The upper leaf(which contributes the sites 1 , . . . ℓ ) carries the label ( m + 1) h ; the lower leaf (whichcontributes sites ℓ + 1 , . . . , n ) carries the label 0.Note that the label of the leaf attached to the top line in an h -AIG of length kh is always kh ; for all other leaves, the label is h times the number of iterations since the line has splitoff.In view of Lemma 3.8, it now transpires that, to obtain a sample of e T kh ( ν ), the types atthe leaves with label mh must be sampled (independently) from T m sel ,h ( ν ), and propagatedforward as before. Since there are no branchings, we only have to consider rule (ii). Indeed,we obtain, in analogy with Lemma 3.5: Proposition 3.11.
When the type of each leaf with label mh of an h - AIG of length kh issampled independently from T m sel ,h ( ν ) , the type of its root has distribution e T kh ( ν ) .Proof. The claim is obvious for k = 1 from the definition of e T h . Now assume it is true for k >
0. By the induction hypothesis, we can obtain e T k +1 h ( ν ) = e T kh (cid:0) e T h ( ν ) (cid:1) , via an AIG oflength kh , replacing the sampling distribution ν for by e T h ( ν ); this means that the type ofa leaf labelled mh is distributed according to T m sel ,h (cid:0) e T h ( ν ) (cid:1) , which is the left-hand side inLemma 3.8. The right-hand side in Lemma 3.8 now tells us that we may instead performone further iteration of the AIG on the leaf of weight m and then sample the types of thedescendant leaves as claimed. Since this is true independently for every leaf of any weight, theroot type may indeed be determined via an AIG of length ( k + 1) h with sampling distribution ν . (cid:3) We finally consider, for any fixed t , the h -AIG of length ⌊ t/h ⌋ , let h →
0, and obtain thefollowing construction (compare the sketch we gave for the proof of Theorem 3.7), which isillustrated in Fig. 6.
Definition 3.12.
The AIG (of length t ) is a random graph of length t that is grown fromright to left, starting with a single line. For ℓ ∈ L , each line is hit at rate ̺ ℓ by an ℓ -splitting PSfrag replacements t t t t t t tt − t t − t t − t t − t t − t t − t Γ t ( ν ) Figure 6.
A realisation of Γ t = Γ (3) t along with the ages of its lines. At itsroot, the distribution Γ t ( ν ) of the type delivered by the AIG; here, Γ t ( ν ) = ϕ t ( ν ) { } ⊗ ϕ t − t ( ν ) { } ⊗ ϕ t − t ( ν ) { } ⊗ ϕ t − t ( ν ) { } . Note that Γ (2) t is obtained byignoring the green elements. The lines that do not carry ancestral material aredashed. For the iteration from k = 2 to k = 3, T = t − t ; see Proposition 4.1.event. Each leaf is labelled by the length (or age) θ of the line it is attached to, measuredfrom the point where that line has split off; if it has not split off from any line (as is the casefor the top line), it has length t .In the following, we will abbreviate the AIG of length t by Γ t . Instead of defining the AIGonly for a fixed length, it is natural to think of it as a Markov process Γ = (cid:0) Γ t ) t > , whichwe simply call the AIG (without specifying its length). As to the sampling, observe thatby Lemma 3.1 (see also Remark 3.2), lim h → T ⌊ θh ⌋ sel ,h = ϕ θ , uniformly in θ , where ϕ = ( ϕ t ) t > denotes the flow semigroup of the pure selection equation; in particular, ϕ t ( ω ) = ω (0) t ofTheorem 2.1. Thus, the type of a leaf of age θ must be sampled from ϕ θ ( ω ) = ω (0) θ . Tosummarise, we have the following analogue of Theorem 3.7. Theorem 3.13.
The random type delivered by Γ t with the type of each leaf sampled indepen-dently from ϕ θ ( ω ) , where θ is the age of the line the leaf is attached to, has distribution ω t . (cid:3) In the following, we write Γ t ( ω ) for the type of the root of Γ t , when the leaves are sampledas in Theorem 3.13. Then, as a corollary, we have that(17) ω t = E (cid:2) Γ t ( ω ) (cid:3) , where the expectation is taken with respect to Γ t .4. proof of Theorem 2.1 We will now use the stochastic representation of ω via the AIG to prove Theorem 2.1.In perfect analogy with the recursion in Theorem 2.1 and for all 0 k < n , we define the ELECTION, RECOMBINATION, AND THE ANCESTRAL INITIATION GRAPH 15 truncated
AIG, Γ ( k ) t , of length t by ignoring all ℓ -splittings in Γ t for ℓ > k ; note that Γ ( k ) t naturally embeds into a Markov process Γ ( k ) . See Fig. 6 for an example. In particular, Γ (0) t consists of a single straight line and Γ ( n − t = Γ t . As for the full AIG Γ t , we denote the typedelivered by Γ ( k ) t , with the types of its leaves sampled as in Theorem 3.13, by Γ ( k ) t ( ω ). Inanalogy to Eq. (17), we have(18) ω ( k ) t = E (cid:2) Γ ( k ) t ( ω ) (cid:3) . To continue, we need the notion of ancestry ; given U ⊆ S nonempty, we say that a linein Γ ( k ) t is ancestral to U if it (or rather, the leaf attached to it) contributes to the type ofthe root the letters at the sites in U . More precisely, in the absence of splitting, Γ ( k ) t consistsof a single line of length t , which is ancestral to S . As t increases, if a line that is currentlyancestral to U is hit by an ℓ -splitting, the continuing (incoming) line will be ancestral to U ∩ [1 : ℓ ] (to U ∩ [ ℓ + 1 : n ]), as long as the intersection is nonempty. If the intersectionis empty, we call that line non-ancestral ; it may be pruned without affecting the type of theroot. For an illustration, see Fig. 6.We now record the basic relation between Γ ( k ) t ( ω ) and Γ ( k − t ( ω ), which follows by de-composing Γ ( k ) t into the ancestral lines of the k leading sites and those of the n − k trailingsites; the reader should convince himself that, in Γ ( k ) t , the ancestral line of site k + 1 is alwaysancestral to k + 2 , . . . , n as well, due to the absence of ℓ -splittings for ℓ > k . Proposition 4.1.
Let t > and < k < n . Then, the following equality holds in distribution. Γ ( k ) t ( ω ) = ( Γ ( k − t ( ω ) , with probability e − ̺ k t , Γ ( k − t ( ω ) [1: k ] ⊗ e Γ ( k − T ( ω ) [ k +1: n ] , with probability − e − ̺ k t , where e Γ ( k − is an independent copy of Γ ( k − and T is an exponential random variable withmean /̺ k , conditioned to be t .Proof. It is clear that we may ignore any k -splittings that occur on lines not ancestral to site k + 1 (see Fig. 6). The reason is that, if the line is ancestral to some U k + 1, then theincoming line is ancestral to U ∩ [ k + 1 , . . . , n −
1] = ∅ and does therefore not contribute tothe type of the root.Now we distinguish two cases. First, assume that there is no k -splitting on the ancestralline of site k + 1; as k -splittings occur on every line with rate ̺ k , this happens with probabilitye − ̺ k t . Then, by ignoring the k -splittings that occur on the remaining lines, we see that inthis case, Γ ( k ) t ( ω ) = Γ ( k − t ( ω ).If, on the other hand, there is a k -splitting on the ancestral line of site k + 1, we observethat the ancestries of the first k and the trailing n − k sites are (conditionally) independent,and, again, we may ignore further k -splittings in the ancestry of the first k sites, whenceΓ ( k ) t ( ω ) [1: k ] = Γ ( k − t ( ω ) [1: k ] . Let T ′ be the time of the leftmost k -splitting on the ancestralline of site k +1. Then, starting at T ′ , this line can be identified with the corresponding line in e Γ ( k − t −T ′ as no further k -splittings are encountered. Thus, by the aforementioned independence of the ancestries of [1 : k ] and [ k + 1 : n ], we see thatΓ ( k ) t ( ω ) = Γ ( k − t ( ω ) [1: k ] ⊗ e Γ ( k − t −T ′ ( ω ) [ k +1: n ] . Finally, by the time reversal symmetry of the Poisson point process, we can replace t − T ′ by T . (cid:3) Now, as the density of T is given (as a function of τ ) by ̺ k e − ̺ k τ [0 ,t ] ( τ ) / (1 − e − ̺ k t ), wecan use Proposition 4.1 to evaluate Eq. (18) and obtain ω ( k ) t = E (cid:2) Γ ( k ) t ( ω ) (cid:3) = e − ̺ k t E (cid:2) Γ ( k − t ( ω )] + (1 − e − ̺ k t ) E (cid:2) Γ ( k − t ( ω ) [1: k ] ⊗ Γ ( k − T ( ω ) [ k +1: n ] ]= e − ̺ k t E (cid:2) Γ ( k − t ( ω ) (cid:3) + E (cid:2) Γ ( k − t ( ω ) (cid:3) [1: k ] ⊗ Z t ̺ k e − ̺ k τ E (cid:2) Γ ( k − τ ( ω ) (cid:3) [ k +1: n ] d τ = e − ̺ k t ω ( k − t + ω ( k − k ] ,t ⊗ Z t ̺ k e − ̺ k τ ω ( k − k +1: n ] ,τ d τ, thus proving Theorem 2.1. (cid:3) Remark 4.2.
Our graphical construction of the solution can easily be adapted to deal withthe case where i ◦ = 1. In a nutshell, if a line that is ancestral to the sites in a set A ⊆ S ishit by an ℓ -splitting, the continuing line (the incoming line) will be ancestral to A ∩ C ( ℓ ) (to A ∩ D ( ℓ ) ); compare Remark 2.2. We will make use of this in the next and final section. Forfurther details, we refer the reader to Alberti and Baake (2020+). ♦ A purely analytic proof.
We now condense the probabilistic argument above into apurely analytic proof of Theorem 2.1. To do this, we generalise the previous definition of theoperators T h and T sel ,h by letting, for 0 k < n and h > T ( k ) h ( ν ) := ν + h Ψ ( k ) ( ν ) , which implies that ψ ( k ) h = T ( k ) h + O ( h ). Note that T (0) h = T sel ,h and T ( n − h = T h . Thefollowing recursion is obvious: for 0 < k < n ,(19) T ( k ) h ( ν ) = T ( k − h ( ν ) + h̺ k (cid:0) R k ( ν ) − ν (cid:1) . Now, we proceed as in the construction of the AIG; observing that T ( k ) h = id + O ( h ) for all k ,we see that (19) entails the approximate recursion T ( k ) h ( ν ) = e − h̺ k T ( k − h ( ν ) + (1 − e − h̺ k ) T ( k − h ( ν ) [1: k ] ⊗ ν [ k +1: n ] + O ( h ) , where we used that e − h̺ k = 1 − h̺ k + O ( h ), together with (12) (which also holds with T sel ,h replaced by T ( k ) ). To summarise, we define a family of operators e T (0) h , . . . , e T ( n − h via e T (0) h := T sel ,h and the (exact) recursion(20) e T ( k ) h ( ν ) := e − h̺ k e T ( k − h ( ν ) + (1 − e − h̺ k ) e T ( k − h ( ν ) [1: k ] ⊗ ν [ k +1: n ] , and conclude(21) ψ ( k ) h = e T ( k ) h + O ( h ), for all 0 k < n ELECTION, RECOMBINATION, AND THE ANCESTRAL INITIATION GRAPH 17 (a straightforward induction shows that the maps ( e T ( k ) h − id) /h are indeed uniformly Lip-schitz). In analogy with the proof of Lemma 3.8, we now first establish the product structure(22) e T ( k ) h (cid:0) ν [1: k ] ⊗ ν ′ [ k +1: n ] (cid:1) = e T ( k ) h ( ν ) [1: k ] ⊗ ν ′ [ k +1: n ] . For k = 0, this is clear because it reduces to (15). For k >
0, observe that the definition of e T ( k ) h immediately gives(23) e T ( k ) h (cid:0) ν [1: k ] (cid:1) = e T ( k − h (cid:0) ν [1: k ] (cid:1) . One then shows inductively that e T ( k ) h (cid:0) ν [1: k ] ⊗ ν ′ [ k +1: n ] (cid:1) = e − h̺ k e T ( k − h (cid:0) ν [1: k ] ⊗ ν ′ [ k +1: n ] (cid:1) + (1 − e − h̺ k ) e T ( k − h (cid:2)(cid:0) ν [1: k ] ⊗ ν ′ [ k +1: n ] (cid:1) [1: k ] ⊗ ν ′ [ k +1: n ] (cid:3) = e − h̺ k e T ( k − h (cid:0) ν (cid:1) [1: k ] ⊗ ν ′ [ k +1: n ] + (1 − e − h̺ k ) e T ( k − h ( ν ) [1: k ] ⊗ ν ′ [ k +1: n ] = e T ( k ) h ( ν ) [1: k ] ⊗ ν ′ [ k +1: n ] , where we have used the definition of e T ( k ) h , (20), in the first step, the induction hypothesis inthe second, and (23) in the last.Now, let 0 < k < n be fixed. In order to evaluate (cid:0) e T ( k ) h (cid:1) m , it will be useful to extend theoperators on the right-hand side of (20) to arbitrary nonnegative measures. We define, for ν = 0, A ( ν ) := k ν k e T ( k − h (cid:16) ν k ν k (cid:17) and, for all i, j ∈ N ,(24) B i,j ( ν ) := k ν k (cid:0) e T ( k − h (cid:1) i (cid:16) ν k ν k (cid:17) [1: k ] ⊗ (cid:0) e T ( k − h (cid:1) j (cid:16) ν k ν k (cid:17) [ k +1: n ] , while setting A (0) := B i,j (0) := 0 for consistency. Obviously, k A ( ν ) k = k ν k = B i,j k ( ν ) k (sothe operators preserve the norm of nonnegative measures); furthermore, A ( αν ) = | α | A ( ν ) forarbitrary α ∈ R , and likewise for B i,j (so the operators are positive homogeneous of degreeone). It is easy to see that Eq. (22) implies the relations(25) AB i,j = B i +1 ,j , B i,j A = B i +1 ,j +1 , B i ,j B , = B i +1 , . Now we are ready to compute (cid:0) T ( k ) h (cid:1) m = (cid:0) e − h̺ k A + (1 − e − h̺ k ) B , (cid:1) m . To this end, consider C C · · · C m with C i ∈ { A, B , } for 1 i m . It is clear from (25) that C C · · · C m = ( A m , C i ≡ A,B m,j , C i = B , for at least one i, where j is the number of trailing A ’s in C C · · · C m . Thus, (cid:0) e − h̺ k A + (1 − e − h̺ k ) B , (cid:1) m = e − mh̺ k A m + m − X j =0 a j B m,j , and it remains to determine the coefficients a j . Since, for every given 0 j m − (cid:0) e − h̺ k A + (1 − e − h̺ k ) B , (cid:1) m = (cid:0) e − h̺ k A + (1 − e − h̺ k ) B , (cid:1) m − j − (cid:0) − e − h̺ k (cid:1) B , e − jh̺ k A j + Y where Y is a linear combination of terms with a number of trailing A ’s different from j , itfollows that a j = (1 − e − h̺ k )e − jh̺ k . Summarising and using our old notation, we have justshown that(26) (cid:0) e T ( k ) h (cid:1) m = e − mh̺ k (cid:0) e T ( k − h (cid:1) m + m − X j =0 (1 − e − h̺ k )e − jh̺ k (cid:0) e T ( k − h (cid:1) m [1: k ] ⊗ (cid:0) e T ( k − h (cid:1) j [ k +1: n ] . Using once more that 1 − e − h̺ k = h̺ k + O ( h ), we see that Eq. (26) implies, together withthe convergence of e T ⌊ th ⌋ h ( ω ) to ω ( k ) t as h → t , that ψ ( k ) t ( ν ) = e − ̺ k t ψ ( k − t ( ν ) + ψ ( k − t ( ν ) [1: k ] ⊗ (cid:16) lim h → h ⌊ t/h ⌋− X j =0 ̺ k e − jh̺ k ψ ( k − jh ( ν ) [ k +1: n ] (cid:17) = e − ̺ k t ω ( k − t + ω ( k − k ] ,t ⊗ Z t ̺ k e − ̺ k τ ω ( k − k +1: n ] ,τ d τ. Application: linkage disequilibria in selective sweeps
We close by showing how our results can explain the effect of a selective sweep on thecorrelation between two neutral sites. A selective sweep (Maynard Smith and Haigh 1974)occurs when a new beneficial mutation at i ◦ becomes prevalent in the population and thusalso increases the frequency of the alleles at the neutral sites that happened to be associatedwith the beneficial mutation when it arose; these alleles thus hitchhike along with the beneficialmutation. Here, we consider the simplest scenario of two neutral sites L and R that are linkedto i ◦ , and drop the assumption that i ◦ = 1. Following Stephan et al. (2006), we thereforetake S = { i ◦ , L, R } = { , , } , where i ◦ ∈ { , , } is given and L, R ∈ S \ i ◦ satisfy L < R ; L and R denote the ‘left’ and the ‘right’ neutral site, respectively, see Fig. 7. We then considerthe correlation function (or linkage disequilibrium ) between sites L and R ,(27) Cor( ω t ) := ω { L,R } ,t { (1 , } − ω { L } ,t { (1) } ω { R } ,t { (1) } . Our goal is to examine how the dynamics of the correlation is affected by the location of i ◦ relative to the neutral sites. Indeed, a somewhat complicated behaviour was observed inFig. 2 of (Stephan et al. 2006), but remained a little mysterious. A partial explanation wasgiven by Pfaffelhuber et al. (2008), which we will complement here.We are interested in a single, rare beneficial mutation that is introduced into a populationthat otherwise consists exclusively of unfit individuals. We model this by picking a singletype x mut ∈ { x ∈ X : x i ◦ = 0 } and set ω ( { x mut } ) := ε (where ε is a small positive number),together with ω ( { x } ) := 0 for all x ∈ { x ∈ X : x = x mut , x i ◦ = 0 } . More specifically, wechoose x mut L = x mut R = 1, in line with Stephan et al. (2006), and adjust the remaining typefrequencies such that ELECTION, RECOMBINATION, AND THE ANCESTRAL INITIATION GRAPH 19
PSfrag replacements ti ◦ = 1 L = 2 R = 3 L = 1 i ◦ = 2 R = 3 L = 1 R = 2 i ◦ = 3 ̺ sep = ̺ , ̺ ns = ̺ ̺ sep = ̺ + ̺ , ̺ ns = 0 ̺ sep = ̺ , ̺ ns = ̺ Figure 7.
The three cases i ◦ = 1 , i ◦ = 2, and i ◦ = 3. The selected siteis represented by a bullet, the other two (neutral) sites by circles. For eacharrangement, the values of ̺ sep and ̺ ns are given in terms of ̺ and ̺ . Notethat ̺ ns = 0 if i ◦ is in the middle. Therefore, the behaviour is fundamentallydifferent in this case compared to when i ◦ is one of the outer sites. See alsoFig. 8 for a comparison of the time-evolution of the correlation function in thecase ̺ ns = 0 versus ̺ sep = 0. ,
000 4 ,
000 6 , . . ̺ ns = 0 ̺ ns = 0 . s̺ ns = 0 . s̺ ns = 0 . st C o r ( ω n s t ) (a) ̺ sep = 0 ,
000 4 ,
000 6 , . . ̺ ns = 0 ̺ ns = 0 . s̺ ns = 0 . s̺ ns = 0 . st C o r ( ω t ) (b) ̺ sep = 0 . s Figure 8.
Time evolution of the correlation under recombination and selec-tion obtained by evaluating the solution formula from Theorem 2.1. In theleft panel, recombination only separates the block { L, R } from the selectedsite, but not L and R from each other. In the right panel, separating recom-bination is added. Parameters: s = 10 − , i ◦ = 1, ω ( { (0 , , } ) = 5 · − = ε, ω ( { (1 , , } ) = 0 . , ω ( { (1 , , } ) = 0 . , ω ( { (1 , , } ) = 0 . ω ( { (1 , , } ) = 0 . • Cor( ω ) >
0, and • for ̺ = ̺ = 0, one has dd t Cor( ω t ) | t =0 > x L , x R ) = (1 ,
1) hitchhikes alongwith x i ◦ = 0).The exact parameter values are given in Fig. 8.It is clear that, for ̺ = ̺ = 0, the correlation eventually decays to zero, after the initialincrease due to hitchhiking; compare Fig. 8 (a). This is because, under the pure selection equation, the single fit mutant type ultimately goes to fixation, that is,(28) ω (0) ∞ := lim t →∞ ω (0) t = δ x mut , and the correlation vanishes for any point measure. Let us now investigate how this behaviourchanges in the presence of recombination. To this end, it is essential to distinguish betweenrecombination events that separate L and R and those that do not. We therefore define ̺ sep as the total recombination rate between the sites L and R , and ̺ ns as the remainingrecombination rate, that is, the recombination rate that keeps { L, R } intact but separates itfrom i ◦ , see Fig. 7: ̺ sep := X ℓ ∈L : L ∈ [1: ℓ ] ,R ∈ [ ℓ +1: n ] ̺ ℓ and ̺ ns := (cid:16) n − X ℓ =1 ̺ ℓ (cid:17) − ̺ sep . (29)We will then see that the dynamics of the linkage disequilibrium (27) may be reparametrisedin terms of ̺ sep and ̺ ns , thus eliminating the need to distinguish the different locations of i ◦ .More precisely, the position of i ◦ affects the dynamics of the correlation only via ̺ sep and ̺ ns ;in particular, i ◦ = 2 implies ̺ ns = 0.Let us first consider ω ns , the solution (with initial condition ω ) of Eq. (1) with all ̺ ℓ forwhich L ∈ [1 : ℓ ] and R ∈ [ ℓ + 1 : n ] (and therefore ̺ sep ) set to 0 (note that, for i ◦ = 2, ̺ ns = 0and therefore ω ns = ω (0) ). We then have Lemma 5.1.
For ̺ ns = 0 , one has Cor( ω ns ∞ ) = ω (0) ∞ = 0 . For ̺ ns > , Cor( ω ns ∞ ) = Z ∞ ̺ ns e − ̺ ns τ ω (0) { L,R } ,τ (1 ,
1) d τ − Z ∞ ̺ ns e − ̺ ns τ ω (0) { L } ,τ (1) d τ Z ∞ ̺ ns e − ̺ ns σ ω (0) { R } ,σ (1) d σ. Proof.
For ̺ ns = 0, one has ω ns = ω (0) and the case is clear due to (28). For ̺ ns > i ◦ = 1, the claim follows immediately from Theorem 2.1 by letting t → ∞ and marginalisation;this carries over to i ◦ = 3 via symmetry. Alternatively, we can argue via the AIG: as t tendsto infinity, the age of the line ancestral to { L, R } , which is never split since ̺ sep = 0, isexponentially distributed with mean 1 /̺ ns ; hence the type of this line is sampled from ω (0) { L,R } evaluated at this exponential time. (cid:3) The numerical solution of (cid:0)
Cor( ω ns t ) (cid:1) t > is shown in Fig. 8 (a) for various values of ̺ ns . Byand large, and as noted by Pfaffelhuber et al. (2008), the correlation is expected to remainpositive in the long run for ̺ ns > x L , x R ) = (0 , , (0 , , (1 , x mut L , x mut R ) = (1 ,
1) in the mutant and thwart its estab-lishment as it is ‘swept’ through the population together with x i ◦ = 0. Consequently, theasymptotic distribution is not a point measure, whence we expect a non-zero limit of thecorrelation as t → ∞ . In more detail, { L, R } initially travels together with i ◦ , so ω ns t is close ELECTION, RECOMBINATION, AND THE ANCESTRAL INITIATION GRAPH 21 to ω (0) t for small t , and likewise for the correlations; note that ω (0) t equals ω ns t with ̺ ns = 0.As soon as { L, R } separates from i ◦ , which happens the sooner the larger ̺ ns , ω ns t deviatesfrom ω (0) t ; in particular, Cor( ω ns t ) ‘decouples’ from Cor( ω (0) t ) and then remains constant since ̺ sep = 0. For ̺ ns = 0 . s , this happens almost instantly, so the correlation is nearly con-stant for all times; this is in agreement with in Lemma 5.1, which tells us that, as ̺ ns → ∞ ,Cor( ω ns ∞ ) approaches Cor( ω ). To the best of our knowledge, the details of this dynamics havenot been described previously. Let us now add the action of ̺ sep . Proposition 5.2.
Let ̺ sep be as in (29) . Then, Cor( ω t ) = e − ̺ sep t Cor( ω ns t ) . Proof.
We argue via the AIG. The line ancestral to { L, R } is hit at rate ̺ sep by a splittingevent that separates the ancestral lines of L and R . Thus, with probability e − ̺ sep t , no suchsplitting occurs on the ancestral line of L and R , and thus the type agrees with the onedelivered by an AIG with ̺ sep = 0; its distribution is ω ns t . On the other hand, if such asplitting has ocurred, the letters at sites L and R are sampled independently, and are thusuncorrelated. (cid:3) The behaviour is shown in Fig. 8 (b) for the values of ̺ ns used in panel (a) — in line withProposition 5.2, it is obtained by multiplying the functions in panel (a) with an exponentiallydecaying factor. The resulting picture resembles Fig. 2 of Stephan et al. (2006). Acknowledgements
This work was supported by the German Research Foundation (DFG), within CRC 1283(project C1).
References
E. Akin,
The Geometry of Population Genetics , Springer, New York (1979).F. Alberti and E. Baake, Solving the selection-recombination equation: Ancestral lines underselection and recombination; arXiv:2003.06831
K. E. Atkinson,
An Introduction to Numerical Analysis , Wiley, New York (1978).E. Baake and M. Baake, An exactly solved model for recombination, mutation and selection,
Can. J. Math. (2003), 3–41; and erratum Can. J. Math. (2008), 264.E. Baake and M. Baake, Haldane linearisation done right: Solving the nonlinear recombinationequation the easy way, Discr. Cont. Dyn. Syst. A (2016), 6645–6656.E. Baake and M. Baake, Ancestral lines under recombination, in: Probabilistic Struc-tures in Evolution , E. Baake and A. Wakolbinger (eds.), EMS Press, Berlin, in press; arXiv:2002.08658 .S. N. Ethier and T. G. Kurtz,
Markov Processes: Characterization and Convergence , Wiley,New York (1986; reprint 2005).
H. Geiringer, On the probability theory of linkage in Mendelian heredity,
Ann. Math. Stat. (1944), 25–57.R.C. Griffiths and P. Marjoram, Ancestral inference from samples of DNA sequences withrecombination, J. Comput. Biol. (1996), 479–502.R. C. Griffiths and P. Marjoram, An ancestral recombination graph. in: Progress in PopulationGenetics and Human Evolution , P. Donnelly and S. Tavar´e (eds.), Springer, New York(1997), pp. 257–270.R.R. Hudson, Properties of an neutral allele model with intragenetic recombination,
Theor.Popul. Biol. (1983), 183–201.H.S. Jennings, The numerical results of diverse systems of breeding, with respect to two pairsof characters, linked or independent, with special relation to the effects of linkage, Genetics (1917), 97–154.J. F. C. Kingman, Poisson Processes , Clarendon Press, Oxford (1993; reprint 2010).S.M. Krone and C. Neuhauser, Ancestral processes with selection,
Theor. Popul. Biol. (1997), 210–237.T. Kurtz, Solutions of ordinary differential equations as limits of pure jump Markov processes, J. Appl. Probab. (1970), 49–58.J. Maynard Smith and J. Haigh, The hitch-hiking effect of a favourable gene, Genet. Res. (1974), 23–35.P. Pfaffelhuber, A. Lehnert and W. Stephan, Linkage disequilibrium under genetic hitchhikingin finite population, Genetics (2008), 527–537.R.B. Robbins, Some applications of mathematics to breeding problems III,
Genetics (1918),375–389.W. Stephan, Y. S. Song and C. H. Langley, The hitchhiking effect on linkage disequilibriumbetween linked neutral loci, Genetics (2006), 2647–2663. { Faculty of Mathematics, Faculty of Technology } , Bielefeld University, Postbox 100131,33501 Bielefeld, Germany Email address : { falberti,ebaake } @math.uni-bielefeld.de Institute of Biometry and Clinical Epidemiology, Charit´e - Universit¨atsmedizin Berlin, Charit´eplatz1, 10117 Berlin, Germany
Email address ::