Predicting folding pathways between RNA conformational structures guided by RNA stacks
 Yuan Li^{1} and
 Shaojie Zhang^{1}Email author
https://doi.org/10.1186/1471210513S3S5
© Li and Zhang; licensee BioMed Central Ltd. 2012
Published: 21 March 2012
Abstract
Background
Accurately predicting low energy barrier folding pathways between conformational secondary structures of an RNA molecule can provide valuable information for understanding its catalytic and regulatory functions. Most existing heuristic algorithms guide the construction of folding pathways by free energies of intermediate structures in the next move during the folding. However due to the size and ruggedness of RNA energy landscape, energyguided search can become trapped in local optima.
Results
In this paper, we propose an algorithm that guides the construction of folding pathways through the formation and destruction of RNA stacks. Guiding the construction of folding pathways by coarse grained movements of RNA stacks can help reduce the search space and make it easier to jump out of local optima. RNAEAPath is able to find lower energy barrier folding pathways between secondary structures of conformational switches and outperforms the existing heuristic algorithms in most test cases.
Conclusions
RNAEAPath provides an alternate approach for predicting lowbarrier folding pathways between RNA conformational secondary structures. The source code of RNAEAPath and the test data sets are available at http://genome.ucf.edu/RNAEAPath.
Keywords
Introduction
RNA molecules play critical roles in the cell. The secondary structures of RNA molecules have been extensively studied because they provide insights into the functionality of RNAs. Native (functional) RNA secondary structures are usually thermodynamically stable and many of them are also the minimum free energy (MFE) structures. Nevertheless, at times, RNA molecules may fold into alternative secondary structures in order to participate in certain biological processes. For example, the SV11 RNA folds into a metastable conformational structure and acts as a template for its own replication using Qβ replicase [1, 2]. Further, RNA conformational switches can transform between alternative secondary structures dynamically in response to various environmental stimuli (such as heat shock and cold shock) [3–6], and carry out RNAmediated biological activities, such as switching on or off downstream gene translation activities [7–9], regulating RNA splicing via multiplestate splicesomal conformations [10], and regulating the life cycles of virus [11].
The conformational transformations between alternative structures involve the folding of an RNA molecule into a series of sequential adjacent intermediate structures [12]. RNA folding pathways provide valuable information for understanding the catalytic and regulatory functions of RNAs (such as hok/sok of plasmid R1 [13]). RNA folding pathways may also impact subsequence biological events (such as formation of tertiary structures). Furthermore, prediction algorithms can help the design of RNA switches by providing prescribed structural alternatives.
In this paper, we present a new approach, RNAEAPath, for computing near optimal direct or indirect folding pathways between two secondary structures of an RNA molecule. We guide the search for low energy barrier folding pathways by integrating a variety of strategies for simulating the formation and destruction of RNA stacks in a flexible framework. Benchmark tests on conformational switches show that RNAEAPath produces lower energy barrier folding pathways and outperforms the existing heuristic approaches in most test cases.
Preliminary
Consider an RNA sequence as a string x = x_{1} ... 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_{ jw }), 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 pseudoknotfree 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_{0}, . . . , 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 thermodynamicallydriven and tends to avoid highenergy 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 NPhard [23]. Many heuristic algorithms have also been proposed following the seminal work of Morgan and Higgs. Flamm et al. [24] used breadthfirst 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 semigreedy 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, semigreedy 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 abovementioned 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 energyguided 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.
Methods
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.
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 Si<{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_{1}, . . . , a_{ m } such that S_{0} = 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 rearranging the order of actions or introducing multiple pairs of complementary actions.
Predicting low energy barrier folding pathways
The initial population of RNAEAPath, ℙ_{0}, 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 ℓ_{1} pathways in ℙ_{ k }_{1} and the best ℓ_{2} 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 ℓ_{2} 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 ℓ_{3} 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 ℓ_{2} offsprings produced by each p ∈ ℙ_{ k }_{1}, no matter how many highqualified 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 ℓ_{1} 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.
Fitness of action chains
The order of folding pathways (valid action chains) is primarily determined by their energy barriers. In case of a tie, the order is determined by the average of energy differences between the initial structure A and intermediate structures. Note that lower energies are preferred in the previous two methods of ordering. If a tie still exists, then shorter action chains are preferred. Action chains are ordered arbitrarily if their relative order can not be determined based on these three criteria.
The initial population of folding pathways
The number of offsprings produced by each mutation strategy
Mutation strategies that have produced more high quality offsprings in the (k  1)^{ st } iteration are allowed to generate more offsprings in the k^{ th } generation. In contrast, mutation strategies that perform poorly in the k  1^{ st } generation, are only allowed to generate a small number ($\mathcal{L}$_{ min }, with default value 3) of offsprings. Note that, the sum of ${\ell}_{{\mathbb{M}}_{y}}^{k}$ for 1 ≤ y ≤ Y may be greater than ℓ.
Stopping conditions
 1.
the energy barrier of OPT_{ k } is equivalent to E(B)  E(A).
 2.
k >γ and the fitness of OPT_{ k } is equivalent to that of OPT_{ k }_{}_{ γ }.
 3.
k ≥ MAX and the fitness of OPT_{ k } is equivalent to that of OPT_{ k }_{1}.
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_{1}, . . . , a_{ m } denote the input pathway A = S_{0}, . . . , 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 1: reordering of actions
As described in the subsection of representation of RNA folding pathways, shuffling the order of actions of the input pathway p can result in a new pathway from A to B. In RNAEAPath, two mutation strategies of this type are employed. ${\mathbb{M}}_{1}$ changes the position of an arbitrary action, and ${\mathbb{M}}_{2}$ swaps the positions of two arbitrary actions.
M_{1}: Let ${\mathbb{M}}_{1}^{{t}_{1},{t}_{2}}\left(p\right)$ denote the sequence of actions obtained by first removing an action ${a}_{{t}_{1}}$ (1 ≤ t_{1} ≤ m) from p and then inserting it after ${a}_{{t}_{2}}$, for all t_{2} ∈ {0,..., t_{1}  1, t_{1} + 1, . . . , m}. Note that the resulting sequence of actions may not necessarily be a valid action chain. For instance, in Figure 1, ${\mathbb{M}}_{1}^{1,4}\left(p\right)={a}_{2},{a}_{3},{a}_{4},{a}_{1},{a}_{5},\dots ,{a}_{8}$ and ${\mathbb{M}}_{1}^{3,2}\left(p\right)=p$ are valid action chains, while ${\mathbb{M}}_{1}^{8,1}\left(p\right)={a}_{1},{a}_{8},{a}_{2},\dots ,{a}_{7}$ is not.
 1.
Choose t_{1} uniformly at random from the interval [1, m].
 2.
Compute the interval [l, u], (t_{1} <l <u <m), where l is the minimum and u is the maximum such that for all t_{2} ∈ [l, u] and t_{2} ≠ t_{1}, ${\mathbb{M}}_{1}^{{t}_{1},{t}_{2}}\left(p\right)$ is a valid action chain.
 3.Choose t_{2} from the interval [l, u].

3.1. If ${a}_{{t}_{1}}$ is an addition operation, for all l ≤ t <t' ≤ u and t ≠ t' ≠ t_{1}, the probability of choosing t is greater than that of t'.

3.2. Otherwise (a deletion operation), for all l ≤ t ≤ t' ≤ u and t ≠ t' ≠ t_{1}, the probability of choosing t is less than that of t'

We do not choose t_{2} (t_{2} ≠ t_{1}) uniformly at random in [l, u], instead, we tend to place addition operations in the front part of p, and deletion operations in the later part of p. This is because adding base pairs early and deleting them late during the folding may help stabilize the intermediate secondary structures. (Please see Additional file 1 for the detailed description of the discrete probability.)
${\mathbb{M}}_{2}$: Let ${\mathbb{M}}_{2}^{{t}_{1},{t}_{2}}\left(p\right)$ denote the sequence of actions obtained by swapping ${a}_{{t}_{1}}$ with ${a}_{{t}_{2}}$. If the resulting sequence of actions is a valid action chain, let it be q; otherwise, restart the process. For example, in Figure 1, ${\mathbb{M}}_{2}^{1,8}\left(p\right)$ is not a valid action chain, while ${\mathbb{M}}_{2}^{2,4}\left(p\right)={a}_{1},{a}_{4},{a}_{3},{a}_{2},{a}_{5},\dots ,{a}_{8}$ is. t_{1} and t_{2} are chosen uniformly at random from {(t_{1}, t_{2}): 1 ≤ t_{1} <t_{2} ≤ m}.
Mutation strategies of type 1 provide methods for shuffling the order of actions of an input pathway and generating slightly different new pathways. However, these strategies are not capable of introducing additional (indirect) base pairs, and the offsprings of a direct pathway produced through type 1 strategies are also direct. In the following, we will describe mutation strategies that are able to construct indirect pathways from a direct pathway.
Type 2: introducing indirect pathways by adding a pair of complementary actions
Morgan and Higgs [16] pointed out that the optimal folding paths are generally indirect pathways. This idea was further described by Dotu et al. [28]. The temporary formation of base pairs, especially those base pairs that do not belong to A ∪ B, may lower the energies of intermediate structures and thus render better folding pathways. Similarly, temporary deletion and reformation of a base pair also can create an indirect pathway.
 1.
Choose t_{1} uniformly at random from the interval [1, m], and obtain the associated intermediate structure ${S}_{{t}_{1}}$.
 2.
Find a set of base pairs that neither conflict with nor clash with ${S}_{{t}_{1}}$ and choose a base pair (i, j) uniformly at random from the set.
 3.
Compute the interval [l, u], (t_{1} <l <u <m), where l is the minimum and u is the maximum such that for all values t_{2} ∈ [l, u] the resulting sequence of actions of ${\mathbb{M}}_{3}^{{t}_{1},{t}_{2},+\left(i,j\right)}\left(p\right)$ is a valid action chain.
 4.
Choose t_{2} from the interval [l, u] with the probability of choosing t greater than that of t' for all t >t'. (This is because (i, j) is not likely to be deleted soon after its formation.)
Mutation strategy ${\mathbb{M}}_{3}$is capable of producing an indirect pathway from a direct pathway. In addition, a proper combination of multiple applications of ${\mathbb{M}}_{3}$ may result in a pathway which simulates the successive formation and deletion of a temporary stack during the folding. Take the pathway p in Figure 1 as an example, we can construct a pathway q that forms a temporary stack consisting of all the GU base pairs via a multiple application of ${\mathbb{M}}_{3},\phantom{\rule{0.5em}{0ex}}q={\mathbb{M}}_{3}^{5,7,+\left(3,14\right)}\left({\mathbb{M}}_{3}^{3,7,+\left(2,15\right)}\left({\mathbb{M}}_{3}^{1,7,+\left(1,16\right)}\left(p\right)\right)\right)$.
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)$.
 1.
Choose t uniformly at random from the interval [1, m], and obtain the associated intermediate structure S_{ t }.
 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.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 }.
 1.
Choose an arbitrary deletion action a_{ t } = del_{ i,j } from p, and obtain the associated intermediate structure S_{ t }.
 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.
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.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}$.

Results and discussion
Benchmark tests
We benchmarked RNAEAPath against existing methods (BARRIEERS [20, 21], PathwayHunter [17], Findpath [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, sbox leader, thiM leader, ms2, HDV, dsrA, ribD leader, amv, alpha operon and HIV1 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].
Benchmarks of BARRIERS, PathwayHunter, Findpath, RNATabuPath, and RNAEAPath for predicting folding pathways between conformational switches on the 18 test cases
Instance  BARRIERS  PathwayHunter  FindPath  RNATabuPath  RNAEAPath  

(n = 1000)  ( n = 1)  ( n = 5)  
rb1      24.04  24.04  23.2  22 
rb2    10  8.2  7.25  6.5  6.5 
rb3      22.4  17.9  17.5  16.7 
rb4      16.9  16.9  16.9  16.9 
rb5      24.54  24.54  21.44  21.44 
hok      28.5  29.66  20.7  20.1 
SL  11.80    13  12.9  13.0  12.9 
attenuator  8.3    8.7  8.6  8.7  8.5 
s15  6.60    7.1  6.6  7.1  7.1 
sbox leader    7.9  5.2  5.2  5.2  5.2 
thiM leader      16.13  14.84  12.3  12.3 
ms2    11.6  6.6  6.6  6.6  6.6 
HDV    23.53  17.4  17.0  16.8  16.8 
dsrA  8.0    8.3  8.2  8.0  8.0 
ribD leader      10.71  9.5  9.5  9.5 
amv    12.2  5.8  5.8  5.74  5.74 
alpha operon    11.8  6.5  6.5  6.1  6.1 
HIV1 leader    14.3  9.3  11.3  8.9  8.9 
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 'pairwiseoptimal' 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 semigreedy modification of the greedy algorithm, a greedy algorithm of Morgan, and Higgs [16] for predicting direct pathways and a variant of the MorganHiggs greedy algorithm capable of producing indirect pathways), that have been shown to perform considerably worse than RNATabuPath [28], are not listed.
Control parameters and performance
In order to evaluate the performance of RNAEAPath with different parameter configurations, we played with several other control parameters, including ℓ_{1}, the number of top offsprings preserved in the next generation, varying from 1 to 16, ℓ_{3}, 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 also evaluated the execution time for each run of RNAEAPath. All the tests were performed on a 32 bit PC with 2.4 GHz Quadprocessor 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.
Conclusions
In conclusion, we have presented a new algorithm, RNAEAPath, for predicting low energy barrier folding pathways between conformational structures. RNAEAPath guides the construction of folding pathways through the destruction and formation of RNA stacks using various types of mutation strategies, and integrates them in a wellestablished computational framework of evolutionary algorithm. These mutation strategies can help reduce the search space and make it easier to jump out of local optima. By analyzing the results, we confirmed that most of the best folding pathways involve the formation of auxiliary stacks, or involve the cooperative formation and disruption of incompatible stacks. The benchmarking results show that RNAEAPath outperforms the existing heuristics on most test cases. We believe that this is because the construction of folding pathways in RNAEAPath captures important biological findings.
Declarations
Acknowledgements
This article has been published as part of BMC Bioinformatics Volume 13 Supplement 3, 2012: ACM Conference on Bioinformatics, Computational Biology and Biomedicine 2011. The full contents of the supplement are available online at http://www.biomedcentral.com/14712105/13/S3.
Authors’ Affiliations
References
 Biebricher CK, Diekmann S, Luce R: Structural analysis of selfreplicating RNA synthesized by Qbeta replicase. J Mol Biol. 1982, 154: 629648. 10.1016/S00222836(82)800193.View ArticlePubMedGoogle Scholar
 Biebricher CK, Luce R: In vitro recombination and terminal elongation of RNA by Q beta replicase. EMBO J. 1992, 11: 51295135.PubMed CentralPubMedGoogle Scholar
 Bukau B: Regulation of the Escherichia coli heatshock response. Mol Microbiol. 1993, 9: 671680. 10.1111/j.13652958.1993.tb01727.x.View ArticlePubMedGoogle Scholar
 Liphardt J, Onoa B, Smith SB, Tinoco I, Bustamante C: Reversible unfolding of single RNA molecules by mechanical force. Science. 2001, 292: 733737. 10.1126/science.1058498.View ArticlePubMedGoogle Scholar
 Narberhaus F: mRNAmediated detection of environmental conditions. Arch Microbiol. 2002, 178: 404410. 10.1007/s0020300204818.View ArticlePubMedGoogle Scholar
 Narberhaus F, Waldminghaus T, Chowdhury S: RNA thermometers. FEMS Microbiol Rev. 2006, 30: 316. 10.1111/j.15746976.2005.004.x.View ArticlePubMedGoogle Scholar
 Montange RK, Batey RT: Riboswitches: emerging themes in RNA structure and function. Annu Rev Biophys. 2008, 37: 117133. 10.1146/annurev.biophys.37.032807.130000.View ArticlePubMedGoogle Scholar
 Roth A, Breaker RR: The structural and functional diversity of metabolitebinding riboswitches. Annu Rev Biochem. 2009, 78: 305334. 10.1146/annurev.biochem.78.070507.135656.View ArticlePubMedGoogle Scholar
 Wakeman CA, Winkler WC, Dann CE: Structural features of metabolitesensing riboswitches. Trends Biochem Sci. 2007, 32: 415424. 10.1016/j.tibs.2007.08.005.PubMed CentralView ArticlePubMedGoogle Scholar
 Smith DJ, Query CC, Konarska MM: Nought may endure but mutability: spliceosome dynamics and the regulation of splicing. Mol Cell. 2008, 30: 657666. 10.1016/j.molcel.2008.04.013.PubMed CentralView ArticlePubMedGoogle Scholar
 Simon AE, Gehrke L: RNA conformational changes in the life cycles of RNA viruses, viroids, and virusassociated RNAs. Biochim Biophys Acta. 2009, 1789: 571583.PubMed CentralView ArticlePubMedGoogle Scholar
 Mahen EM, Watson PY, Cottrell JW, Fedor MJ: mRNA secondary structures fold sequentially but exchange rapidly in vivo. PLoS Biol. 2010, 8: e100030710.1371/journal.pbio.1000307.PubMed CentralView ArticlePubMedGoogle Scholar
 Gerdes K, Gultyaev AP, Franch T, Pedersen K, Mikkelsen ND: Antisense RNAregulated programmed cell death. Annu Rev Genet. 1997, 31: 131. 10.1146/annurev.genet.31.1.1.View ArticlePubMedGoogle Scholar
 Nussinov R, Jacobson AB: Fast algorithm for predicting the secondary structure of singlestranded RNA. Proc Natl Acad Sci USA. 1980, 77: 63096313. 10.1073/pnas.77.11.6309.PubMed CentralView ArticlePubMedGoogle Scholar
 Nussinov R, Pieczenik G, Griggs JR, Kleitman DJ: Algorithms for loop matchings. SIAM Journal on Applied Mathematics. 1978, 35: 6882. 10.1137/0135006.View ArticleGoogle Scholar
 Morgan S, Higgs P: Barrier heights between ground states in a model of RNA secondary structure. J Phys A. 1998, 31 (14): 315310.1088/03054470/31/14/005.View ArticleGoogle Scholar
 Thachuk C, Manuch J, Rafiey A, Mathieson LA, Stacho L, Condon A: An algorithm for the energy barrier problem without pseudoknots and temporary arcs. Pac Symp Biocomput. 2010, 108119.Google Scholar
 Mathews DH, Sabina J, Zuker M, Turner DH: Expanded sequence dependence of thermodynamic parameters improves prediction of RNA secondary structure. J Mol Biol. 1999, 288: 911940. 10.1006/jmbi.1999.2700.View ArticlePubMedGoogle Scholar
 Turner DH, Sugimoto N, Freier SM: RNA structure prediction. Annu Rev Biophys Biophys Chem. 1988, 17: 167192. 10.1146/annurev.bb.17.060188.001123.View ArticlePubMedGoogle Scholar
 Flamm C, Fontana W, Hofacker IL, Schuster P: RNA folding at elementary step resolution. RNA. 2000, 6: 325338. 10.1017/S1355838200992161.PubMed CentralView ArticlePubMedGoogle Scholar
 Flamm C, Hofacker IL, Stadler P, Wolfinger M: Barrier trees of degenerate landscapes. Z Phys Chem. 2002, 216: 155173. 10.1524/zpch.2002.216.2.155.View ArticleGoogle Scholar
 Wuchty S, Fontana W, Hofacker IL, Schuster P: Complete suboptimal folding of RNA and the stability of secondary structures. Biopolymers. 1999, 49: 145165.View ArticlePubMedGoogle Scholar
 Manuch J, Thachuk C, Stacho L, Condon A: NP completeness of the direct energy barrier problem without pseudoknots. 15th International Conference DNA Computing and Molecular Programming. Edited by: Deaton R, Suyama A. 2009, Berlin, Heidelberg: SpringerVerlag, 106115. [http://dx.doi.org/10.1007/9783642106040_11]View ArticleGoogle Scholar
 Flamm C, Hofacker IL, MaurerStroh S, Stadler PF, Zehl M: Design of multistable RNA molecules. RNA. 2001, 7: 254265. 10.1017/S1355838201000863.PubMed CentralView ArticlePubMedGoogle Scholar
 Hofacker IL: Vienna RNA secondary structure server. Nucleic Acids Res. 2003, 31: 34293431. 10.1093/nar/gkg599.PubMed CentralView ArticlePubMedGoogle Scholar
 Voss B, Meyer C, Giegerich R: Evaluating the predictability of conformational switching in RNA. Bioinformatics. 2004, 20: 15731582. 10.1093/bioinformatics/bth129.View ArticlePubMedGoogle Scholar
 Geis M, Flamm C, Wolfinger MT, Tanzer A, Hofacker IL, Middendorf M, Mandl C, Stadler PF, Thurner C: Folding kinetics of large RNAs. J Mol Biol. 2008, 379: 160173. 10.1016/j.jmb.2008.02.064.View ArticlePubMedGoogle Scholar
 Dotu I, Lorenz WA, Van Hentenryck P, Clote P: Computing folding pathways between RNA secondary structures. Nucleic Acids Res. 2010, 38: 17111722. 10.1093/nar/gkp1054.PubMed CentralView ArticlePubMedGoogle Scholar
 Bogomolov S, Mann M, Voß B, Podelski A, Backofen R: Shapebased barrier estimation for RNAs. German Conference on Bioinformatics. 2010, 4150.Google Scholar
 Lorenz R, Flamm C, Hofacker IL: 2D projections of RNA folding landscapes. German Conference on Bioinformatics. 2009, 1120.Google Scholar
 Flamm C, Hofacker IL: Beyond energy minimization: Approaches to the kinetic folding of RNA. Monatsh Chem. 2008, 139: 447457. 10.1007/s0070600808953.View ArticleGoogle Scholar
 Shcherbakova I, Mitra S, Laederach A, Brenowitz M: Energy barriers, pathways, and dynamics during folding of large, multidomain RNAs. Curr Opin Chem Biol. 2008, 12: 655666. 10.1016/j.cbpa.2008.09.017.PubMed CentralView ArticlePubMedGoogle Scholar
 Ndifon W: A complex adaptive systems approach to the kinetic folding of RNA. Biosystems. 2005, 82: 257265. 10.1016/j.biosystems.2005.08.004.View ArticlePubMedGoogle Scholar
 Sponer J, Leszczynski J, Hobza P: Nature of nucleic acidbase stacking: nonempirical ab initio and empirical potential characterizaion of 10 stacked base dimers. Comparison of stacked and Hbonded base pairs. J Phys Chem. 1996, 100 (13): 55905596. 10.1021/jp953306e.View ArticleGoogle Scholar
 Yakovchuk P, Protozanova E, FrankKamenetskii MD: Basestacking and basepairing contributions into thermal stability of the DNA double helix. Nucleic Acids Res. 2006, 34: 564574. 10.1093/nar/gkj454.PubMed CentralView ArticlePubMedGoogle Scholar
 Zhao P, Zhang WB, Chen SJ: Predicting secondary structural folding kinetics for nucleic acids. Biophys J. 2010, 98: 16171625. 10.1016/j.bpj.2009.12.4319.PubMed CentralView ArticlePubMedGoogle Scholar
 Eiben A: Evolutionary computing: the most powerful problem solver in the universe?. Dutch Mathematical Archive. 2002, 5: 126131.Google Scholar
 Bafna V, Tang H, Zhang S: Consensus folding of unaligned RNA sequences revisited. J Comput Biol. 2006, 13: 283295. 10.1089/cmb.2006.13.283.View ArticlePubMedGoogle Scholar
 Mandal M, Boese B, Barrick JE, Winkler WC, Breaker RR: Riboswitches control fundamental biochemical pathways in Bacillus subtilis and other bacteria. Cell. 2003, 113: 577586. 10.1016/S00928674(03)00391X.View ArticlePubMedGoogle Scholar
 paRNAss. [http://bibiserv.techfak.unibielefeld.de/parnass/examples.html]
 Freyhult E, Moulton V, Clote P: Boltzmann probability of RNA structural neighbors and riboswitch detection. Bioinformatics. 2007, 23: 20542062. 10.1093/bioinformatics/btm314.View ArticlePubMedGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.