Predicting folding pathways between RNA conformational structures guided by RNA stacks

Preliminary

Consider an RNA sequence as a string x = x₁ ... x _n of n letters over alphabet ∑ = {A, U, G, C}. A pair of complementary nucleotides x _i and x _j , can form hydrogen bonds and interact with each other, denoted by x _i · x _j . In this paper, we only consider the canonical base pairings (A · U and G · C) and the wobble base pairing (G · U). A secondary structure S of the RNA sequence x is a set of disjoint paired bases (i, j), where 1 ≤ i <j ≤ n. S may be represented by a length n string of dots and brackets, where dots represent unpaired bases and brackets represent paired bases. An RNA structure can comprise of stacks which are lists of consecutive base pairs ({(i, j), (i + 1, j - 1), . . . , (i + w, j - w)} such that x _i · x _j , . . . , x _i+w · x _j-w ), and unstacking base pairs. A secondary structure is pseudoknotted if it contains two base pairs (i, j) and (i', j') with i <i' <j <j'. In this paper, we only consider pseudoknot-free structures. A base pair is compatible with a secondary structure if the base pair can be added to the structure without leading to a pseudoknotted structure or pairing a base with more than one partner. A stack is compatible with S if each base pair in the stack is either in S or is compatible with S.

The free energy of a secondary structure S is denoted by E(S). The set of neighboring structures of S consists of all structures that differ from S by an addition or deletion of exactly one base pair. For two secondary structures A and B, the distance between A and B is the number of base pairs in A not in B plus the number of base pairs in B not in A (i.e. |(A - B) ∪ (B - A)|). A folding pathway from A to B is a sequence of intermediate structures A = S₀, . . . , S _m = B such that for all 0 ≤ i <m, intermediate structure S _{i +1} is a neighboring structure of S _i . A folding pathway is direct if the intermediate structures contain only base pairs in A and B (i.e. S _i ⊆ A ∪ B for 1 ≤ i <m) and otherwise is indirect. The saddle point of a pathway is an intermediate structure with the highest energy, and the energy barrier of a pathway is the energy difference between its saddle point and the initial structure. Since the folding of RNA structures is thermodynamically-driven and tends to avoid high-energy intermediate structures, current computational methods aim to find RNA folding pathways with the lowest energy barriers.

Previous studies

A lot of research has been done on predicting low energy barrier folding pathways. Morgan and Higgs proposed a greedy algorithm that employs the Nussinov model [14, 15] for computing direct folding pathways with minimum energy barrier. They also described a heuristic that samples low energy structures from the partition function and glues them together by direct pathways [16]. The Nussinov model is simple and easy to implement, in which base stacking and loop entropies have no energetic contributions. Based on this model, Thachuk et al. [17] developed an exact algorithm, PathwayHunter, which exploits elegant properties of bipartite graphs for finding the globally optimal direct pathways. However, the Nussinov model is not as accurate as the Turner energy model [18, 19] for approximating RNA thermodynamics. An exact solution based on the Turner energy model is also available. BARRIERS [20, 21], exactly computes the globally optimal folding pathways between any two locally optimal secondary structures. BARRIERS reads an energy sorted list of RNA secondary structural conformations produced by RNAsubopt [22] and is able to compute both direct and indirect low energy barrier pathways.

Nevertheless, the above exact solutions are all exponential in time, because the problem itself is NP-hard [23]. Many heuristic algorithms have also been proposed following the seminal work of Morgan and Higgs. Flamm et al. [24] used breadth-first search in their heuristics (in Vienna RNA Package [25]) and kept the best k candidates at each step to bound the search. Voss et al. [26] devised a straightforward strategy for greedily searching direct pathways. Geis et al. [27] described a greedy heuristic to explore the search space of direct pathways and they also integrated look ahead techniques to diminish the search space. Recently, Dotu et al. [28] developed RNATabuPath, a fast heuristic that employs a TABU semi-greedy search to construct near optimal (both direct and indirect) folding trajectories. In addition, other heuristic approaches, by splitting the pathways into shorter pathways and solving each individually, have also been proposed [29, 30]. There are also other formula presented for the prediction of RNA folding kinetics (see Flamm and Hofacker's review [31] for a systematic discussion).

Many of the existing heuristic algorithms start from an initial structure A, and, at each single step i, walk from the intermediate structure S _i to one of its neighbors S _{i +1} until finally the end structure B is reached. The definition of neighborhood relationships as well as the fitness functions can be different. The fitness function of S _i is usually defined on the free energy of S _i , or the distance from S _i to B, or a function of both. In general, greedy algorithms select the 'best' neighbor structure that has the best fitness. In contrast, semi-greedy algorithms may select any one from the top k structures for randomization. RNATabuPath, which is more sophisticated and outperforms other methods [28], keeps a tabu list for saving recently taken moves such that they can not be applied in certain steps until being removed from the tabu list. In general, during the construction of a folding pathway, these heuristic algorithms select the next intermediate structures from a set of neighboring structures that have the top lowest free energy or have the top shortest distance to B (or the combination of both).

Motivations

However, using energy to guide the construction of folding pathways in the above-mentioned heuristic algorithms has its downsides. The RNA energy landscapes can be extremely large and rugged [11, 32] and the ruggedness of RNA energy landscape may cause the energy-guided search to become trapped in a local optimum. Similar to using structural rearrangements for modeling RNA folding kinetics [33], we want to construct candidate folding pathways in a manner that make it easier to jump out of local optima. It has been revealed that stacking base pairs contribute significantly to the stabilization of RNA secondary structures [34, 35]. The dominant RNA folding pathways involve the formation and destruction of the stacks, and the cooperative formation of a stack along with the partial melting of an incompatible stack [36]. In this paper, we propose to guide the construction of pathways by the formation and destruction of stacks (not by free energy or by distance to the end structure). We still select the constructed folding pathways according to their energy barriers. Although the construction of folding pathways is not driven by thermodynamics, the selection of folding pathways is based on energy barriers. Guiding the construction of folding pathways by coarse grained movements of RNA stacks may help reduce the search space and makes it easier to jump out of local optima. In the rest of this paper, the Methods section describes the representation of folding pathways and the detailed strategies employed by RNAEAPath. The Results and Discussion section presents benchmarking results of RNAEAPath against existing methods followed by concluding remarks in the Conclusions section.

Representation of RNA folding pathways

Given an initial structure A and an end structure B, we use a sequence of actions successively applied to A, rather than a sequence of intermediate structures, to represent a folding pathway from A to B. Representing a pathway by an action chain can avoid cyclic additions and deletions of base pairs and make it easy to simulate the formation and deletion of RNA stacks. A similar representation has also been employed in the previous work of Thachuk et al. [17].

We use two types of actions, add _i,j and del _i,j in the representation of RNA folding pathways. For an intermediate secondary structure S of an RNA sequence x, the action add _i,j denotes the 'add'ition of base pair (i, j) to S (i.e. add _i,j (S) = S ∪ {(i, j)}) and del _i,j denotes the 'del'etion of base pair (i, j) from S (i.e. del _i,j (S) = S - {(i, j)}). An action is direct if it concerns a base pair in A ∪ B and indirect otherwise. The simplest direct pathways from A to B concern sequential deletions of all base pairs in A - B followed by additions of all base pairs in B - A.

Consider an example sequence x = GGGGAAAACCCCUUUU with initial and final structures shown in Figure 1. This simple pathway is obtained by first deleting all GC pairs from A until the RNA is single stranded, and then adding all AU pairs until B is obtained. Note that each intermediate structure S _i differs from both its successor and predecessor by exactly one base pair. The actions in the example are all direct actions and the energy barrier is 5.50 - (-6.60) = 12.10 kcal/mol.

An addition action add _i,j (S) conflicts with S if either x _i or x _j is already paired in S, and it clashes with S if there exists a base pair $\left\{\left(x_i^{\prime },x_j^{\prime }\right)\in S|i<i^{\prime }<j<j^{\prime }\mathsf{\text{ or }}{i}^{\prime }<i<j^{\prime }<j\right\}$. A deletion action del _i,j (S) conflicts with S if (x _i , x _j ) ∉ S. An addition or deletion action is valid and can be applied to S properly if it neither conflicts with nor clashes with S.

A pathway from A to B can be represented by an action chain, which is a sequence of valid actions a₁, . . . , a _m such that S₀ = A, S _t = a _t (S _{t -1}) for 1 ≤ t ≤ m and S _m = B. Note that an action chain for A to B implies a sequence of valid actions that can be successively applied to A without introducing conflicts or clashes and produce B. We use the term "action chain" when the sequence is certified to be valid, and the term "sequence of actions" if its validity is not guaranteed.

This representation of a pathway p from A to B has the following important properties. First, every folding pathway can be represented by a unique action chain and every action chain represents a unique folding pathway (note that it is not necessarily true for a sequence of actions). Second, rearranging the order of actions in p results in a new sequence of actions which represents a new folding pathway from A to B when it is valid. (It is an action chain that can be successively applied to A properly and obtain B.) Third, introducing a pair of complementary actions (e.g. add _i,j and del _i,j ) to p results in a new sequence of actions which also represents a new folding pathway from A to B if it is valid.

In RNAEAPath, folding pathways are represented in the form of action chains, instead of a sequence of intermediate structures. This representation makes the life cycle of a folding pathway transparent to the algorithm and also makes it easier for us to simulate the cooperative formation and destruction of RNA stacks by re-arranging the order of actions or introducing multiple pairs of complementary actions.

Predicting low energy barrier folding pathways

Given an RNA sequence x, an initial structure A and a final structure B, RNAEAPath computes a near optimal low energy barrier folding pathway from A to B in an evolutionary algorithm framework [37]. Figure 2 elucidates the overall paradigm for RNAEAPath. In this algorithm, the population of each generation is comprised of folding pathways ordered by their fitness. The functions ${\mathbb{M}}_y\left(p\right)$ are mutation strategies, each of which takes in a pathway p and produces a set of offspring pathways. These mutation strategies are central to the effectiveness of RNAEAPath and will be discussed in the Mutation strategies subsection. ℓ₁, ℓ₂, ℓ₃, MAX and γ are positive integer control parameters.

The initial population of RNAEAPath, ℙ₀, is filled with a set of simple pathways. Then, the algorithm goes through several iterations. ℙ _{k -1} is the population of the k - 1 ^st iteration. In the k ^th iteration, the algorithm produces ${\mathbb{O}}_k$ (an ordered list of pathways) and ℙ _k (the population of the k ^th iteration) from ℙ _{k -1}. ${\mathbb{O}}_k$ stores the best ℓ₁ pathways in ℙ _{k -1} and the best ℓ₂ pathways produced by each p ∈ ℙ _{k -1}. More specifically, each pathway p ∈ ℙ _{k -1} produces $t_y^k$ offsprings through every mutation strategy ${\mathbb{M}}_y\left(1\le y\le Y\right)$. The resulting offsprings produced by p are stored in a temporary list $\mathbb{T}$, and the top ℓ₂ pathways are added to ${\mathbb{O}}_k$. Finally, the best solution of the k ^th iteration, termed as OPT _k , is the best pathway in ${\mathbb{O}}_k$. And, ℙ _k (the population of the k ^th iteration) is composed of the best ℓ₃ pathways of ${\mathbb{O}}_k$ and will be used in the next iteration to produce ℙ _{k +1}. This helps keep the diversity of the population large, since ℙ _k contains at most ℓ₂ offsprings produced by each p ∈ ℙ _{k -1}, no matter how many high-qualified offsprings are produced by each pathway. The algorithm terminates when a stopping condition is met, and it returns the best solution of the last iteration. Since ${\mathbb{O}}_k$ retains the best ℓ₁ pathways from ℙ _{k -1} in each iteration, the best one ever encountered by the algorithm is retained in lists ${\mathbb{O}}_k$ and ℙ _k , and stored in OPT _k . So, OPT _k has no worse fitness when compared to OPT _{k -1}, and RNAEAPath always returns the best action chain it ever discovered.

In the remaining of this section, we discuss details regarding fitness evaluation, initialization of the population, stopping conditions and mutation strategies of RNAEAPath.

The initial population of folding pathways

The initial population, ℙ₀, contains 4 simple pathways from A to B formed by first deleting all base pairs in A - B and then adding those in B - A, similar to the pathway shown in Figure 1. Although we can also arrange base pair deletions and additions in an arbitrary order, we tailor them in a manner that simulates successive degradation and formation of RNA stacks. This is because random deletions and additions of base pairs tend to form additional unpaired loop regions that introduce entropic penalties (see Figure 3 for an illustration). We can degrade or form each stack either from the outmost base pair to the innermost base pair or vice verse. Usually, it yields a lower energy barrier if we degrade a stack from the outmost base pair to the innermost base pair and form a stack from the innermost base pair to the outmost base pair. However, for the sake of simplicity and generosity, we construct 4 simple pathways in ℙ₀, which degrade all the stacks from the same direction and form all the stacks from the same direction. These simple pathways constitute a diversified and unbiased initial population for the algorithm to start from.

Mutation strategies

In RNAEAPath, the mutation strategies employed to evolve folding pathways can be categorized into three types: (1) rearranging the order of actions, (2) introducing indirect pathways and (3) formation of a single stack or cooperative conversion of a pair of incompatible stacks. In this section, let ${\mathbb{M}}_1,\dots ,{\mathbb{M}}_Y$ denote the mutation strategies and let p = a₁, . . . , a _m denote the input pathway A = S₀, . . . , S _m = B. For each mutation strategy ${\mathbb{M}}_y\left(p\right)$, we describe the process for generating one new pathway q using each mutation strategy when given p.

Type 3: formation of a single stack or simultaneous formation and deletion of a pair of incompatible stacks

In this section, we will introduce mutation strategies for producing pathways that involve with formation and deletion of stacks. To perform this type of strategies, we first need to find all possible stacks in an RNA sequence x. We use the algorithm of Bafna et al. [38] to find the set of all possible stacks with more than 3 consecutive base pairs, and denote it by STA(x). There are two strategies in Type 3: formation of a single stack $\left({\mathbb{M}}_4\right)$ and simultaneous formation and destruction of a pair of incompatible stacks $\left({\mathbb{M}}_5\right)$.

${\mathbb{M}}_4$: Let ${\mathbb{M}}_4^{t,h}\left(p\right)$ denote the sequence of actions obtained by forcing the formation of a stack stack _h ∈ STA after action a _t , where stack _h is compatible with S _t . The following describes the procedure for computing ${\mathbb{M}}_4^{t,h}\left(p\right)$.

1. 1.

   Choose t uniformly at random from the interval [1, m], and obtain the associated intermediate structure S _t .

    
2. 2.

   Find a set of stacks that neither conflict with nor clash with S _t , and pick up a stack stack _h uniformly at random from the set.

    
3. 3.
   Ensure that each base pair (i, j) in {stack _h - S _t } is sequentially (from the innermost base pair to the outmost base pair) formed after a _t .

   • 3.1. If an action add _i,j appears in {a_t+1, . . . , a _m }, move it up and place it after a _t using strategy ${\mathbb{M}}_1$.

   • 3.2. Otherwise, introduce a pair of complementary actions add _i,j and del _i,j to p after a _t using strategy ${\mathbb{M}}_3$.

    

We can introduce additional stacks that are compatible with S _t using ${\mathbb{M}}_4$ by forcing a sequence of addition actions successively forming base pairs in {stack _h - S _t }, after a _t .

${\mathbb{M}}_5$: Let ${\mathbb{M}}_5^{t,h}\left(p\right)$ denote the sequence of actions obtained by forcing the formation of a stack stack _h ∈ STA which is incompatible with S _t , after action a _t . Shown on the right side of Figure 4 is a folding pathway which simultaneously destructs and forms a pair of incompatible stacks. Shown on the left side is a simple folding pathway which has exactly the same start and end structures, while it folds into a single stranded structure during the folding. Usually, the pathway on the right has lower energy barrier than the one on the left because it never folds into a single stranded structure. The folding pathway on the right side of Figure 4 can be introduced using strategy ${\mathbb{M}}_5$. And, the procedure for computing ${\mathbb{M}}_5^{t,h}\left(p\right)$ is as follows:

1. 1.

   Choose an arbitrary deletion action a _t = del _i,j from p, and obtain the associated intermediate structure S _t .

    
2. 2.

   Find a set of stacks which either conflicts with or clashes with S _t , and choose a stack stack _h uniformly at random from the set.

    
3. 3.

   For each base pair (i', j') in {stack _h - S _t } that is compatible with S _t , place add _{i ', j '} to p after a _t using strategy ${\mathbb{M}}_4$.

    
4. 4.
   For each base pair (i', j') in {stack _h - S _t } that is incompatible with S _t ,

   • 4.1. Find all the base pairs (i*, j*) in S _t that are incompatible with (i', j'), and ensure that each base pair (i*, j*) is deleted before the action add _{i ', j'} .

   • 4.2. If a action $\mathsf{\text{de}}{\mathsf{\text{l}}}_{i^{*},j^{*}}$ appears in {a _{t +1}, . . . , a _m }, move it up before add _{i '}, _{j '}using strategy ${\mathbb{M}}_1$.

   • 4.3. Otherwise, introduce a pair of complementary actions $\mathsf{\text{de}}{\mathsf{\text{l}}}_{i^{*},j^{*}}$ and $\mathsf{\text{ad}}{\mathsf{\text{d}}}_{i^{*},j^{*}}$using strategy ${\mathbb{M}}_3$.

    

Using ${\mathbb{M}}_5$, we can introduce the simultaneous formation of a stack stack _h , which is incompatible with S _t , and destruction of existent stacks (or base pairs) that hamper the formation of stack _h . Since cooperative formation and destruction of stacks may contribute additional stacking energies for stabilizing the intermediate structures, better folding pathways with lower energy barriers may be rendered.

Benchmark tests

We benchmarked RNAEAPath against existing methods (BARRIEERS [20, 21], PathwayHunter [17], Find-path [24], and RNATabuPath [28]) by predicting low energy barrier folding pathways between two designated RNA secondary structures of 18 conformational switches. All the conformational switches were taken from the work of Dotu et. al [28]. Five of them are riboswitches, including rb1, rb2, rb3, rb4, and rb5. The metastable structures of these riboswitches have been experimentally determined by inline probing [9, 39]. The thirteen remaining cases concern conformational switches, including hok, SL (Spliced leader RNA), s15, s-box leader, thiM leader, ms2, HDV, dsrA, ribD leader, amv, alpha operon and HIV-1 leader. Sequences of these conformational switches can also be obtained from paRNAss web site [40], and some of the metastable secondary structures were computationally determined using RNAbor [41].

BARRIERS is the only exact solution that produces indirect pathways based on the Turner model. BARRIERS has already been compared with existing heuristic algorithms on the same test cases in the work of Dotu et al. [28]. We put the results of BARRIERS in the table just for the sake of comparison. It has been pointed out that BARRIERS gives provably globally optimal pathways in 4 out of 18 cases (i.e. SL, attenuator, s15 and dsrA). BARRIERS can not be directly applied to 5 cases because either the initial or the end structure is not locally optimal (i.e. rb2, sbox leader, ms2, amv and alpha operon), and can not converge in the remaining cases. Possibly due to the fact that both the number of RNA secondary conformations to consider and the computational resources required increase exponentially with the growing length of the RNA sequence and the growing range of energy barrier. PathwayHunter is an exact algorithm capable of producing the optimal direct folding pathways based on the Nussinov model. PathwayHunter can not be directly applied to 10 cases, because it requires the pair of input structures being able to form a 'pairwise-optimal' bipartite conflicting graph (see the work of Thachuk et al. [17] for details). It is not surprising that the performance of the exact algorithm, PathwayHunter, evaluated by free energy (in kcal/mol), is worse than the heuristic algorithms. This is because PathwayHunter is optimized based on the Nussinov model and only produces direct pathways, while the optimal direct pathways predicted based on the Nussinov model may not be the optimal pathways (considering both direct and indirect pathways) based on the Turner model. All the remaining three methods are heuristics capable of producing both direct and indirect pathways based on the Turner model. Findpath produces folding pathways very quickly, however it performs worse than both RNATabuPath and RNAEAPath in most cases. RNATabuPath performs better than Findpath, but produces less optimal pathways than RNAEAPath. The energy barriers predicted by RNAEAPath over 5 runs are exactly the same as RNATabuPath in 5 cases, worse in 1 case, and better in all the remaining 12 cases.

Other heuristic algorithms (including a greedy algorithm of Voss et al. [26], a semi-greedy modification of the greedy algorithm, a greedy algorithm of Morgan, and Higgs [16] for predicting direct pathways and a variant of the Morgan-Higgs greedy algorithm capable of producing indirect pathways), that have been shown to perform considerably worse than RNATabuPath [28], are not listed.

By analyzing the best folding pathways produced by RNAEAPath, we found that most high-quality pathways involve the melting of stacks in the initial structure, the (possibly simultaneous) construction of stacks in the final structure, and the formation of auxiliary temporary stacks for obtaining folding pathways with lower energy barriers. We may take the lowest energy barrier folding pathway of rb2 found by RNAEAPath, shown in Figure 5 as an example. The stack colored in red is an auxiliary temporary stack introducing intermediate structures with lower free energies (which is constructed using ${\mathbb{M}}_4$). Some of the stacks in the initial structure (in blue) are gradually melting, while at the same time, an incompatible stack (in green) is being formed (which is constructed using ${\mathbb{M}}_5$). The stack colored in red is an auxiliary temporary stack introducing intermediate structures with lower free energies. This example convinces us that the advantages of RNAEAPath mainly come from employing mutation strategies that guide the construction of folding pathways by the formation and destruction of stacks and introducing additional stacking interactions that are important for stabilizing the intermediate structures. Detailed low energy barrier folding pathways for all the test cases are available on RNAEAPath web site.

Control parameters and performance

In order to evaluate the performance of RNAEAPath with different parameter configurations, we played with several other control parameters, including ℓ₁, the number of top offsprings preserved in the next generation, varying from 1 to 16, ℓ₃, the size of population in each generation, varying from 80 to 120 and $\mathcal{L}$, the total number of offsprings each individual is expected to produce, varying from 80 to 120. The detailed results are shown in Additional file 1. In general, RNAEAPath produces pathways of roughly the same quality for most test cases with different control parameters, among which the default parameter setting is the best.

We explored the relationship between the performance of RNAEAPath and the number of generations completed by plotting energy barriers of the best folding pathways produced by RNAEAPath with the default parameters in each generation, as shown in Figure 6. In general, the energy barriers decrease dramatically in the first one or two generations, and then the decrements slow down and finally plateau within 10 generations. For instance, in the case of rb3, the predicted energy barriers of folding pathways in the initial population is 27.3 kcal/mol. It decreases by 7.2 kcal/mol (24.9%) through the first two generations and decreases by 2.5 kcal/mol (9.2%) through the next three generations. Through all the remaining generations, no further improvement is made.

We also evaluated the execution time for each run of RNAEAPath. All the tests were performed on a 32 bit PC with 2.4 GHz Quad-processor and 3.2 GB memory, running Fedora 11. With the default control parameters, RNAEAPath terminates in 1 minute in the best case (rb4), 445 minutes in the worst case (hok), and 43 minutes on average. The detailed running times are shown in Additional file 1. We did not perform direct comparisons between the running time of RNATabuPath and that of RNAEAPath, since RNATabuPath is only accessible via web server.