Markov Chain-Based Sampling for Exploring RNA Secondary Structure under the Nearest Neighbor Thermodynamic Model and Extended Applications.

Ribonucleic acid (RNA) secondary structures and branching properties are important for determining functional ramifications in biology. While energy minimization of the Nearest Neighbor Thermodynamic Model (NNTM) is commonly used to identify such properties (number of hairpins, maximum ladder distance, etc.), it is difficult to know whether the resultant values fall within expected dispersion thresholds for a given energy function. The goal of this study was to construct a Markov chain capable of examining the dispersion of RNA secondary structures and branching properties obtained from NNTM energy function minimization independent of a specific nucleotide sequence. Plane trees are studied as a model for RNA secondary structure, with energy assigned to each tree based on the NNTM, and a corresponding Gibbs distribution is defined on the trees. Through a bijection between plane trees and 2-Motzkin paths, a Markov chain converging to the Gibbs distribution is constructed, and fast mixing time is established by estimating the spectral gap of the chain. The spectral gap estimate is obtained through a series of decompositions of the chain and also by building on known mixing time results for other chains on Dyck paths. The resulting algorithm can be used as a tool for exploring the branching structure of RNA, especially for long sequences, and to examine branching structure dependence on energy model parameters. Full exposition is provided for the mathematical techniques used with the expectation that these techniques will prove useful in bioinformatics, computational biology, and additional extended applications.

Introduction

Computational and mathematical applications play a critical role in the analysis of the structure and function of biological molecules, including ribonucleic acid (RNA). RNA is an essential biological polymer with many roles including information transfer, regulation of gene expression, and catalysis of chemical reactions. The primary structure of an RNA molecule may be understood as a sequence of amino acids: arginine, urasil, guanine, and cytosine. As is standard, we frequently abbreviate these as A, U, G, and C, respectively. RNA molecules are single-stranded and may therefore interact with themselves, forming A–U, G–U, and G–C bonds. The secondary structure of an RNA molecule is a set of such bonds.

The determination of secondary structure is an important step to understanding an RNA molecule’s full shape and therefore its function [1,2]. Accordingly, secondary structure information is commonly used in tertiary structure prediction algorithms, see, e.g., [3-6]. Identifying the secondary structure of RNA is crucial to understanding its function and mechanism in a cell [7]. Thus, the structure of RNA is critical to the development of biological and pharmaceutical therapeutics. Biologists use inexpensive and expedient means to sequence RNA, but the experimental determination of structure is more difficult and time-consuming. Therefore, computational methods are the primary means to determine possible RNA secondary structures.

For decades, one of the main computational approaches for examining RNA structure and branching properties has been thermodynamic free energy minimization using Nearest Neighbor Thermodynamics Modeling (NNTM) [8-10]. This free energy is in turn used in algorithms to predict secondary structure given an RNA sequence, see, e.g., [11-13]. Under the NNTM, the free energy of a structure is computed as the sum of the free energy of its various substructures. Many common programs (e.g., mFold, RNAFold, RNA Structure, sFold, Vienna RNA, etc.) intake a single sequence to produce secondary structures based on NNTM energy minimizations performed via dynamic programming. Nearest neighbor parameter sets include both a set of rules, referred to as equations or features, and a set of parameter values used by the equations. Separate rules exist for predicting stabilities of helices, hairpin loops, small internal loops, large internal loops, bulge loops, multi-branch loops, and exterior loops. Other branching properties of interest include, but are not limited to, average ladder distance, maximum ladder distance, maximum branching degree, average contact distance, average branching degree, degree of branching at the exterior loop, number of multi-loops with n braches, etc. The online nearest neighbor database (NNDB) archives and stores complete nearest neighbor sets, including rules and corresponding parameter values [14].

A common challenge is inferring whether the predicted results of NNTM for a set of RNA structural features or branching properties are within expected dispersion thresholds for a given energy model. For example, is the number of hairpins more than 2–3 standard deviations greater than the expected mean for a given energy model? This challenge is particularly vexing if the sequence is relatively long (greater than 1000 nucleotides). If structural features or branching properties are determined to exceed expected energy model dispersion thresholds, it relays potential scientific and/or mechanistic insight. Continuing with our hairpin example, what if an NNTM model produces a result where the number of hairpins seems rather large for the given sequence length? If the number of hairpins exceeds the expected dispersion of the NNTM model, it might be inferred that the greater number of hairpins is evidence of natural selection.

The primary objective of the present study is to enable mathematical determination of the dispersion of RNA secondary structural features for a given sequence length. We present a Markov-based algorithm to provide samples of the branching structure under the NNTM and Gibbs distribution, but without reference to a particular sequence of nucleotides. The algorithm enables the determination of where the predicted feature or branching property for an actual sequence falls within this distribution, which in turn enables the determination of whether the predicted NNTM feature or branching property is within expected dispersion limits.

In particular, this work investigates RNA substructures called multi-loops, the places where three or more helices join. Though multi-loops are crucial to the overall shape of a secondary structure, the models used to predict them algorithmically do not produce accurate results [15]. This investigation builds on an existing model of RNA branching [16] and provides a theoretical grounding for a Markov chain which may be used to algorithmically investigate branching properties of secondary structure models. The investigational foundation is a model for RNA secondary structure developed by Hower and Heitsch [16], in which secondary structures are in bijection with plane trees and the minimal energy structures of the model have been previously characterized. The present study characterizes the full Gibbs distribution of possible structures. Notably, Bakhtin and Heitsch [17] analyzed a very similar model and determined degree sequence properties of the distribution of plane trees asymptotically. However, the present study utilizes a Markov chain-based sampling algorithm to investigate the Gibbs distribution in the finite case. A full explanation of the plane tree model as well as the derivation of the energy functions is provided in Section 2.1.

Methods

The methods are divided into an overview of the RNA secondary structure NNTM plane tree model and energy functions (Section 2.1) and an all-encompassing explanation of the mathematical preliminaries that lay the foundation for the derived results and corresponding algorithms (Section 2.2).

Derivation of Energy Functions

The energy function studied here is derived from the Nearest Neighbor Thermodynamic Model (NNTM). The numerical parameters from the NNTM can be found in the NNDB [14]. In calculating energy functions for the sequences, we consider thermodynamic parameter values published by Turner in 1989 [8], 1999 [9], and 2004 [10].

The plane trees that we study in this paper come from two combinatorial RNA sequences, both of the form A4(Y5ZA4YZ5A4). The sequences of interest have (Y, Z) = (C, G) or (Y, Z) = (G, C). For both of these sequences, the set of maximally-paired secondary structures is in bijection with the set of plane trees of size n [18]. Figure 1 shows one example of a secondary structure and corresponding plane tree.

Figure 1.

A ribonucleic acid (RNA) secondary structure for one of the combinatorial RNA sequences used in this work and its corresponding plane tree. The ordering of the edges in the plane tree is derived from the 3’ to 5’ ordering of the RNA sequence. Note that the exterior loop corresponds to the root of the plane tree. The diagram in (a) was generated by ViennaRNA [19]. (a) A maximally-paired secondary structure for A4(C5GA4CG5A4)4 has 4 helices; (b) The corresponding plane tree has 4 edges and encodes the branching pattern seen in the secondary structure.

These specific combinatorial sequences are chosen because they allow for the study of the relationship between NNTM multiloop parameters and the branching behavior of secondary structures without interference from the energy contributions have specific base pairing combinations. In particular, the only places where the free energy differs between different secondary structures (for the same sequence) is in the type and number of multi-loops, the branching at the exterior loop, the number of hairpins, and the number of internal nodes. All of these energies directly relate to branching, not to specific base pairs. This simplification achieved by focusing only on multi-loops and branching both creates a model that is more amenable to theoretical analysis and speed computation.

Note that these secondary structures should not be considered representative of naturally occurring secondary structures. Instead, the only properties of interest in these structures are branching-related properties.

Three constants determine the free energy contribution of multiloops under NNTM, a, b, and c. The value of a encodes the energy penalty per multiloop. The constant b specifies the energy penalty per single-stranded nucleotide in a multiloop. The value of c gives the energy penalty for each helix branching from a multiloop.

In addition to the multiloop parameters a, b, c discussed above, we must account for the energy contributions of stacking base pairs, hairpins, interior loops, and dangling energy contributions. The energy of one helix is given by h. The energy associated with a hairpin is f, and the energy contribution of an interior loop is i. Finally, the parameter g encodes the dangling energy contributions. All of these values can be computed directly from the parameters found in the NNTM.

We wish to compute the energy of the structure corresponding to plane tree t having (down) degree sequence d0, d1, …, d and root degree r. Note that the down degree of a node x is equal to the number of children of x, and, in the down degree sequence, d is the number of non-root nodes with exactly i children. The energy contribution of all hairpin loops will be d0 f, and similarly, the total energy of all interior loops will be d1i. For a multi-loop having down degree j, the energy contribution will be a + 4b(j + 1) + c(j + 1) + (j + 1)g, and so the contribution of all multi-loops is given by . Te root vertex of the tree corresponds to the exterior loop and has energy contribution gr. Finally, our structure has n helices, each with energy h. Summing all of these components gives the total energy. where we have used the facts and .

Set α = f − a − 4b − c − g, β = i − a − 8b − 2c − 2g, γ = −4b − c, and δ = a + 8b + 2c + h + 2g. Then, the energy function is αd0 + βd1 +γr +δn. Since n will be fixed, we disregard the term δn, giving

Though we study these energy functions for arbitrary values of (α, β, γ), numerical values for both the input energy parameters from NNTM and the resulting energy function coefficients are given in Table 1.

Table 1.

Nearest Neighbor Thermodynamic Model (NNTM) parameters and resulting energy functions. Energy functions are of the form αd + βd1 + γr.

Y	Z	Turner	a	b	c	h	f	i	g	α	β	γ
C	G	89	4.6	0.4	0.1	−10.9	3.8	3.0	−1.6	−0.9	−1.8	−1.7
G	C	89	4.6	0.4	0.1	−16.5	3.5	3.0	−1.9	−0.9	−1.2	−1.7
C	G	99	3.4	0	0.4	−12.9	4.5	2.3	−1.6	2.3	1.3	−0.4
G	C	99	3.4	0	0.4	−16.9	4.1	2.3	−1.9	2.2	1.9	−0.4
C	G	04	9.3	0	−0.9	−12.9	4.5	2.3	−1.1	−2.8	−3.0	0.9
G	C	04	9.3	0	−0.9	−16.9	4.1	2.3	−1.5	−2.8	−2.2	0.9

Mathematical Preliminaries

Section 2.2 of this manuscript provides the necessary mathematical background, including a formal introduction of combinatorial objects and a review of the relevant Markov chain mixing results used to construct our resultant sampling Markov chain and corresponding mixing time proof in Section 3.

Combinatorial Objects

A plane tree is a rooted, ordered tree. We will use to denote the set of plane trees with n edges. It is known that is given by the nth Catalan number . In a plane tree, a leaf is a node with down degree 0, and an internal node is a non-root node with down degree 1. For a given plane tree t, we will use d0(t) to denote the number of leaves and d1(t) to denote the number of internal nodes.

For a plane tree t, the energy of the tree is given by where α and β are real parameters of the energy function. Note that this function is a simplification of the model due to Hower and Heitsch [16] discussed in Section 2.1. Making this simplification effectively disregards the energy contribution of the exterior loop, which is small in comparison to the total energy of a structure, especially for the longer sequences that are of interest to us. Other authors have made similar simplifications, e.g., [17].

For our purposes, we consider α and β to be arbitrary but fixed. We will consider a Gibbs distribution g on the set , where the weight of each tree t is given by where is a normalizing constant.

A Motzkin path of length n is a lattice path from (0, 0) to (n, 0), which consists of steps along the vectors U = (1, 1), H = (1, 0), and D = (1, −1) and never crosses below the x-axis. We can also represent Motzkin paths as strings from the alphabet {U, H, D} where, in any prefix, the number of Us is greater than or equal to the number of Ds. The number of Motzkin paths of length n is given by the Motzkin numbers M where

Motzkin numbers and Motzkin paths have been well-studied in the combinatorics literature, see, e.g., [20-24].

A Dyck path is a Motzkin path with no H steps. It is easy to see that a Dyck path must have even length, so we will use to denote the set of Dyck paths on length 2n. It is well known that (see, e.g., [25]).

A 2-Motzkin path is a Motzkin path in which (1, 0) steps are given one of two distinguishable colors. Let be the set of all 2-Motzkin paths of length m. We can also represent 2-Motzkin paths as strings from the alphabet {U, H, I, D}, where, as before, the number of Ds never exceeds the number of Us in any prefix. In a such a string x, we denote by |x| the number of times the symbol a appears in x, where a ∈ {U, H, I, D}. Notice that we always have |x| = |x|. For any and k ∈ {1, · · ·, n}, let x(k) denote the symbol at index k in the string representation of x. Additionally, the skeleton of a 2-Motzkin path x is the Dyck path of Us and Ds which results from removing all Hs and Is from x. We will denote the skeleton of x by σ(x).

A Bijection Between and

We will use the particular bijection between plane trees and 2-Motzkin paths from Deutsch [26], which neatly encodes information about d0 and d1. For clarity, we will overview the bijection here.

For a given plane tree t with n edges, assign a label from the set {U, H, I, D} to each edge e according to the following rules:

If e is the leftmost edge off a non-root node of down degree at least 2, assign the label U.

If e is the rightmost edge off a non-root node of down degree at least 2, assign the label D.

If e is the only edge off a non-root node of degree 1, assign the label I.

If e is an edge off the root node, or if e is neither the leftmost nor the rightmost edge off its parent node, assign the label H.

Now, if we traverse t in a preorder reading off these labels, we get a 2-Motzkin path of length n. However, this path will always begin with H, so we define Φ(t) to be the 2-Motzkin path of length n − 1 after this initial H is removed. Figure 2 gives an example of this labeling process. From Deutsch, we know not only that Φ is a bijection, but also that if x = Φ(t) then |x| = d1(t) and |x| + |x| + 1 = d0(t).

Figure 2.

A plane tree with edges labeled according to the bijection Φ, along with its corresponding 2-Motzkin path.

Using this bijection, it is natural to extend our energy function to 2-Motzkin paths. We define the energy of a 2-Motzkin path x to be and we extend our definition of the distribution g to accordingly. We note that, while this energy function does not capture all possible weightings on 2-Motzkin paths, it does capture all weightings possible under our simplification of the model due to Hower and Heitsch [16] after applying the bijection due to Deutsch [26].

Markov Chains

A Markov chain is a sequence of random variables X0, X1, X2, · · · taking values in a state space Ω subject to the condition that

All Markov chains that we consider will be implicitly time-homogeneous (meaning Pr(X = y | X = x) does not depend on t) and finite (meaning |Ω| < ∞). The transition matrix of a time-homogeneous Markov chain is the matrix P: Ω × Ω → [0, 1] given by

It is easy to see that if X0 has distribution vector x, then X has distribution vector Px.

A finite Markov chain with transition matrix P is said to be ergodic if it has the following two properties.

Irreducibility: For any x, y ∈ Ω, there is some integer for which P(x, y) > 0.

Aperiodicity: For any state x ∈ Ω, we have .

It is well-known that if is ergodic, then there exists a unique distribution vector π, the stationary distribution, such that Pπ = π, and lim P(x, y) = π(y) for any states x, y ∈ Ω. Additionally, we call reversible if for all states x, y ∈ Ω, we have π(x)P(x, y) = π(y)P(y, x).

For ϵ > 0, the mixing time τ(ϵ) of is given by

Intuitively, the mixing time gives a measure of the number of steps required for to get sufficiently close to its stationary distribution from any starting state.

Let be a finite ergodic Markov chain over a state space Ω with transition matrix P. Let the eigenvalues of P be λ0, λ1, …, λ|Ω|−1 such that 1 = λ0 > |λ1| ≥ … ≥ |λ|Ω|−1|. The spectral gap of is given by . As is standard, it will be convenient to denote the inverse of the spectral gap by relaxation time .

Additionally, the spectral gap is given by the following functional definition [27]. where the infimum is taken over all non-constant functions . A direct consequence of this definition of the spectral gap is the following lemma.

Lemma 1.

Let and be ergodic Markov chains over Ω with the same stationary distribution. Let P1 and P2 be the transition matrices of and respectively. If for all x, y ∈ Ω and for some constant c > 0 we have P1(x, y) ≤ cP2(x, y), then .

Additionally, spectral gap is related to the mixing time by the following lemma [28].

Lemma 2.

Let be an ergodic Markov chain with state space Ω, and let λ0, λ1, …, λ|Ω|−1 be the eigenvalues of the transition matrix P as defined above. Then, for all ϵ > 0 and x ∈ Ω, we have

We say that a Markov chain , whose state space depends on a variable , is rapidly mixing if τ(ϵ) is bounded above by some polynomial in n and log(ϵ−1). For the specific chains studied in this manuscript, we will show that is bounded by a polynomial in n and log(ϵ−1) if and only if is bounded by a polynomial in n and log(ϵ−1). Our next lemma presents sufficient conditions.

Lemma 3.

Let be an ergodic Markov chain with state space Ω and let λ0, λ1, …, λ|Ω|−1 be the eigenvalues of its transition matrix. Let ϵ > 0. If τ(ϵ) is bounded by a polynomial in n and log(ϵ−1), then τ log(ϵ−1). Further, suppose we have log(1/π(x)) bounded by some polynomial q(n) for all x ∈ Ω. Then, being bounded by a polynomial in n and log(ϵ−1) implies that τ(ϵ) is also bounded by some polynomial in n and log(ϵ−1).

Proof.

Suppose that τ(ϵ) ≤ p(n, log(ϵ−1)), where p is a polynomial. Beginning with the left hand side of Lemma 2, note that

Then, applying Lemma 2 and the bound on τ(ϵ), where p′ is again a polynomial in n and log(ϵ−1).

Turning now to converse, suppose that we have τ ≤ p(n, log(ϵ−1)), for some polynomial p. Additionally suppose log(1/π(x)) ≤ q(n) for all x ∈ Ω, for some polynomial q.

Applying Lemma 2, where p′ is some polynomial. □

Coupling

A coupling of a Markov chain on Ω is a chain on Ω × Ω for which the following properties hold.

Each chain and , when viewed in isolation, is a copy of , given initial states X0 = x and Y0 = y.

Whenever X = Y, we have X = Y.

Formally, the first item above requires that the joint distribution of (X, Y) (given (X, Y)) should satisfy the property that the marginal of X (and also Y) is consistent with the probability transitions of . We define the coupling time T to be

The following lemma [29] is useful in bounding the coupling time T.

Lemma 4.

Suppose that is a coupling of a Markov chain M. Let φ be an integer-valued distance function on Ω × Ω taking values in the range [0, B], and suppose that φ(x, y) = 0 if and only if x = y. Let φ(t) = φ(x, y). Suppose that the coupling satisfies E (φ(t + 1) − φ(t)|X, Y) ≤ 0. Additionally, suppose that whenever φ(t) > 0, E |φ(t + 1) − φ(t)|2|X, Y ≥ V. Then, the expected coupling time satisfies E (Tx,y) ≤ φ(0)(2B − φ(0))/V.

Coupling time and mixing time are then related by the following theorem [28].

Theorem 1.

A Markov chain M with coupling time T has mixing time τ(ϵ) bounded by

Decomposition

We use two disjoint decomposition methods for bounding the spectral gap, one developed by Martin and Randall [30], and a very recent one given by Hermon and Salez [31], building on the work by Jerrum, Son, Tetali and Vigoda [32]. We use both theorems because, while the latter gives better bounds, the former has more relaxed conditions, which is necessary in one of our applications. The setup for both methods is the same.

Let be an ergodic, reversible Markov chain over a state space Ω with transition matrix P and stationary distribution π. Suppose Ω can be partitioned into disjoint subsets Ω1, …, Ω. For each i ∈ [m], let be the restriction of to Ω, which is obtained by rejecting any transition that would leave Ω. Let P be the transition matrix of Additionally, we define to be the projection chain of over the state space [m] as follows. Let the transition matrix of be given by

One can check that is reversible and has stationary distribution while each has stationary distribution

With this notation, we have the following theorem by Martin and Randall [30].

Theorem 2.

Defining and as above, we have

The theorem due to Hermon and Salez obtains better bounds if, for each pair (i, j) ∈ [m] × [m] with , we can find an effective joint distribution (often referred to as a “coupling”) κ : Ω × Ω → [0, 1] of the distributions π and π. In other words, we must have The quality of the joint distribution κ is defined as where the minimum is taken over all (x, y, i, j) with x ∈ Ω, y ∈ Ω for which and κ(x, y) > 0. Hermon and Salez [31] prove the following.

Theorem 3.

With P, , Pi, and χ defined as above,

The utility of these decomposition theorems is that they allow us to break down a more complicated Markov chain into pieces that are easier to analyze. If we can show that the pieces rapidly mix, and the projection chain rapidly mixes, then we may conclude that the original chain rapidly mixes as well.

Additionally, to aid with the analysis of some projection chains, we will need another lemma from [30].

Let be the Markov chain on [m] with Metropolis transitions whenever , where Δ is the maximum degree of vertices in the transition graph of . Let ∂(Ω) = {y ∈ Ω : ∃x ∈ Ω with P(x, y) > 0}. Then we have the following

Lemma 5.

With as defined above, suppose there exist constants a > 0 and b > 0 with

P(x, y) ≥ a for all x, y such that P(x, y) > 0.

π(∂(Ω)) ≥ bπ(Ω) for all i, j with

Then .

In order to help analyze the mixing time of , we will also require the following two lemmas. Note that Lemma 6 is used only in the proof of Lemma 7.

Lemma 6.

Let be a log concave sequence, with a > 0 for all 1 ≤ i ≤ m. Then, for all 1 ≤ i ≤ j ≤ m.

In order to use induction, we will slightly reframe the statement. We will prove for all i + k ≤ n.

We now proceed by induction on k. The base case, k = 0, is trivial.

Now fix l > 0 and suppose that the induction hypothesis is true for k = l − 1, that is, By log concavity , or, equivalently, Therefore, where the first inequality follows from the induction hypothesis, and the second inequality follows from log concavity. □

Lemma 7.

Let π be a probability distribution on [m]. Let be a Markov chain on [m] with the transition probabilities and the appropriate self-loop probabilities P(i, i). If π(i) is log-concave in i, then has mixing time (and hence also relaxation time) O(m2).

We define a coupling (X, Y) on as follows. If X ≠ Y, then at time step t + 1, flip a fair coin.

If heads, set Y = Y. Let l be either 1 or −1, each with probability 1/2. If possible, let X = X + l with probability . Otherwise, let X = X.

If tails, set X = X, and update Y the same way as we did for X in the previous case.

Now, suppose that for some t we have X = i and Y = j for i ≠ j. WLOG, assume that i < j. Let φ(t) = φ(X, Y) = j − i, and let Δφ(t) = φ(t) − φ(t − 1). Note that we have two moves, with probabilities P(i, i − 1) and P(j, j + 1), which will increase the distance φ by 1 and similarly two moves, with probabilities P(i, i + 1) and P(j, j + 1), will decrease the distance by 1. Then we have

By the log-concavity of π(i) and Lemma 6, we have P(i, i + 1) ≥ P(j, j + 1) and P(i, i − 1) ≤ P(j, j − 1). Therefore, the expected change in φ(t) is always non-positive. We also have

We claim that . Suppose, for contradiction, that the expectation is less than . Then, for each of the minimum functions in the above expression, 1 must be the larger argument. Equivalently, π(i − 1) < π(i), π(i) > π(i + 1), π(j − 1) < π(j), and π(j) > π(j + 1).

Therefore, π(i) is not unimodal in i and is therefore also not log-concave in i, contradicting our hypothesis. Therefore we have , as desired. □

Results

Here we present the constructed Markov chain and corresponding algorithms devised for the sampling task and the proof of an upper bound on the relaxation time—that the chain mixes rapidly. Collectively, the results illustrate an analytical approach to calculate the dispersion of the secondary structure and corresponding branching properties of RNA based on the NNTM energy function minimization and without reference to a specific nucleotide sequence.

Our Markov Chain on

We define a Markov chain on to sample 2-Motzkin paths as a representation of plane trees. Here, we use m = n − 1 to denote the length of the 2-Motzkin paths corresponding to plane trees with n edges.

We define each step of as follows. First, pick a random element l uniformly from {1, 2, 3, 4}. Now choose y as follows.

If l = 1, pick a random pair of consecutive symbols in X, and call this pair s. If s is UD or HH, let s′ be either UD or HH with probabilities and respectively. Let y be the string X with s replaced by s′. Otherwise, let y = X.

If l = 2, pick i uniformly from {1, · · ·, m}. If X(i) is H or I, choose a symbol c to be either H or I with probabilities and respectively. Let y be the 2-Motzkin path given by changing the symbol in X(j) to c. Otherwise, we let y = X.

If l = 3, pick i and j each uniformly from {1, · · ·, m}. If each of X(i) and X(j) are either U or D, let y be the string X with the symbols at indices i and j swapped. Otherwise, let y = X.

If l = 4, pick a random pair of consecutive symbols in X, and call this pair s. If s is of the form ab or ba for some a ∈ {U, D} and b ∈ {H, I}, let s′ be the reverse of s, and let y be the string X with s replaced by s′. Otherwise, let y = X.

If y is a valid 2-Motzkin path, set X = y with probability . Otherwise, set X = X.

One can see that is irreducible by noting that every path can be transformed to the path consisting of all H’s. To make this transformation, first use the l = 4 rule to move all H’s and I’s to the end of the path. If there are any U’s in the path, we must now have at least one consecutive pair UD. Use the l = 1 rule to convert the UD to a HH. From here we can repeat, again moving all H’s to the end and replacing UD with HH, until only H’s and I’s remain. Finally, we can use the l = 2 rule to convert all I’s to H’s. Since all of these steps can also be taken in reverse, this gives a procedure to move between two arbitrary paths, demonstrating irreducibility. We can also conclude that is aperiodic, due to the existence of self-loops. Combined with irreducibility, this establishes that is ergodic.

We claim that is reversible with respect to the stationary distribution , where . This can be easily verified by considering the four move types listed above. For example, for the first move type given above (transforming UD to HH and vice versa), let x and y be the states of interest. Suppose that y has the consecutive symbols HH where x contains UD. Then,

One can verify that similar computations hold for the remaining three types of moves. Therefore, we conclude that the chain has stationary distribution .

The Markov chain can be implemented in pseudocode as in Algorithm 1. Here, the Ber(p) function returns true with probability p, and false otherwise. We also use the addition of strings to denote concatenation.

Algorithm 1:

The main Markov chain algorithm. This pseudocode calculates X given X0.

Require: X0 is a valid 2-Motzkin path of length m.

x ← X0

for s = 1 → t do

y ← x

l ← randInt(1, 4)

if l = 1 then

i ← randInt(1, m − 1)

if x[i: i + 1] = UD and Ber(e^−α2(1+e^−α)) then

y[i: i + 1] ← HH

else if x[i: i + 1] = HH and Ber(12(1+e^−α)) then

y[i: i + 1] ← UD

else if l = 2 then

i ← randInt(1, m)

if x[i] = I and Ber(e^−α2(e^−α+e^−β)) then

y[i] ← H

else if x(i) = H and Ber(e^−β2(e^−α+e^−β)) then

y[i] ← I

else if l = 3 then

i ← randInt(1, m)

j ← randInt(1, m)

if (x[i] ∈ {U, D} and x[j] ∈ {U, D}) and Ber(12) then

y[i] ← x[j]

y[j] ← x[i]

if y is not a valid 2-Motzkin path then

y ← x

else if l = 4 then

i ← randInt(1, m − 1)

if (x[i] ∈ {U, D} and x[j + 1] ∈ {H, I}) or (x[i] ∈ {H, I} and x[j + 1] ∈ {U, D}) and Ber(12) then

y[i: i + 1] ← x[j + 1] + x[j]

x ← y

return x

Additionally, in order to convert the 2-Motzkin path X into a plane tree, we use the algorithm in Algorithm 2, which assumes the existence of a Node object with children and parent attributes.

Algorithm 2:

Algorithm to convert a sampled 2-Motzkin path to a plan tree. The pseudocode calculates φ−1(x).

Require: x is a valid 2-Motzkin path of length m.

root ← new Node()

// u will be where a new node will be added for an H or D symbol

u ← root

// v will be always the last node added

v ← new Node()

// the stack will keep track of previous values of u

stack = new Stack()

root.children.append(v)

for i = 1 → m do

node ← new Node()

if x[i] = U then

v.children.append(node)

stack.push(u)

u ← v

else if x[i] = I then

v.children.append(node)

else if x[i] = H then

u.children.append(node)

else if x[i] = D then

u.children.append(node)

u ← stack.pop()

v ← node

return root

Mixing Time Results

Our main result is to prove the rapid mixing of the Markov chain defined in Section 3.1. An upper bound on the relaxation time is achieved by bounding the spectral gap from below. A spectral gap bound for the complex chain at hand is obtained through the use of multiple decomposition theorems, which give bounds on the spectral gap of the complex chain in terms of the spectral gaps of multiple simpler chains. The disjoint decomposition theorem due to Martin and Randall [30] provides a flexible approach to the decomposition of Markov chains. Very recent work by Hermon and Salez [31], building on the work of Jerrum, Son, Tetali, and Vigoda [32], proves a decomposition theorem with tighter bounds but stronger hypotheses.

Since this proof involves multiple decomposition steps, we provide an overview here. The primary tools used in this proof are the two decomposition theorems presented in Section 2.2.5. We first partition the state space of all 2-Motzkin paths by the number of Us in the path. The projection chain from this first decomposition is linear and is proved to be rapidly mixing using a result of Martin and Randall [30] (Lemma 8). Each of the restriction chains are decomposed again, this time by the pattern of H and I symbols. The projection chains for this second decomposition are shown to be rapidly mixing by coupling (Lemma 9). The restriction chains are decomposed a third time, this time according to the skeleton of U and D steps. The projection chains for this third decomposition are shown to be rapidly mixing by comparison to the classic mountain valley moves chain on Dyck paths (Lemma 10). This last set of restriction chains are found to be rapidly mixing by isomorphism to the chain consisting of adjacent transpositions on binary strings (Lemma 11). Finally, starting from the most restricted chains, we use the decomposition theorems to obtain a bound on the spectral gap of the original chain (Theorem 4).

We now proceed with a formal presentation. We will use a series of decompositions of . We will first decompose our state space into , where Let denote the Markov chain restricted to the set S, and let be the projection chain over this decomposition as outlined for Theorem 2.

Additionally, we will decompose each S into the sets {T : q ∈ (H + I)}, where (H + I) denotes the set of strings with length m − 2k from the alphabet {H, I}. We define T to be the set of 2-Motzkin paths x ∈ S such that the substring of H and I symbols in x is q. Let denote the chain restricted to T, and let be the projection chain of over this decomposition.

Finally, we decompose each T into the partition based on the skeletons of the 2-Motzkin paths. For each , we define

As before, we let be the Markov chain restricted to U, and let be the appropriate projection chain. For clarity, this four-level decomposition is summarized in Figure 3.

Figure 3.

The four level decomposition of (left), and the projection chains corresponding to each decomposition (right).

Lemma 8.

has relaxation time .

Proof.

The chain is a linear chain with states k in {0, …, ⌊m/2⌋}, and with stationary distribution where is defined as in Section 2.2.5. Notice that transitions in which move between the S sets are those which change a HH substring into a UD or DU substring, or vise versa. Thus, the transitions in only increase or decrease k by at most 1. We seek to apply Lemma 5. To choose a, notice that for x ∈ S and y ∈ S with P(x, y) > 0, we have

Note that the factor 1/4 comes from the choice l = 4, and the factor 1/(m − 1) comes from the fact that there are m − 1 adjacent pairs to pick from. Then,

Thus, we pick .

To pick b, we let for k ∈ {1, · · ·, ⌊m/2⌋}, and we let for k ∈ {0, · · ·, ⌊m/2⌋ − 1}.

Additionally, let A for k ∈ {1, · · ·, ⌊m/2⌋} be the subset of S consisting of the 2-Motzkin paths in which the first D symbol appears immediately after a U. Let B for k ∈ {0, · · ·, ⌊m/2⌋ − 1} be the subset of S consisting of the 2-Motzkin paths in which a pair of adjacent H symbols occur before all other H or I symbols. It is easy to see that A ⊂ ∂−(S) and B ⊂ ∂+(S). We have as there are C ways to arrange the U and D symbols and ways to insert m − 2k H or I symbols (treating H and I as being identical for now) without placing anything between the first D and the U immediately before it. The energy contribution of the U and D symbols is given by e−, and the energy contribution of the H and I symbols is (e− + e−). The required normalizing constant is Z. Similarly, we also get because there are C ways to arrange the U and D symbols and ways insert m − 2k − 1 H or I symbols (treating the initial pair of H’s as a single symbol gives us only m − 2k − 1 symbols to insert). The energy contribution of the U’s, D’s, and the initial two H’s is given by , and the energy contribution of the remaining H’s and I’s is (e− + e−). Finally, Z is again a normalizing constant.

Hence combining these two results, we have and

Thus, we may let .

Applying Lemma 5, we get that . Additionally, one can check that is log concave in i. Hence, using Lemma 7, we get , and in turn , as claimed. □

Lemma 9.

has mixing time , for all k.

Notice that appears as a chain with states q in the set Q = (H + I). Additionally, transitions in only occur between strings in Q that differ at only one index. The stationary distribution of is given by , where we have intentionally used the constant of proportionality to remove all dependence on k, which we consider, in this context, to be fixed.

Additionally, for q1, q2 ∈ Q which differ at exactly one index, we have the transition probability

We may show that rapidly mixes by a simple coupling argument. Let be our coupled Markov chain on Q × Q. We define one step in this coupled chain as follows.

With probability , set (X, Y) = (X, Y).

Otherwise, pick a random index j ∈ [m − 2k]. Let a ∈ {H, I} be a random symbol such that and . Now let X and Y be X and Y respectively, each with the jth symbol changed to a.

One can check that each of (X) and (Y) are indeed copies of . Additionally, notice that we will have X = Y after all m − 2k possible indices j have been updated. By the Coupon Collector Theorem, we have the coupling time of this chain to be . Thus, using Theorem 1, we have the mixing time (and the relaxation time) also O(m log m). □

Lemma 10.

has relaxation time , for all pairs (k, q).

Notice that all x ∈ T have equal energy, and that for all s. Thus, has a uniform stationary distribution. If we represent each set U by the Dyck path s, we can think of as a chain over . Since all the transitions in that move between the U sets are moves that exchange the positions of a U and a D, the transitions in are simply the moves on elements of which exchange a U with a D. We call these moves on the elements of , transposition moves.

For each that differ by a transposition move, the transition probabilities in our projection chain are given by

The last equality above relies on counting the number of terms in the sum. Notice that for each , there is a unique for which P(x, y) > 0. Therefore, the number of terms is simply . Compare this chain to the traditional mountain valley Markov chain on , which we will denote by . The transition probabilities of are given by for each pair (s1, s2) which differ by a mountain-valley move. It is known from Cohen [33] that . Thus, applying Lemma 1 to and , we see that . □

Lemma 11.

has relaxation time , for all valid triples (k, q, s).

Notice that transitions in consist only of moves which involve swapping an H or an I with an adjacent U or D. Additionally, all 2-Motzkin paths in U have equal energy, so for all x, y ∈ U such that P(x, y) > 0, we have .

To determine the mixing time of , consider an isomorphic chain. Let U′ be the set of all binary strings of length m with 2k zeros and m − 2k ones. Let be the Markov chain on U′ where each step does nothing with probability 7/8 and swaps a random pair of adjacent (potentially identical) digits with probability 1/8. From Wilson [34], we know that the spectral gap of is . □

Finally, we can combine our bounds on the spectral gaps of all of these chains to prove our main result.

Theorem 4.

The Markov chain has relaxation time , for all .

We use Lemmas 11 and 10 with Theorem 3 to obtain a bound on . We define a coupling for each pair with . For each such pair, notice that the set of pairs with P(x, y) > 0 is a perfect matching. Thus, we may set

To compute χ, we begin by observing π(x) = π(y) for all . Note also for all skeletons s of length 2k. Before computing χ, we start by finding .

We now proceed with the calculation of χ. Recall that the minimum is taken over all tuples x, y, s1, s2 where and .

Theorem 3 then gives

Similarly, we define a coupling for each pair (q1, q2) ∈ (H + I) × (H + I) with to apply Theorem 3 to . Notice that once again, the set of pairs for which P(x, y) > 0 forms a perfect matching. Thus, we take

To compute χ for this coupling, we again begin with a few preliminary computations. In all of the following, let , with . Note that q1 and q2 have the same length and differ at only one index. We will show the computations for the case where q1 has a I where q2 has a H. The computations for the other case are nearly identical.

Note that . Note also and

Now we can compute

Applying Theorem 3 then gives

Unfortunately, we have not been able to find a useful coupling for , so for the last step of our decomposition, we apply Theorem 2. Since and for all k, we have establishing Theorem 4. □

Finally, an application of Lemma 3 allows us to conclude that the mixing time is also polynomially-bounded.

Corollary 1.

is rapidly mixing.

In order to apply Lemma 3, we need to obtain a polynomial bound on log(1/π(x)) for all x ∈ Ω. Let t ∈ Ω have maximum energy among all elements of Ω. For any x ∈ Ω,

This gives us the required polynomial bound, and therefore Lemma 3 implies that is rapidly mixing. □

Discussion and Conclusions

The goal of this work was to identify a Markov chain and construct a corresponding algorithm by which to examine the non-uniform distribution and dispersion properties of NNTM RNA secondary structures and branching properties independent of a specific nucleotide sequence. This study successfully identifies the existence of a Markov chain, with a provably polynomial mixing time, which generates a Gibbs distribution on plane trees. This stationary probability distribution models branching characteristics of RNA secondary structure under the NNTM. While the exploration of sampled structures obtained from this algorithm are beyond the scope of the presented results, pseudocode (see Section 3.1) is provided to facilitate future work in this area. Below we discuss the direct applications and implications of this work to RNA modeling, the possibility of implementing a dynamic programming approach, the possibility of an approach using stochastic context-free grammars, other biological applications of this work, contributions of this work towards independent mathematical research interests, and limitations and future directions of the present work.

Applications to RNA Modeling

The most straightforward application of this work is in understanding the background distribution of the branching behavior for secondary structures predicted under the NNTM. While the NNTM is widely used to predict secondary structures from sequence data, little is known about the general branching characteristics of the predicted structures, independent of a specific input sequence. Quantities such as the number of hairpins, the maximum branching in a multiloop, the average branching in a multiloop, and the maximum ladder distance of the structure [7,35] help to characterize the branching behavior and could be computed from samples obtained from this algorithm. These quantities also have been studied in native structures and/or could be easily obtained from databases such as RNA STRAND [36]. The parameter values of α, β, and γ corresponding to various revisions of the NNTM are given in Table 1 in Section 2.1. The Markov chain and corresponding algorithms presented will enable biologists to calculate the dispersion of key branching properties for a specific energy function. As described with the detailed hairpin dispersion example in the Introduction (Section 1), knowing whether branching properties fall within acceptable dispersion limits is crucial for deducing potential functional insight or hypothesizing other scientific ramifications.

Another key application to RNA modeling of the presented algorithms is the ability to explore the parameter space of possible values for α and β. While the various revisions of the NNTM correspond to specific values for these parameters, in principle any real-valued parameters could be used. Finding values for these parameters that approximate reality remains an open question. Yet, determination of how differences in parameter values change the distribution of NNTM branching properties, such as maximum ladder distance, is crucial. Moreover, parameter space exploration is necessary to identify and further explore the phase transitions that exist. The presented Markov chain and corresponding algorithms expedite such future computational experimentation. Therefore, collectively, the presented algorithm enables exploration that will greatly improve understanding of NNTM-based RNA secondary structures and branching properties, as well as identify potential limitations or specific branching structures where the NNTM models do not sufficiently emulate reality. For example, NNTM-based free energy minimization algorithms achieved an accuracy of at least 60% in only 9% of 16S secondary structures analyzed by Doshi et. al. [15].

The algorithm presented here can only sample under an energy function of the form αd0 + βd1, and this does not capture the entirety of the model presented in [16], which considers energy functions of the form αd0 + βd1 + γr. However, the missing term, γr, represents the energy contribution of the exterior loop, and the exterior loop contributes less of the total free energy as sequence length increases. Therefore, when interested in sequences of at least moderate length, this algorithm may be able to provide insight, as long as information about the exterior loop is not the specific object of study. Note that other authors have made similar simplifications with respect to the exterior loop, e.g., [17].

Possibility of a Dynamic Programming Approach

This sampling problem to calculate the dispersion of NNTM RNA secondary structure and properties utilized Markov chain techniques. However, is it possible to utilize a dynamic programming algorithm? It is straightforward to sample Dyck paths under a uniform probability distribution using dynamic programming techniques. However, it is not clear whether a similar technique could be used for the Gibbs distribution we define here, due to the complexity of the energy function. In particular, large numeric computations may be required to handle the variable k, the number of U steps in a path. While Alonso presents a way to sample from the unweighted distribution time without large computations [37], it is unclear if a similar method may be used for the present application.

Possibility of an SCFG Approach

Stochastic context-free grammars (SCFGs) have been widely used in the field of RNA secondary structure prediction, e.g., [38-41]. Most commonly, the probabilities for production rules in an SCFG are determined by training on a set of known secondary structures, often including covariance information from homologous structures. These approaches are not immediately applicable to the problem we study here, as they do not give any insight into the NNTM multiloop energy parameters.

However, some authors have constructed SCFGs based on the NNTM. In particular, Nebel and Scheid [38] construct an SCFG with 29 distinct production rules to mirror the NNTM features. They also present a sampling algorithm allowing for sampling structures of a fixed size using the grammar. However, they do not actually compute probabilities for the production rules that would allow one to sample from a Gibbs distribution (with NNTM energy) and instead rely on training on a set of known structures. Indeed, it is not clear from the paper whether such a set of probabilities must exist.

Even in the case of the simplified model we present in this manuscript, it is not clear how to assign probabilities to production rules in an SCFG so that the probability of obtaining a given structure matches the Gibbs probability under the NNTM. See Supplement 1: SCFG for more details.

Even if a suitable SCFG could be formulated, the SCFG approach is not necessarily superior. The sampling algorithm presented by Nebel and Scheid has time complexity O(n3) and space complexity O(n2). While the algorithm we present does have large time complexity, it only requires linear space, which may be an advantage for some applications.

Even though we cannot easily formulate a SCFG, it is reasonable to consider whether a context-free grammar (such as that presented in Supplement 1) could nonetheless be used as the basis for a dynamic programming algorithm. In fact, this is possible. The key idea is to create a table for each non-terminal symbol X and then populate entry k of the table with where the sum is taken over all trees which can be derived from symbol X.

Once the tables been populated with these (non-normalized) probabilities, a stochastic backtracking procedure can be used to obtain samples.

However, as in Section 4.2, an assumption that each arithmetic operation can be performed in unit time is not appropriate here. Because the elements of our dynamic programming tables are in fact parts of the partition function, we can conclude that the numbers involved could have up to O(n) digits. Each arithmetic operation, therefore, becomes much more expensive. While a polynomial-time dynamic programming algorithm based on a context-free grammar is possible, an efficient dynamic programming algorithm would require substantially more work.

Extended Applications

The Markov chain mixing analysis techniques explored in this manuscript have the potential for useful application in a variety of fields. Markov chain Monte Carlo algorithms are widely used in several fields including, machine learning [42], econometrics [43], and Bayesian Statistics [44]. In virtually all applications, an understanding of mixing time increases confidence in the results. In some situations, an understanding of mixing time may also allow for more efficient algorithm selection and implementation.

While many Markov chains with nonuniform stationary distributions have been used for biological applications (e.g., [45-48]), theoretical guarantees on the mixing time are generally not known. Instead, researchers must rely on convergence heuristics, and in fact, many introductions to Markov chain Monte Carlo written for biologists explain such heuristic techniques [49-52]. Of course, heuristics can be misleading, and rigorous mixing time guarantees would be significantly preferable. The same techniques used in this work might be used to generate algorithms with rigorous mixing time bounds for other biological problems concerning a nonuniform distribution.

The mathematical techniques used in this manuscript have been widely used in mathematics, physics, and computer science, demonstrating their broader applicability. For numerous examples, we direct the reader to the books of Levin, Peres, and Wilmer [53]; Montenegro and Tetali [54]; and Jerrum [55].

As an example where similar techniques have found utility in biological applications, it is interesting to briefly consider the study of cladograms, which arise from phylogenetic trees. Mathematically, a cladogram is a binary tree with n labeled leaves and n − 2 unlabeled internal nodes. While an explicit formula is known for the exact number of cladograms of a given size, mixing time under certain dynamics has also been studied. For example, Aldous [56] studied a Markov chain where a leaf is removed at random and then attached to a random edge in the tree, obtaining a proof that the mixing time is bounded below by O(n2) and bounded above by O(n3). Further work by Schweinsberg [57] later proved an upper bound of O(n2), closing the gap between the upper and lower bounds.

Independent Mathematical Research Interests

The plane trees examined as a model for RNA secondary structure are of independent mathematical interest. As Catalan objects, they have been studied combinatorially (see, for example, [25,58]), and Markov chains on Catalan objects have received significant attention over the years [33,34,59-61], but with very few results providing tight estimates on the corresponding mixing times; most commonly these are discussed in the language of Dyck paths. Cohen’s thesis [33] gives an overview of the known mixing time results for chains on Catalan objects. All of the chains surveyed there have a uniform distribution over the Catalan-sized state space as their stationary distribution. Among these, essentially the only known chain with tight bounds (upper and lower bounds differing by a small multiplicative constant) is due to Wilson [34] and gives the relaxation time of O(n3) for the walk consisting of adjacent transpositions on Dyck paths. In comparison, in [59] the chain using all (allowed) transpositions has been shown to have relaxation time of O(n2), and further conjectured to have O(n) as the relaxation time, in analogy with the random transposition shuffle of n cards.

Judging from the lack of progress on several of these chains, it is evident that determining mixing or relaxation time for these chains is typically a challenging problem, even in the case where the stationary distribution is uniform.

In the current work, the RNA secondary structure modeling naturally leads to a state-space on Catalan objects with a nonuniform distribution, making the corresponding mixing time analysis even more challenging. Another example where mixing times are estimated for Markov chains on Catalan objects with nonuniform stationary distribution is the work of Martin and Randall [30], which examines a Gibbs distribution on Dyck paths weighted by the number of returns to the x-axis.

Limitations and Future Directions

While the mixing time proved here is polynomial, it is almost certainly too large to allow for any practical computational sampling experiments. However, we conjecture the actual mixing time to be much smaller, and future work may provide a better bound. Even without additional theoretical results, interesting work is possible using the algorithm we present and heuristic methods for evaluating Markov chain mixing. See ([62], Ch. 8) for a discussion of heuristic methods for monitoring Markov chain convergence.

The results of this study provide an important mathematical foundation for examining the dispersion of RNA secondary structures and branching properties using a Markov chain. However, more work is necessary to optimize the developed computational application for incorporation into the software utilized by biologists that study RNA. Example questions that strongly compel further investigation include:

Can the mixing time bound in our main result be improved?

Is there a rapidly mixing chain, with the same stationary distribution studied here, whose transitions correspond naturally to moves on the set plane trees? Mixing time bounds on the chain of matching exchange moves, as defined in [63], would be especially interesting, as such a chain may relate to RNA folding kinetics.

Is there a rapidly mixing chain converging to the Gibbs distribution using the full energy function for the utilized NNTM model [16]? The chain presented here uses only the parameters α and β, setting γ = 0.

Is there a stochastic context-free grammar which generates secondary structures (in our simplified model or using the full NNTM) according to a Gibbs distribution with NNTM energy?