Two-locus coalescent

InductiveStates of the two-locus ancestral process#✓fully proved

The states of the two-locus ancestral process of a sample of two sequences.

inductive TwoLocusState
  /-- Two chromosomes, each ancestral at both loci. -/
  | linked
  /-- One chromosome ancestral at both loci, one at `A` only and one at `B` only. -/
  | partlySplit
  /-- Four chromosomes: two ancestral at `A` only and two at `B` only. -/
  | split
  /-- Locus `A` has coalesced; two lineages remain at locus `B`. -/
  | aCoalesced
  /-- Locus `B` has coalesced; two lineages remain at locus `A`. -/
  | bCoalesced
  /-- Both loci have coalesced. -/
  | bothCoalesced
  deriving DecidableEq
Used by 17

TheoremSum over the two-locus ancestral states#✓fully proved

A sum over the six states.

theorem sum_univ {M : Type*} [AddCommMonoid M] (f : TwoLocusState → M) :
    ∑ i, f i = f linked + f partlySplit + f split + f aCoalesced + f bCoalesced +
      f bothCoalesced
Uses 1 · Used by 6

DefinitionCoalescence indicator at locus A#✓fully proved

Whether locus A has found its most recent common ancestor.

def coalescedA : TwoLocusState → Bool
  | aCoalesced | bothCoalesced => true
  | _ => false
Uses 1 · Used by 9

DefinitionCoalescence indicator at locus B#✓fully proved

Whether locus B has found its most recent common ancestor.

def coalescedB : TwoLocusState → Bool
  | bCoalesced | bothCoalesced => true
  | _ => false
Uses 1 · Used by 8

DefinitionRate matrix of the two-locus coalescent#✓fully proved

The rate matrix of the two-locus ancestral process of a sample of two, with time in units of 2 N generations and ρ = 4 N r: coalescence at rate 1 for every pair of ancestral lineages, recombination at rate ρ / 2 for every lineage ancestral at both loci.

MathlibExt.Biology.PopulationGenetics.TwoLocusCoalescent.twoLocusRate · 19 lines
noncomputable def twoLocusRate (ρ : ℝ) : Matrix TwoLocusState TwoLocusState ℝ := fun i j =>
  match i, j with
  | linked, linked => -(1 + ρ)
  | linked, partlySplit => ρ
  | linked, bothCoalesced => 1
  | partlySplit, linked => 1
  | partlySplit, partlySplit => -(3 + ρ / 2)
  | partlySplit, split => ρ / 2
  | partlySplit, aCoalesced => 1
  | partlySplit, bCoalesced => 1
  | split, partlySplit => 4
  | split, split => -6
  | split, aCoalesced => 1
  | split, bCoalesced => 1
  | aCoalesced, aCoalesced => -1
  | aCoalesced, bothCoalesced => 1
  | bCoalesced, bCoalesced => -1
  | bCoalesced, bothCoalesced => 1
  | _, _ => 0
Uses 1 · Used by 10

TheoremThe two-locus coalescent rate matrix is a rate matrix#✓fully proved

For ρ ≥ 0, twoLocusRate ρ is a rate matrix.

theorem isRateMatrix_twoLocusRate {ρ : ℝ} (hρ : 0 ≤ ρ) : IsRateMatrix (twoLocusRate ρ)
Uses 4 · Used by 4

TheoremPositivity of the two-locus moment denominator#✓fully proved

The denominator ρ² + 13 ρ + 18 is positive for ρ ≥ 0.

theorem twoLocusDenom_pos {ρ : ℝ} (hρ : 0 ≤ ρ) : 0 < twoLocusDenom ρ
Uses 1 · Used by 10

TheoremMean coalescence time at locus A#✓fully proved

For ρ ≥ 0, the backward equations of E_i T_A (m i = 0 where A has coalesced, (Q m) i = -1 elsewhere) force m i = 1 in every state where A has not coalesced: locus A alone is Kingman's coalescent of two genes.

theorem eq_one_of_meanTimeA {ρ : ℝ} (hρ : 0 ≤ ρ) {m : TwoLocusState → ℝ}
    (h0 : ∀ i, coalescedA i = true → m i = 0)
    (h : ∀ i, coalescedA i = false → (twoLocusRate ρ *ᵥ m) i = -1) :
    ∀ i, coalescedA i = false → m i = 1
Uses 6 · Used by 4

TheoremMean coalescence time at locus B#✓fully proved

Mean coalescence time at locus B (the mirror image of eq_one_of_meanTimeA).

theorem eq_one_of_meanTimeB {ρ : ℝ} (hρ : 0 ≤ ρ) {m : TwoLocusState → ℝ}
    (h0 : ∀ i, coalescedB i = true → m i = 0)
    (h : ∀ i, coalescedB i = false → (twoLocusRate ρ *ᵥ m) i = -1) :
    ∀ i, coalescedB i = false → m i = 1
Uses 6 · Used by 3

TheoremSecond moment of T_A#✓fully proved

If m solves the mean equations of eq_one_of_meanTimeA and g the second-moment equations (g i = 0 where A has coalesced, (Q g) i = -2 m i elsewhere), then g i = 2 where A has not coalesced, so the variance read off from these equations is 2 - 1² = 1.

theorem eq_two_of_secondMomentA {ρ : ℝ} (hρ : 0 ≤ ρ) {m g : TwoLocusState → ℝ}
    (hm0 : ∀ i, coalescedA i = true → m i = 0)
    (hm : ∀ i, coalescedA i = false → (twoLocusRate ρ *ᵥ m) i = -1)
    (hg0 : ∀ i, coalescedA i = true → g i = 0)
    (hg : ∀ i, coalescedA i = false → (twoLocusRate ρ *ᵥ g) i = -2 * m i) :
    ∀ i, coalescedA i = false → g i = 2
Uses 7

DefinitionCovariance of coalescence times at two linked loci#✓fully proved

The covariance Cov(T_A, T_B) = (ρ + 18) / (ρ² + 13 ρ + 18) of the coalescence times at two loci of two sequences sampled linked (Hudson 1983; Kaplan & Hudson 1985).

noncomputable def twoLocusCov (ρ : ℝ) : ℝ := (ρ + 18) / twoLocusDenom ρ
Uses 1 · Used by 10

TheoremCovariance of coalescence times#✓fully proved

Covariance of coalescence times (Hudson 1983; Kaplan & Hudson 1985). Let mA, mB solve the mean equations of the two loci and f the product-moment equations (f i = 0 where either locus has coalesced, (Q f) i = -(mA i + mB i) where neither has). Then f - mA mB is (ρ + 18) / (ρ² + 13ρ + 18) from linked, 6 / (ρ² + 13ρ + 18) from partlySplit and 4 / (ρ² + 13ρ + 18) from split. This is a statement about the backward equations; that their solutions are the path-space moments of the chain is proved in TwoLocusCoalescentPath.

theorem covariance_of_productMoment {ρ : ℝ} (hρ : 0 ≤ ρ) {mA mB f : TwoLocusState → ℝ}
    (hA0 : ∀ i, coalescedA i = true → mA i = 0)
    (hA : ∀ i, coalescedA i = false → (twoLocusRate ρ *ᵥ mA) i = -1)
    (hB0 : ∀ i, coalescedB i = true → mB i = 0)
    (hB : ∀ i, coalescedB i = false → (twoLocusRate ρ *ᵥ mB) i = -1)
    (hf0 : ∀ i, (coalescedA i || coalescedB i) = true → f i = 0)
    (hf : ∀ i, (coalescedA i || coalescedB i) = false →
      (twoLocusRate ρ *ᵥ f) i = -(mA i + mB i)) :
    f linked - mA linked * mB linked = twoLocusCov ρ ∧
      f partlySplit - mA partlySplit * mB partlySplit = 6 / twoLocusDenom ρ ∧
      f split - mA split * mB split = 4 / twoLocusDenom ρ
Uses 10 · Used by 1

DefinitionSolution of the two-locus product-moment equations#✓fully proved

The explicit solution of the product-moment equations: 1 + Cov in the three states where neither locus has coalesced, 0 elsewhere.

noncomputable def productMoment (ρ : ℝ) : TwoLocusState → ℝ
  | linked => 1 + twoLocusCov ρ
  | partlySplit => 1 + 6 / twoLocusDenom ρ
  | split => 1 + 4 / twoLocusDenom ρ
  | _ => 0
Uses 3 · Used by 4

DefinitionMean coalescence time at locus A#✓fully proved

The mean coalescence time at locus A: 1 before A coalesces, 0 after.

def meanTimeA (i : TwoLocusState) : ℝ := if coalescedA i then 0 else 1
Uses 2 · Used by 4

DefinitionMean coalescence time at locus B#✓fully proved

The mean coalescence time at locus B: 1 before B coalesces, 0 after.

def meanTimeB (i : TwoLocusState) : ℝ := if coalescedB i then 0 else 1
Uses 2 · Used by 4

TheoremThe explicit two-locus moments solve the moment equations#✓fully proved

The systems of covariance_of_productMoment are solvable: meanTimeA, meanTimeB and productMoment ρ satisfy them.

theorem productMoment_solution {ρ : ℝ} (hρ : 0 ≤ ρ) :
    (∀ i, coalescedA i = false → (twoLocusRate ρ *ᵥ meanTimeA) i = -1) ∧
    (∀ i, coalescedB i = false → (twoLocusRate ρ *ᵥ meanTimeB) i = -1) ∧
    (∀ i, (coalescedA i || coalescedB i) = true → productMoment ρ i = 0) ∧
    (∀ i, (coalescedA i || coalescedB i) = false →
      (twoLocusRate ρ *ᵥ productMoment ρ) i = -(meanTimeA i + meanTimeB i))
Uses 11 · Used by 3

TheoremTwo-locus covariance under complete linkage#✓fully proved

Complete linkage: Cov(T_A, T_B) = 1 = Var T_A.

theorem twoLocusCov_zero : twoLocusCov 0 = 1
Uses 2

TheoremPositivity of the two-locus covariance#✓fully proved

The covariance is positive for ρ ≥ 0.

theorem twoLocusCov_pos {ρ : ℝ} (hρ : 0 ≤ ρ) : 0 < twoLocusCov ρ
Uses 3 · Used by 3

TheoremTwo-locus covariance is at most one#✓fully proved

The covariance is at most Var T_A = 1 for ρ ≥ 0 (outside this domain the formula can exceed 1; e.g. twoLocusCov (-1) = 17 / 6).

theorem twoLocusCov_le_one {ρ : ℝ} (hρ : 0 ≤ ρ) : twoLocusCov ρ ≤ 1
Uses 3

TheoremTwo-locus covariance decreases with recombination rate#✓fully proved

The covariance decreases with the recombination rate.

theorem antitoneOn_twoLocusCov : AntitoneOn twoLocusCov (Set.Ici 0)
Uses 3

TheoremTwo-locus covariance vanishes under free recombination#✓fully proved

Free recombination: Cov(T_A, T_B) → 0 as ρ → ∞.

theorem tendsto_twoLocusCov_atTop : Tendsto twoLocusCov atTop (𝓝 0)
Uses 4

DefinitionLinkage disequilibrium polynomial#✓fully proved

The linkage disequilibrium D = x_AB x_ab - x_Ab x_aB as a polynomial.

noncomputable def diseqPoly : MvPolynomial (Bool × Bool) ℝ :=
  X (true, true) * X (false, false) - X (true, false) * X (false, true)
Used by 7

DefinitionFrequency of allele A as a polynomial#✓fully proved

The frequency p = x_AB + x_Ab of allele A as a polynomial.

noncomputable def freqFstPoly : MvPolynomial (Bool × Bool) ℝ := X (true, true) + X (true, false)
Used by 2

DefinitionFrequency of allele B as a polynomial#✓fully proved

The frequency q = x_AB + x_aB of allele B as a polynomial.

noncomputable def freqSndPoly : MvPolynomial (Bool × Bool) ℝ := X (true, true) + X (false, true)
Used by 2

DefinitionWeighted disequilibrium polynomial \(D(1-2p)(1-2q)\)#✓fully proved

D (1 - 2p)(1 - 2q) as a polynomial.

noncomputable def diseqZPoly : MvPolynomial (Bool × Bool) ℝ :=
  diseqPoly * (1 - 2 * freqFstPoly) * (1 - 2 * freqSndPoly)
Uses 3 · Used by 2

DefinitionProduct of heterozygosities polynomial#✓fully proved

Π = p (1 - p) q (1 - q) as a polynomial.

noncomputable def hetProdPoly : MvPolynomial (Bool × Bool) ℝ :=
  freqFstPoly * (1 - freqFstPoly) * freqSndPoly * (1 - freqSndPoly)
Uses 2 · Used by 1

DefinitionGenerator of the two-locus Wright–Fisher diffusion#✓fully proved

The generator of the two-locus Wright–Fisher diffusion on polynomials in the gamete frequencies, time in units of 2 N generations, ρ = 4 N r, symmetric mutation at rate ν per locus in each direction (ν = θ / 2, θ = 4 N u): L f = ½ ∑ x_g (δ_{gh} - x_h) ∂_g ∂_h f + ∑ (-(ρ/2) ε_g D + ν (x_{g'} + x_{g''} - 2 x_g)) ∂_g f.

noncomputable def twoLocusGenerator (ρ ν : ℝ) (f : MvPolynomial (Bool × Bool) ℝ) :
    MvPolynomial (Bool × Bool) ℝ :=
  ∑ g, ∑ h, C (1 / 2 : ℝ) * X g * ((if g = h then 1 else 0) - X h) * pderiv g (pderiv h f) +
    ∑ g, (C (-(ρ / 2) * (if g.1 = g.2 then 1 else -1)) * diseqPoly +
      C ν * (X (!g.1, g.2) + X (g.1, !g.2) - 2 * X g)) * pderiv g f
Uses 1 · Used by 6

TheoremDrift of gamete frequencies matches recombination#✓fully proved

Without mutation the drift of a gamete frequency is 2 N times its per-generation change under Linkage.recombine when ρ = 4 N r.

theorem twoLocusGenerator_X {x : Bool × Bool → ℝ} (hx : ∑ g, x g = 1) {N r ρ : ℝ}
    (hρ : ρ = 4 * N * r) (g : Bool × Bool) :
    eval x (twoLocusGenerator ρ 0 (X g)) = 2 * N * (recombine r x g - x g)
Uses 5

TheoremHill–Robertson / Ohta–Kimura, first row#✓fully proved

On the simplex, L D² = -(3 + ρ + 8ν) D² + D (1 - 2p)(1 - 2q) + p (1 - p) q (1 - q) (for all real ρ, ν).

theorem twoLocusGenerator_diseqSq {x : Bool × Bool → ℝ} (hx : ∑ g, x g = 1) (ρ ν : ℝ) :
    eval x (twoLocusGenerator ρ ν (diseqPoly ^ 2)) =
      -(3 + ρ + 8 * ν) * linkageDiseq x ^ 2 +
        linkageDiseq x * (1 - 2 * alleleFreqFst x true) * (1 - 2 * alleleFreqSnd x true) +
        alleleFreqFst x true * (1 - alleleFreqFst x true) * alleleFreqSnd x true *
          (1 - alleleFreqSnd x true)
Uses 5 · Used by 1

TheoremHill–Robertson / Ohta–Kimura, second row#✓fully proved

On the simplex, L (D z) = 4 D² - (5 + ρ/2 + 8ν) D z with z = (1 - 2p)(1 - 2q).

theorem twoLocusGenerator_diseqZ {x : Bool × Bool → ℝ} (hx : ∑ g, x g = 1) (ρ ν : ℝ) :
    eval x (twoLocusGenerator ρ ν diseqZPoly) =
      4 * linkageDiseq x ^ 2 - (5 + ρ / 2 + 8 * ν) *
        (linkageDiseq x * (1 - 2 * alleleFreqFst x true) * (1 - 2 * alleleFreqSnd x true))
Uses 6 · Used by 1

TheoremHill–Robertson / Ohta–Kimura, third row#✓fully proved

On the simplex, L Π = D z - (2 + 8ν) Π + ν (p (1 - p) + q (1 - q)); for ν = 0 the three rows close.

theorem twoLocusGenerator_hetProd {x : Bool × Bool → ℝ} (hx : ∑ g, x g = 1) (ρ ν : ℝ) :
    eval x (twoLocusGenerator ρ ν hetProdPoly) =
      linkageDiseq x * (1 - 2 * alleleFreqFst x true) * (1 - 2 * alleleFreqSnd x true) -
        (2 + 8 * ν) * (alleleFreqFst x true * (1 - alleleFreqFst x true) * alleleFreqSnd x true *
          (1 - alleleFreqSnd x true)) +
        ν * (alleleFreqFst x true * (1 - alleleFreqFst x true) +
          alleleFreqSnd x true * (1 - alleleFreqSnd x true))
Uses 6

StructureStationary law of the two-locus diffusion#✓fully proved

A stationary law of the two-locus diffusion in the weak (basic adjoint) sense: a probability measure on the gamete-frequency simplex with ∫ L f dμ = 0 for every polynomial f.

structure IsStationaryLaw (ρ ν : ℝ) (μ : Measure (Bool × Bool → ℝ)) : Prop where
  /-- `μ` is a probability measure. -/
  isProbabilityMeasure : IsProbabilityMeasure μ
  /-- `μ` lives on the simplex of gamete frequencies. -/
  ae_mem_simplex : ∀ᵐ x ∂μ, (∀ g, 0 ≤ x g) ∧ ∑ g, x g = 1
  /-- The expected generator vanishes on polynomials. -/
  integral_generator : ∀ f, ∫ x, eval x (twoLocusGenerator ρ ν f) ∂μ = 0
Uses 1 · Used by 1

DefinitionStationary standardized linkage disequilibrium \(\sigma_d^2\)#✓fully proved

The stationary standardized disequilibrium σ_d² = (10 + ρ + 16ν) / (22 + 13ρ + ρ² + 128ν + 24ρν + 128ν²).

noncomputable def sigmaDSq (ρ ν : ℝ) : ℝ :=
  (10 + ρ + 16 * ν) / (22 + 13 * ρ + ρ ^ 2 + 128 * ν + 24 * ρ * ν + 128 * ν ^ 2)
Used by 3

TheoremOhta–Kimura σ_d²#✓fully proved

For a stationary law with ρ, ν ≥ 0 and E Π > 0, E D² / E[p (1 - p) q (1 - q)] = sigmaDSq ρ ν.

theorem IsStationaryLaw.sigmaDSq_eq {ρ ν : ℝ} (hρ : 0 ≤ ρ) (hν : 0 ≤ ν)
    {μ : Measure (Bool × Bool → ℝ)} (hμ : IsStationaryLaw ρ ν μ)
    (hpos : 0 < ∫ x, alleleFreqFst x true * (1 - alleleFreqFst x true) * alleleFreqSnd x true *
      (1 - alleleFreqSnd x true) ∂μ) :
    (∫ x, linkageDiseq x ^ 2 ∂μ) / (∫ x, alleleFreqFst x true * (1 - alleleFreqFst x true) *
      alleleFreqSnd x true * (1 - alleleFreqSnd x true) ∂μ) = sigmaDSq ρ ν
Uses 10

TheoremOhta–Kimura formula for \(\sigma_d^2\) without mutation#✓fully proved

The low-mutation value σ_d² = (10 + ρ) / (22 + 13ρ + ρ²) (Ohta & Kimura 1971).

theorem sigmaDSq_zero (ρ : ℝ) : sigmaDSq ρ 0 = (10 + ρ) / (22 + 13 * ρ + ρ ^ 2)
Uses 1 · Used by 1

TheoremLow-mutation limit of \(\sigma_d^2\)#✓fully proved

σ_d² → (10 + ρ) / (22 + 13ρ + ρ²) as the mutation rate tends to 0, for ρ ≥ 0.

theorem tendsto_sigmaDSq_zero {ρ : ℝ} (hρ : 0 ≤ ρ) :
    Tendsto (sigmaDSq ρ) (𝓝 0) (𝓝 ((10 + ρ) / (22 + 13 * ρ + ρ ^ 2)))
Uses 2

TheoremMean coalescence time at locus A, path form#✓fully proved

For the two-locus ancestral process started from two linked sequences, the hitting time T_A of the states where locus A has coalesced has mean 1.

theorem lintegral_coalescenceTimeA {ρ : ℝ} (hρ : 0 ≤ ρ) :
    ∫⁻ ω, (⨅ (t : ℝ≥0) (_ : coalescedA (ctmcPath (twoLocusRate ρ) linked ω t) = true),
        (t : ℝ≥0∞)) ∂ctmcInputMeasure (isRateMatrix_twoLocusRate hρ) = 1
Uses 14 · Used by 1

TheoremMean coalescence time at locus B, path form#✓fully proved

From two linked sequences, the hitting time T_B of the states where locus B has coalesced has mean 1.

theorem lintegral_coalescenceTimeB {ρ : ℝ} (hρ : 0 ≤ ρ) :
    ∫⁻ ω, (⨅ (t : ℝ≥0) (_ : coalescedB (ctmcPath (twoLocusRate ρ) linked ω t) = true),
        (t : ℝ≥0∞)) ∂ctmcInputMeasure (isRateMatrix_twoLocusRate hρ) = 1
Uses 14 · Used by 1

TheoremProduct moment of the two coalescence times, path form#✓fully proved

Product moment of the two coalescence times, path form (Hudson 1983; Kaplan & Hudson 1985): from two linked sequences, E (T_A T_B) = 1 + (ρ + 18) / (ρ² + 13ρ + 18).

theorem lintegral_coalescenceTimeA_mul_coalescenceTimeB {ρ : ℝ} (hρ : 0 ≤ ρ) :
    ∫⁻ ω, (⨅ (t : ℝ≥0) (_ : coalescedA (ctmcPath (twoLocusRate ρ) linked ω t) = true),
        (t : ℝ≥0∞)) *
      (⨅ (t : ℝ≥0) (_ : coalescedB (ctmcPath (twoLocusRate ρ) linked ω t) = true),
        (t : ℝ≥0∞)) ∂ctmcInputMeasure (isRateMatrix_twoLocusRate hρ) =
      ENNReal.ofReal (1 + twoLocusCov ρ)
Uses 21 · Used by 1

TheoremCovariance of the coalescence times, path form#✓fully proved

Covariance of the coalescence times, path form (Hudson 1983; Kaplan & Hudson 1985): for the two-locus ancestral process started from two linked sequences, E (T_A T_B) - E T_A E T_B = (ρ + 18) / (ρ² + 13ρ + 18), with all three expectations finite lower Lebesgue integrals of the hitting times.

theorem covariance_coalescenceTimes {ρ : ℝ} (hρ : 0 ≤ ρ) :
    (∫⁻ ω, (⨅ (t : ℝ≥0) (_ : coalescedA (ctmcPath (twoLocusRate ρ) linked ω t) = true),
        (t : ℝ≥0∞)) *
      (⨅ (t : ℝ≥0) (_ : coalescedB (ctmcPath (twoLocusRate ρ) linked ω t) = true),
        (t : ℝ≥0∞)) ∂ctmcInputMeasure (isRateMatrix_twoLocusRate hρ)).toReal -
      (∫⁻ ω, (⨅ (t : ℝ≥0) (_ : coalescedA (ctmcPath (twoLocusRate ρ) linked ω t) = true),
        (t : ℝ≥0∞)) ∂ctmcInputMeasure (isRateMatrix_twoLocusRate hρ)).toReal *
      (∫⁻ ω, (⨅ (t : ℝ≥0) (_ : coalescedB (ctmcPath (twoLocusRate ρ) linked ω t) = true),
        (t : ℝ≥0∞)) ∂ctmcInputMeasure (isRateMatrix_twoLocusRate hρ)).toReal =
      twoLocusCov ρ
Uses 12
41 fully proved

The dependency explorer needs JavaScript and its graph data.

Download the complete graph data.