Zuker's minimum free energy algorithm

StructureNearest-neighbor loop energy model#✓fully proved

A loop (nearest-neighbor) energy model: energies of hairpin loops closed by (i, j), of stacks, bulges and interior loops closed by (i, j) with inner pair (k, l), and the affine multiloop parameters a (per multiloop), b (per branch) and c (per unpaired base). The exterior loop has energy 0. Units: free energy, e.g. kcal/mol.

structure LoopEnergy where
  /-- Energy of the hairpin loop closed by the pair `(i, j)`. -/
  hairpin : ℕ → ℕ → ℝ
  /-- Energy of the stack, bulge or interior loop closed by `(i, j)` with inner pair `(k, l)`. -/
  interior : ℕ → ℕ → ℕ → ℕ → ℝ
  /-- Multiloop initiation energy `a`. -/
  multiA : ℝ
  /-- Multiloop energy per branch `b`. -/
  multiB : ℝ
  /-- Multiloop energy per unpaired base `c`. -/
  multiC : ℝ
Used by 10

DefinitionStrict enclosure of one base pair by another#✓fully proved

q lies strictly inside the pair p: p.1 < q.1 and q.2 < p.2.

def Encloses (p q : ℕ × ℕ) : Prop :=
  p.1 < q.1 ∧ q.2 < p.2
Used by 5

DefinitionBranches of a loop closed by a base pair#✓fully proved

The branches of the loop closed by p in P: the pairs of P strictly inside p that are not strictly inside another pair of P strictly inside p.

def branches (P : Finset (ℕ × ℕ)) (p : ℕ × ℕ) : Finset (ℕ × ℕ) :=
  P.filter fun q => Encloses p q ∧ ∀ r ∈ P, Encloses p r → ¬ Encloses r q
Uses 1 · Used by 7

DefinitionExterior branches of a secondary structure#✓fully proved

The exterior branches of P: the pairs not strictly inside any pair of P.

def exteriorBranches (P : Finset (ℕ × ℕ)) : Finset (ℕ × ℕ) :=
  P.filter fun q => ∀ r ∈ P, ¬ Encloses r q
Uses 1 · Used by 3

TheoremUnique loop decomposition#✓fully proved

In a secondary structure every pair is either an exterior branch or a branch of exactly one pair of the structure (and not both): each pair is shared by the loop it closes and exactly one enclosing loop, possibly the exterior loop.

theorem existsUnique_closing_pair (P : Finset (ℕ × ℕ)) (hP : IsStructure pair m s P)
    (q : ℕ × ℕ) (hq : q ∈ P) :
    (q ∈ exteriorBranches P ∧ ∀ p ∈ P, q ∉ branches P p) ∨
      (q ∉ exteriorBranches P ∧ ∃! p, p ∈ P ∧ q ∈ branches P p)
Uses 4

DefinitionEnergy of the loop closed by a base pair#✓fully proved

The energy of the loop closed by p = (i, j) in P: a hairpin with no branch, a stack, bulge or interior loop with one branch (k, l), and the affine multiloop energy a + b · r + c · u with r ≥ 2 branches and u = (j - i - 1) - ∑ (l - k + 1) unpaired bases.

noncomputable def loopEnergy (E : LoopEnergy) (P : Finset (ℕ × ℕ)) (p : ℕ × ℕ) : ℝ :=
  let br := branches P p
  if br.card = 0 then E.hairpin p.1 p.2
  else if h : br.card = 1 then
    let q := (Finset.card_eq_one.mp h).choose
    E.interior p.1 p.2 q.1 q.2
  else
    E.multiA + E.multiB * br.card +
      E.multiC * ((p.2 - p.1 - 1 - ∑ q ∈ br, (q.2 - q.1 + 1) : ℕ) : ℝ)
Uses 2 · Used by 5

DefinitionFree energy of a structure in the loop model#✓fully proved

The free energy of a structure in the loop model: the sum of the energies of the loops closed by its pairs (the exterior loop contributes 0).

noncomputable def freeEnergy (E : LoopEnergy) (P : Finset (ℕ × ℕ)) : ℝ :=
  ∑ p ∈ P, loopEnergy E P p
Uses 2 · Used by 2

TheoremA pair with no branch closes a hairpin loop#✓fully proved

A pair with no branch closes a hairpin loop.

theorem loopEnergy_eq_hairpin (E : LoopEnergy) {P : Finset (ℕ × ℕ)} {p : ℕ × ℕ}
    (h : branches P p = ∅) : loopEnergy E P p = E.hairpin p.1 p.2
Uses 3 · Used by 1

TheoremA pair with one branch closes a stack, bulge or interior loop#✓fully proved

A pair with exactly one branch q closes a stack, bulge or interior loop with inner pair q.

theorem loopEnergy_eq_interior (E : LoopEnergy) {P : Finset (ℕ × ℕ)} {p q : ℕ × ℕ}
    (h : branches P p = {q}) : loopEnergy E P p = E.interior p.1 p.2 q.1 q.2
Uses 3 · Used by 1

TheoremA pair with two or more branches closes a multiloop#✓fully proved

A pair with at least two branches closes a multiloop with the affine energy.

theorem loopEnergy_eq_multiloop (E : LoopEnergy) {P : Finset (ℕ × ℕ)} {p : ℕ × ℕ}
    (h : 2 ≤ (branches P p).card) : loopEnergy E P p = E.multiA + E.multiB * (branches P p).card +
      E.multiC * ((p.2 - p.1 - 1 - ∑ q ∈ branches P p, (q.2 - q.1 + 1) : ℕ) : ℝ)
Uses 3 · Used by 1

TheoremEvery position has a unique innermost strict container#✓fully proved

Every position has a unique innermost strict container (not loop membership): a position x either lies strictly inside no pair, or there is a unique pair p with p.1 < x < p.2 such that no branch of p has x strictly inside it. This is an innermost-container assignment, not the conventional loop decomposition of positions: endpoints are excluded from the loop they close (a position x = q.1 or x = q.2 is assigned to the loop of which q is a branch, never to the loop closed by q), and no x < s.length hypothesis is assumed, so positions outside the sequence fall in the first (exterior) case. In particular do not read this as saying loops uniquely decompose positions; see existsUnique_closing_pair for the closing-pair side of the convention.

theorem existsUnique_loop_of_pos (P : Finset (ℕ × ℕ)) (hP : IsStructure pair m s P) (x : ℕ) :
    (∀ p ∈ P, ¬ (p.1 < x ∧ x < p.2)) ∨
      ∃! p, p ∈ P ∧ (p.1 < x ∧ x < p.2) ∧ ∀ q ∈ branches P p, ¬ (q.1 < x ∧ x < q.2)
Uses 4

DefinitionZuker's recursion#✓fully proved

Zuker's recursion on the interval [i, i + len), with j = i + len - 1, returning (V, WM): * V = ⊤ unless (i, j) can close a loop (i + m < j, CanPair pair s i j); otherwise V is the least of the hairpin energy hairpin i j, the interior-loop energies interior i j k l + V(k, l) over i < k < l < j, and the multiloop energies a + WM(i + 1, u) + WM(u + 1, j - 1) over i < u < j - 1; * WM (a multiloop fragment with at least one branch, charging b per branch and c per unpaired base) is the least of V + b, WM(i + 1, j) + c (base i unpaired), WM(i, j - 1) + c (base j unpaired) and WM(i, u) + WM(u + 1, j) over i ≤ u < j. Intervals are written by start and length.

Cslib.Bioinformatics.RNA.Zuker.zukerVW · 30 lines
noncomputable def zukerVW (E : LoopEnergy) : ℕ → ℕ → WithTop ℝ × WithTop ℝ
  | i, len =>
    let j := i + len - 1
    let V : WithTop ℝ :=
      if i + m < j ∧ CanPair pair s i j then
        min (E.hairpin i j : WithTop ℝ)
          (min
            (Finset.univ.inf fun a : Fin len => Finset.univ.inf fun b : Fin len =>
              if 1 ≤ a.val ∧ a.val + b.val < len ∧ 1 ≤ b.val then
                (E.interior i j (i + a.val) (i + a.val + b.val - 1) : WithTop ℝ) +
                  (zukerVW E (i + a.val) b.val).1
              else ⊤)
            (Finset.univ.inf fun a : Fin len =>
              if 1 ≤ a.val ∧ a.val + 2 < len then
                (E.multiA : WithTop ℝ) + (zukerVW E (i + 1) a.val).2 +
                  (zukerVW E (i + a.val + 1) (len - a.val - 2)).2
              else ⊤))
      else ⊤
    let WM : WithTop ℝ :=
      if len = 0 then ⊤
      else
        min (V + (E.multiB : WithTop ℝ))
          (min ((zukerVW E (i + 1) (len - 1)).2 + (E.multiC : WithTop ℝ))
            (min ((zukerVW E i (len - 1)).2 + (E.multiC : WithTop ℝ))
              (Finset.univ.inf fun a : Fin len =>
                if 1 ≤ a.val then (zukerVW E i a.val).2 + (zukerVW E (i + a.val) (len - a.val)).2
                else ⊤)))
    (V, WM)
  termination_by _ len => len
  decreasing_by all_goals omega
Uses 2 · Used by 2

DefinitionZuker exterior-loop recursion#✓fully proved

The exterior-loop recursion on the prefix [0, len): W(0) = 0 and W(len + 1) = min (W(len), min_{i ≤ len} W(i) + V(i, len + 1 - i)) (the last base is unpaired, or closes an exterior branch (i, len)).

noncomputable def zukerExterior (E : LoopEnergy) : ℕ → WithTop ℝ
  | 0 => 0
  | len + 1 =>
    min (zukerExterior E len)
      (Finset.univ.inf fun i : Fin (len + 1) =>
        zukerExterior E i.val + (zukerVW pair m s E i.val (len + 1 - i.val)).1)
  termination_by len => len
  decreasing_by all_goals omega
Uses 2 · Used by 2

DefinitionZuker minimum free energy#✓fully proved

Zuker's minimum free energy of s: the exterior recursion on the whole sequence.

noncomputable def zukerMFE (E : LoopEnergy) : WithTop ℝ :=
  zukerExterior pair m s E s.length
Uses 2 · Used by 2

TheoremZuker's algorithm#✓fully proved

Zuker's algorithm (Zuker and Stiegler 1981; Zuker and Sankoff 1984), here with the affine multiloop energy of Jaeger, Turner and Zuker (1989): the Zuker recursion computes the minimum free energy over all secondary structures on s with minimum hairpin m, for any real loop energies (of either sign). The minimum ranges over a finite set that contains the empty structure, so it is never ⊤.

theorem zukerMFE_eq_inf (E : LoopEnergy) :
    zukerMFE pair m s E =
      (structures pair m s).inf fun P => ((freeEnergy E P : ℝ) : WithTop ℝ)
Uses 16 · Used by 1

TheoremMinimum free energy is at most zero#✓fully proved

The minimum free energy is at most 0, the energy of the open chain (the empty structure).

theorem zukerMFE_le_zero (E : LoopEnergy) : zukerMFE pair m s E ≤ 0
Uses 7
16 fully proved

The dependency explorer needs JavaScript and its graph data.

Download the complete graph data.