in this article , we describe a new computational tool to determine nearly optimal folding pathways between two given secondary structures of an rna sequence .
our tool , rnatabupath , and related web server , rnapathfinder , have potential applications in synthetic biology ; in particular , our work can be used to help engineer bistable conformational switches with reasonable folding kinetics [ see abfalter et al .
( 1 ) and flamm et al . ( 2 ) for methods to computationally design bistable switches . ) .
folding pathways play an important role in various biological processes , including the hok / sok ( host killing / suppression of killing ] system ( 3 ) and transition between two metastable structures , as in the conformational switch in spliced leader ( sl ) rna from leptomonas collosoma ( 4 ) . in the hok / sok system ,
the hok gene of escherichia coli codes a small ( 52 amino acids ) toxin causing irreversible damage to the cell membrane .
while the very stable hok - mrna is constitutively expressed from a weak promoter , the highly unstable ( rapidly degraded ) sok - rna is constitutively expressed from a strong promoter .
dalgarno sequence ; however , slow exonucleolytic processing digests the last 40 nt of the 3-end of hok - mrna , thus transforming the molecule into its active form in which the shine 
dalgarno sequence is no longer sequestered . if r1 plasmids of e. coli are present in sufficient copy number , then a portion of the 64 nt sok - rna , which is complementary to hok - mrna leader region , binds to the active conformation of hok - mrna , thus causing degradation of the complex by rnaseiii ( 3 ) . if plasmids are not present in sufficient copy number , then the cell is killed by hok toxin . in this fashion , efficient plasmid stabilization is ensured in the population .
[ see ( 3 ) for a review of the hok / sok system . ] in the case of sl rna from certain trypanosomes and nematodes , a portion of the 5 exon is donated to another mrna by trans splicing .
intermediate structures may be important for the process of splicing , as shown by lecuyer and crothers ( 4 ) , who performed stopped - flow rapid - mixing and temperature - jump measurements of the kinetics for the structural transition between two low energy structures of sl rna from l. collosoma .
conformational switches are thought not only to play a role in such trans splicing but also in transcriptional and translational regulation , protein synthesis and mrna splicing . as indicated by the examples of hok / sok and sl rna ,
for that reason , this problem has been considered by a number of authors , both in the context of rna secondary and tertiary structure .
nudge elastic band ( neb ) method in amber to sample low energy paths for rna conformational changes at the three - dimensional atomic scale .
 crick / hoogsteen gg non - canonical pairs , where one g is syn around the glycosidic bond while the other g is anti .
since prior nmr - constrained modeling had demonstrated that the gg pairs change from ( syn)g-(anti)g to ( anti)g-(syn)g on the millisecond timescale , such atomic - level simulations using amber were feasible . due to large structural transitions between metastable structures in conformational switches , it seems clear that three - dimensional atomic scale simulations using molecular dynamics can not adequately address the general problem posed in this article .
for that reason , it is important to develop efficient algorithms to determine optimal and suboptimal folding pathways between rna secondary structures .
intermediate structures from such low energy pathways can then be further investigated using atomic scale methods such as neb .
morgan and higgs ( 6 ) appear to be among the first to have considered the problem of determining an optimal  folding pathway between two given secondary structures a , b of a given rna sequence .
if a , b are secondary structures for a given rna sequence s , then a  folding pathway from a to b is a sequence a=0 , 1 ,  , n = b such that each intermediate structure i differs from the next structure i+1 by exactly one base pair .
direct if every intermediate structure i+1 is obtained from the preceding structure i by either adding a base pair that belongs to b but not a or by removing a base pair that belongs to a but not b. if a pathway is not direct , then it is 
indirect. the  saddle point in a pathway a = 0 , 1 ,  , n = b is the intermediate structure i of highest energy ( in case there is more than one intermediate structure having maximum energy along the path , we define the saddle point to be the first such structure , having smallest index . ) .
barrier energy of a pathway from a to b is the energy difference e( )  e(a ) , where  is the saddle point of the pathway . clearly , the barrier energy is of fundamental importance in folding kinetics .
morgan and higgs describe both a  greedy algorithm to construct a direct pathway , as described in the  materials and methods section , as well as an algorithm to construct an indirect pathway by gluing together greedy direct pathways between low energy structures sampled from the partition function ( in ( 6 ) , morgan and higgs compute the partition function z = s exp(e(s)/rt ) , where the sum is over all secondary structures of a given rna sequence , and e(s ) is the nussinov energy model ( 7 ) .
since the partition function is inductively computed , it is simple to the stochastically sampled structures from the low energy boltzmann ensemble .
later , ding and lawrence ( 8) describe the same stochastic sampling algorithm , sfold , with the exception that turner energy model ( 9 ) is used in the place of the nussinov energy model . ) .
while morgan and higgs had worked with the nussinov energy model ( 7 ) , which ascribes 1 per base pair , with no energetic contribution due to base stacking or loop entropies , our implementation of the direct and indirect pathway algorithms of morgan and higgs uses the turner nearest neighbor energy model ( 914 ) , whose parameters have been obtained by uv absorption ( optical melting ) experiments . since the pioneering work of morgan and higgs , other groups have developed methods to compute folding pathways between secondary structures .
flamm et al . ( 2,15 ) describe an exact algorithm , barriers ( barriers is available at http://www.tbi.univie.ac.at/ivo / rna / barriers/. ) , that computes optimal ( possibly indirect ) folding pathways between any two locally optimal secondary structures ( a locally optimal secondary structure is one in which the energy is not lowered if a single base pair is either added or removed . sometimes such structures are called 
metastable ) . while most biologically important examples of pathway computation concern metastable or locally optimal structures , there are nevertheless important exceptions , such as conformational switches ( incompletely ) determined by experimental methods ; the adenine riboswitch from vibrio vulnificus ( rb2 ) ( 16 ) is indeed one such example .
barriers relies on the vienna rna package program rnasubopt ( 17 ) that exhaustively generates all secondary structures within a user - specified energy upper bound . for this reason , although barriers is the only exact algorithm , it is generally limited to relatively small rna sequences or those for which the energy of the saddle point between a and b is not too large . in ( 18 ) , flamm et al . describe a breadth - first search algorithm with bounded look - ahead , to compute nearly optimal
the algorithm is implemented in the program findpath.c , now part of the vienna rna package . finally , as part of the method parnass , voss et al .
our new algorithm , rnatabupath , produces ( possibly indirect ) almost optimal folding pathways by using a heuristic from combinatorial optimization theory known as 
tabu search. tabu search , described in the text by f. w. glover and m. laguna ( 20 ) , is a meta - heuristic to avoid being trapped in local optima in local search algorithms .
one of its key components , which we use in this article , is a short memory , often called the 
tabu search then selects the best configuration in the neighborhood which is not in the tabu list . this neighbor may in fact degrade the value of the objective function
tabu search has been a very effective technique in combinatorial optimization for a wide variety of problems and is an integral part of the repertoire of optimization techniques . to fix ideas , figure 1
depicts three folding pathways for a toy 12 nt rna sequence ggggggcccccc , with structures a = .
( ( ..... ) ) .. having free energy of 1.40 kcal / mol and b = .. ( ( ( ... ) ) ) . having free energy of 1.70 kcal / mol .
structure b is not locally optimal , since by adding the base pair ( 1 , 11 ) to a and by adding base pair ( 2 , 12 ) , one obtains structures ( ( ( ..... ) ) ) . and .
( ( ( ( ... ) ) ) ) having free energies 4.70 and 4.20 kcal / mol , respectively .
it follows that barriers can not be applied ( in such cases , following the suggestion of an anonymous referee , one could first determine locally optimal structures  ,  that , respectively , contain a , b , apply barriers to find an optimal path between  , . this yields a near - optimal path between a , b. ) . the left
, middle and right panels displays the path computed by our implementation of the morgan 
figure 1.three folding pathways for the ( toy ) rna sequence s = ggggggcccccc , between the secondary structure a = .
( ( ..... ) ) .. with free energy 1.40 kcal / mol and the structure b = .. ( ( ( ... ) ) ) . with 1.70 kcal / mol .
the left panel of this figure depicts a ( direct ) folding pathway from a to b produced by our implementation of the morgan 
the middle panel depicts the indirect folding pathway produced by our implementation of the extension of morgan 
( ......... ) contains the base pair ( 2 , 12 ) which is present in neither a nor b. the right panel depicts a folding pathway from a to b produced by our rnatabupath algorithm .
although rnatabupath often yields indirect pathways , in this case , the pathway returned by rnatabupath is direct .
note that the last three structures proposed by rnatabupath are ... ( ..... ) .. , .. ( ( ..... ) ) .
loop is energetically more favorable than the alternative ( not proposed by rnatabupath ) , given by structures .... ( ... ) ... , ... ( ( ... ) ) .. , .. ( ( ( ... ) ) ) .
, respectively , having free energy of 4.90 , 1.60 and 1.70 kcal / mol .
secondary structures are indicated in the familiar ( vienna ) dot bracket notation , while free energy in kcal / mol appears to the right of each structure .
free energies are determined by the program rnaeval from the vienna rna package ( 27 ) .
three folding pathways for the ( toy ) rna sequence s = ggggggcccccc , between the secondary structure a = .
( ( ..... ) ) .. with free energy 1.40 kcal / mol and the structure b = .. ( ( ( ... ) ) ) . with 1.70 kcal / mol .
the left panel of this figure depicts a ( direct ) folding pathway from a to b produced by our implementation of the morgan 
the middle panel depicts the indirect folding pathway produced by our implementation of the extension of morgan 
( ......... ) contains the base pair ( 2 , 12 ) which is present in neither a nor b. the right panel depicts a folding pathway from a to b produced by our rnatabupath algorithm .
although rnatabupath often yields indirect pathways , in this case , the pathway returned by rnatabupath is direct .
note that the last three structures proposed by rnatabupath are ... ( ..... ) .. , .. ( ( ..... ) ) .
loop is energetically more favorable than the alternative ( not proposed by rnatabupath ) , given by structures .... ( ... ) ... , ... ( ( ... ) ) .. , .. ( ( ( ... ) ) ) .
, respectively , having free energy of 4.90 , 1.60 and 1.70 kcal / mol .
secondary structures are indicated in the familiar ( vienna ) dot bracket notation , while free energy in kcal / mol appears to the right of each structure .
free energies are determined by the program rnaeval from the vienna rna package ( 27 ) .
figures 24 depict examples where indirect pathways may have ( provably ) lower barrier energies than every direct pathway , while figure 5 displays the two meta - stable structures of host killing ( hok ) rna .
pseudobase ( 21 ) , furnish canonical examples where direct pathways are likely to have greater barrier energies than even naive indirect pathways .
type - h pseudoknots admit a planar representation where certain base pairs are depicted above the horizontal line corresponding to the rna sequence , while others are depicted below the line  see figure 2 for illustration .
define structure a [ respectively b ] to consist of those base pairs above [ respectively below ] the line .
clearly any direct path from a to b must proceed by removal of all base pairs from a , resulting in the empty structures , followed by addition of all base pairs from b. it follows that e(a ) is a lower bound for the barrier energy of every direct path from a to b , where a , b are indicated in figure 2 . 

figure 2.consider the 46 nt rna sequence cgcgacggcu acgcgacggc aaugccguug cgaagccguc gcgauc , with secondary structures a = ( ( ( ( ( ( ( ( ( .............. ) ) ) ) ) ) ) ) ) .............. having free energy 16.04 kcal / mol and b = ........... ( ( ( ( ( ( ( ( ( .............. ) ) ) ) ) ) ) ) ) ... having free energy 18.14 kcal / mol .
the structure a consists of the base pairs lying above the line in this figure , while the structure b consists of the base pairs lying below the line .
program barriers can not be used , since neither a nor b is locally optimal . 

figure 3.a manually designed indirect folding pathway for the 46 nt rna sequence cgcgacggcu acgcgacggc aaugccguug cgaagccguc gcgauc , proceeding from locally optimal secondary structure a = ( ( ( ( ( ( ( ( ( .............. ) ) ) ) ) ) ) ) ) .............. having free energy 16.04 kcal / mol to locally optimal structure b = ........... ( ( ( ( ( ( ( ( ( .............. ) ) ) ) ) ) ) ) ) ... having free energy 18.14 kcal / mol . intuitively , this pathway can be visualized as repeatedly moving the remaining rightmost right - parenthesis to the right , then repeatedly moving the rightmost left - parenthesis to the right . in this manner ,
the barrier energy of this indirect path is 13.68 kcal / mol , while every direct path must have a barrier energy of at least 16.04 kcal / mol , since the empty structure must be an intermediate structure in every direct path in this example . indeed , due to nucleation energy required to start a hairpin in the empty structure , the barrier energy of every direct path must properly exceed 16.04 . in this case , vienna package program findpath.c with look - ahead 100 returns a barrier energy of 18.27 , while rnatabupath returns a barrier energy of 16.84 . 

figure 4.this figure depicts a folding intermediate in a low energy  indirect path from a to b ( unexplained notation taken from figure 2 ) .
clearly , every direct path from a to b must have the empty structure as an intermediate structure , hence the lowest barrier energy of a direct folding pathway must be at least 16.04 kcal / mol ( in fact even larger due to nucleation energy ) .
however the indirect folding pathway depicted in table 3 has a barrier energy of 13.68 kcal / mol . 

figure 5.two secondary structures of host killing ( hok ) rna , taken from figure 8 of shapiro et al .
the left panel depicts the secondary structure of 396 nt hok - rna , presumably based on figure 1b of franch et al .
the right panel depicts the secondary structure of the 361 nt truncated hok - rna , after 3 processing .
since rnatabupath requires two structures a , b , of the same length for a given rna sequence , we have extended the secondary structure of truncated hok - rna to consist of unpaired nucleotides .
free energy of tructure a is 186.1 kcal / mol , while that of structure b is 142.8 kcal / mol . in this case , findpath ( 18 ) obtained the best barrier energy . consider the 46 nt rna sequence cgcgacggcu acgcgacggc aaugccguug cgaagccguc gcgauc , with secondary structures a = ( ( ( ( ( ( ( ( ( .............. ) ) ) ) ) ) ) ) ) .............. having free energy 16.04 kcal / mol and b = ........... ( ( ( ( ( ( ( ( ( .............. ) ) ) ) ) ) ) ) ) ... having free energy 18.14 kcal / mol .
the structure a consists of the base pairs lying above the line in this figure , while the structure b consists of the base pairs lying below the line .
program barriers can not be used , since neither a nor b is locally optimal . a manually designed indirect folding pathway for the 46 nt rna sequence cgcgacggcu acgcgacggc aaugccguug cgaagccguc gcgauc , proceeding from locally optimal secondary structure a = ( ( ( ( ( ( ( ( ( .............. ) ) ) ) ) ) ) ) ) .............. having free energy 16.04 kcal / mol to locally optimal structure b = ........... ( ( ( ( ( ( ( ( ( .............. ) ) ) ) ) ) ) ) ) ... having free energy 18.14 kcal / mol .
intuitively , this pathway can be visualized as repeatedly moving the remaining rightmost right - parenthesis to the right , then repeatedly moving the rightmost left - parenthesis to the right . in this manner ,
the barrier energy of this indirect path is 13.68 kcal / mol , while every direct path must have a barrier energy of at least 16.04 kcal / mol , since the empty structure must be an intermediate structure in every direct path in this example .
indeed , due to nucleation energy required to start a hairpin in the empty structure , the barrier energy of every direct path must properly exceed 16.04 . in this case , vienna package program findpath.c with look - ahead 100 returns a barrier energy of 18.27 , while rnatabupath returns a barrier energy of 16.84 .
clearly , every direct path from a to b must have the empty structure as an intermediate structure , hence the lowest barrier energy of a direct folding pathway must be at least 16.04 kcal / mol ( in fact even larger due to nucleation energy )
. however the indirect folding pathway depicted in table 3 has a barrier energy of 13.68 kcal / mol .
two secondary structures of host killing ( hok ) rna , taken from figure 8 of shapiro et al .
the left panel depicts the secondary structure of 396 nt hok - rna , presumably based on figure 1b of franch et al .
the right panel depicts the secondary structure of the 361 nt truncated hok - rna , after 3 processing .
since rnatabupath requires two structures a , b , of the same length for a given rna sequence , we have extended the secondary structure of truncated hok - rna to consist of unpaired nucleotides .
free energy of tructure a is 186.1 kcal / mol , while that of structure b is 142.8 kcal / mol . in this case , findpath ( 18 ) obtained the best barrier energy . given the combinatorial difficulty of determining optimal pathways for the turner energy model and the inherent exponential time complexity of the program barriers , it is perhaps not surprising that the problem of computing the minimum energy path between two given rna structures has recently been announced to be np - complete .
the np - completeness of computing an optimal pathway is proven in the pre - print , manuch , j . ,
( 2009 )  np - completeness of the direct energy barrier problem without pseudoknots , in 15th international meeting on dna computing and molecular programming , june 811 , fayetteville , arkansas .
in this section , we survey several known heuristics for determining folding pathways between two secondary structures , as well as present our novel semi - greedy and rnatabupath methods . to explain the morgan 
higgs greedy direct pathway algorithm , we first define the notion of a base pair 
clashing with another base pair  base pair ( i , j ) is said to clash with base pair ( x , y ) if either x  i  y  j or i  x  j  y. more generally , a base pair ( i , j ) clashes with a secondary structure a if there exists ( x , y )  a such that ( i , j ) clashes with ( x , y ) .
the set of base pairs ( x , y )  a such that ( i , j ) clashes with ( x , y ) is denoted clash(i , j , a ) ; i.e. 

 with this definition , the morgan  higgs greedy algorithm repeatedly performs the following steps : ( i ) determine the base pair ( i , j ) belonging to b but not a which has minimum size clash set c , ( ii ) remove base pairs from c , and ( iii ) add base pairs in b that do not induce any new clashes .
figure 6.morganhiggs greedy algorithm ( 6 ) to construct a greedy direct pathway from secondary structure a to b. morgan 
higgs greedy algorithm ( 6 ) to construct a greedy direct pathway from secondary structure a to b. the morgan 
higgs algorithm to compute a nearly optimal ( possibly ) indirect pathway between secondary structures a , b proceeds as follows . by sampling ,
if either a or b does not belong to  , then add the missing structure to . define a complete , weighted , undirected graph g = ( v , e ) , where the set v of vertices consists of all structures in  , and the edge weight between any two structures  ,  is defined to be energy barrier max{e(i )  e( ) : 1  i  n } , where  = 0 ,  , n =  is the greedy direct pathway from  to  , as determined by the morgan 
morgan and higgs then apply  single link cluster ( slc ) algorithm , as described in ( 22 ) , in order to determine an optimal pathway , starting from structure a , proceeding by hopping from one low energy structure in  to another via a greedy direct pathway , and terminating by the structure b. in our implementation of the morgan 
higgs indirect algorithm , we sample low energy structures with respect to the turner energy model by applying the ding  lawrence algorithm ( 8) , as implemented in rnasubopt - p from the vienna rna package ( step 2 of figure 7 ) . in place of the slc algorithm , we apply a modified form of dijkstra 's single source shortest path algorithm ( 23 ) , in order to determine a sequence a=0 , 1 ,  , n = b of structures , where each i   , then concatenate the direct pathways between successive i to i+1 , as determined by the morgan 
figure 7.morganhiggs algorithm ( 6 ) to construct an indirect pathway from secondary structure a to b. in line 2 of this algorithm , we use stochastic sampling of ding and lawrence ( 8) , as implemented in rnasubopt - p , and applied a modified version of dijkstra 's single source shortest path algorithm to determine low energy structures , whose greedy direct paths can be glued together for a pathway from a to b. morgan  higgs algorithm ( 6 ) to construct an indirect pathway from secondary structure a to b. in line 2 of this algorithm , we use stochastic sampling of ding and lawrence ( 8) , as implemented in rnasubopt - p , and applied a modified version of dijkstra 's single source shortest path algorithm to determine low energy structures , whose greedy direct paths can be glued together for a pathway from a to b. perhaps the simplest possible algorithm to find a nearly optimal ( direct ) pathway between a and b is to apply a greedy approach , where at each step we choose to remove a base pair belonging to a but not b , or add a base pair belonging to b but not a , where the choice of base pair to be removed or added is made so as to ensure the lowest energy next structure . pseudocode for this method , described by voss et al . ( 19 ) , is depicted in figure 8 . 

figure 8.greedy method to determine direct pathway between a and b , as described by voss et al .
secondary structures a , b can be considered to be sets of base pairs , so the requirement that   a  b means that every base pair of  belongs to either a or b. this condition ensures that the pathway produced is direct .
the notation dbp( ,  ) = 1 means that the base pair distance ( 30 ) between  ,  is 1 ; i.e.  ,  differ by one base pair .
moreover , since dbp( , b ) = dbp( , b ) + 1 , each iteration in the while loop ensures advancement by one base pair to the target structure b. it follows that the while loop involves dbp(a , b ) iterations .
greedy method to determine direct pathway between a and b , as described by voss et al .
secondary structures a , b can be considered to be sets of base pairs , so the requirement that   a  b means that every base pair of  belongs to either a or b. this condition ensures that the pathway produced is direct .
the notation dbp( ,  ) = 1 means that the base pair distance ( 30 ) between  ,  is 1 ; i.e.  ,  differ by one base pair .
moreover , since dbp( , b ) = dbp( , b ) + 1 , each iteration in the while loop ensures advancement by one base pair to the target structure b. it follows that the while loop involves dbp(a , b ) iterations .
however , we can benefit from a more randomized approach in which we randomly add or remove a valid base pair that yields a structure that is among the k lowest energy structures
since the result is clearly dependent on a parameter k we can iterate the same approach for several values of k and return the route with the lowest energy barrier . 

figure 9.semi-greedy method to determine direct pathway between a and b. the only difference between the greedy and semi - greedy method is that the latter randomly selects one of the k lowest energy neighbors ( step 5 ) rather than the minimum energy neighbor .
benchmarking indicates that the semi - greedy method generally outperforms the greedy method when determining low energy pathways between conformers of a riboswitch .
semi - greedy method to determine direct pathway between a and b. the only difference between the greedy and semi - greedy method is that the latter randomly selects one of the k lowest energy neighbors ( step 5 ) rather than the minimum energy neighbor .
benchmarking indicates that the semi - greedy method generally outperforms the greedy method when determining low energy pathways between conformers of a riboswitch .
indirect routes present more opportunities and challenges , since the space of possible routes increases considerably .
also , a purely greedy approach is not possible since the algorithm would not be able to escape from cycles . indeed , suppose that the structure a is the minimum free energy structure for the given rna sequence ; then the first step would add or remove a base pair , yielding a structure that is no longer the minimum free energy structure . in the next step , the added ( respectively removed ) base pair would then be removed ( respectively added ) , in order to return to the minimum free energy structure . for that reason
tabu search ( 20 ) is a well - studied combinatorial optimization method that entails a greedy strategy where a list of recently taken moves is placed temporarily on a 
tabu list , and can not be applied until removed from the tabu list . in figure 10 , we present pseudocode for a tabu semi - greedy algorithm ( rnatabupath ) to find nearly optimal , possibly indirect pathways between designated secondary structures a and b. the algorithm starts with the initial structure . at each successive step in the execution of the algorithm ,
we choose to add or remove that base pair resulting in the lowest energy ( greedy ) , after which the base pair is placed in the tabu list , hence can not be added or removed for a certain number of steps .
figure 10.tabu semi - greedy algorithm to compute near - optimal folding pathway between two designated structures a , b for a given rna sequence . in line 11 , we assume that  is obtained from  without using a base pair in the tabu list . the tabu list contains base pairs that were recently added or removed from an intermediate structure . when added to the tabu list , a base pair
fitness f of a structure is defined by f = e + w  bp , where e is energy of current structures , w is weight defined in line 2 and bp is  incremental  distance toward the target , i.e. 1 .
tabu semi - greedy algorithm to compute near - optimal folding pathway between two designated structures a , b for a given rna sequence . in line 11 , we assume that  is obtained from  without using a base pair in the tabu list .
the tabu list contains base pairs that were recently added or removed from an intermediate structure . when added to the tabu list , a base pair
fitness f of a structure is defined by f = e + w  bp , where e is energy of current structures , w is weight defined in line 2 and bp is  incremental  distance toward the target , i.e. 1 . as in every optimization algorithm we need to define the fitness function , f. the fitness function is a measure of quality of each state . in the case at hand ,
a state is a secondary structure  , and the fitness function must account for the free energy e( ) as well as the distance dbp( , b ) from  to the target structure b. hence , the fitness f( ) of secondary structure  is defined by 

 where w represents a weight that regulates the importance of reaching the target structure versus choosing a low energy structure . a low weight has the potential of driving the algorithm to structures that are too far away from the target , b , while a higher weight can quickly converge to the target structure at the expense of including higher energy intermediate structures in the path produced .
an intermediate value for weight w will tend to cause the algorithm to behave in a manner similar to that of the greedy algorithm for direct pathways . in order to avoid the latter , we have developed a  weight oscillation strategy that can be explained in the following steps : start with a given initial weight w0.increase the value of w when the distance to the target has not been improved for a number of iterations and restart from the structure found to be closest to the target.decrease the weight when the distance to the target is improved.if the weight reaches a certain value wmax , increase the value of w0 and restart the search ( with w = w0 ) . 
 
the distance to the target has not been improved for a number of iterations and restart from the structure found to be closest to the target .
if the weight reaches a certain value wmax , increase the value of w0 and restart the search ( with w = w0 ) .
our tabu strategy starts with the initial structure a , and in each step either adds or removes the base pair that minimizes the fitness function f. the base pair that has just been added or removed will be kept in a tabu list for a certain number of steps during which time it can not be added or removed to any structure in the pathway being constructed .
the fitness function f is adaptive , since it depends on the weight oscillation scheme .
additionally , the algorithm introduces an aspiration criterion for which a base pair can be changed ( even if it is tabu ) when the resulting structure reduces the best distance to the target found so far , provided that its free energy does not exceed that of the maximum energy of a structure in the pathway constructed so far .
additionally , we introduce two stochastic aspects to the tabu algorithm : the time a base pair remains on the tabu list , and the way to break ties when choosing the best base pair
consequently , we can start with a given value and iterate the algorithm using different values while maintaining the best pathway so far found .
note that we assume that in line 10 of the pseudocode of figure 10 , we assume that  is obtained from  without using a base pair in the tabu list unless the aspiration criterion just mentioned has been applied , and that the tabu list is updated .
, we found that by adding a semi - greedy component to tabu search , the resulting algorithm was substantially improved .
( 19 ) , described in figure 8 , is improved by adding a semi - greedy component for the search .
clearly , one could apply monte carlo and simulated annealing strategies to sample low energy folding pathways , as well as envision a genetic algorithm , that permits the crossover between folding pathways having a common source a and target b. nevertheless , the tabu semi - global approach of rnatabupath appears to be a very fast method to quickly determine near - optimal folding pathways .
the web site for rnapathfinder includes additional tools to determine the frequency of occurrence of secondary startures in ( say ) 1000 low energy folding pathways , and to determine the similarity between two pathways . in this section ,
we survey several known heuristics for determining folding pathways between two secondary structures , as well as present our novel semi - greedy and rnatabupath methods .
to explain the morgan  higgs greedy direct pathway algorithm , we first define the notion of a base pair 
clashing with another base pair  base pair ( i , j ) is said to clash with base pair ( x , y ) if either x  i  y  j or i  x  j  y. more generally , a base pair ( i , j ) clashes with a secondary structure a if there exists ( x , y )  a such that ( i , j ) clashes with ( x , y ) . the set of base pairs ( x , y )  a such that ( i , j ) clashes with ( x , y ) is denoted clash(i , j , a ) ; i.e. 

 with this definition , the morgan  higgs greedy algorithm repeatedly performs the following steps : ( i ) determine the base pair ( i , j ) belonging to b but not a which has minimum size clash set c , ( ii ) remove base pairs from c , and ( iii ) add base pairs in b that do not induce any new clashes .
figure 6.morganhiggs greedy algorithm ( 6 ) to construct a greedy direct pathway from secondary structure a to b. morgan 
higgs greedy algorithm ( 6 ) to construct a greedy direct pathway from secondary structure a to b. the morgan 
higgs algorithm to compute a nearly optimal ( possibly ) indirect pathway between secondary structures a , b proceeds as follows . by sampling ,
if either a or b does not belong to  , then add the missing structure to . define a complete , weighted , undirected graph g = ( v , e ) , where the set v of vertices consists of all structures in  , and the edge weight between any two structures  ,  is defined to be energy barrier max{e(i )  e( ) : 1  i  n } , where  = 0 ,  , n =  is the greedy direct pathway from  to  , as determined by the morgan 
morgan and higgs then apply  single link cluster ( slc ) algorithm , as described in ( 22 ) , in order to determine an optimal pathway , starting from structure a , proceeding by hopping from one low energy structure in  to another via a greedy direct pathway , and terminating by the structure b. in our implementation of the morgan 
higgs indirect algorithm , we sample low energy structures with respect to the turner energy model by applying the ding 
lawrence algorithm ( 8) , as implemented in rnasubopt - p from the vienna rna package ( step 2 of figure 7 ) . in place of the slc algorithm , we apply a modified form of dijkstra 's single source shortest path algorithm ( 23 ) , in order to determine a sequence a=0 , 1 ,  , n = b of structures , where each i   , then concatenate the direct pathways between successive i to i+1 , as determined by the morgan 
figure 7.morganhiggs algorithm ( 6 ) to construct an indirect pathway from secondary structure a to b. in line 2 of this algorithm , we use stochastic sampling of ding and lawrence ( 8) , as implemented in rnasubopt - p , and applied a modified version of dijkstra 's single source shortest path algorithm to determine low energy structures , whose greedy direct paths can be glued together for a pathway from a to b. morgan 
higgs algorithm ( 6 ) to construct an indirect pathway from secondary structure a to b. in line 2 of this algorithm , we use stochastic sampling of ding and lawrence ( 8) , as implemented in rnasubopt - p , and applied a modified version of dijkstra 's single source shortest path algorithm to determine low energy structures , whose greedy direct paths can be glued together for a pathway from a to b.
perhaps the simplest possible algorithm to find a nearly optimal ( direct ) pathway between a and b is to apply a greedy approach , where at each step we choose to remove a base pair belonging to a but not b , or add a base pair belonging to b but not a , where the choice of base pair to be removed or added is made so as to ensure the lowest energy next structure .
figure 8.greedy method to determine direct pathway between a and b , as described by voss et al .
secondary structures a , b can be considered to be sets of base pairs , so the requirement that   a  b means that every base pair of  belongs to either a or b. this condition ensures that the pathway produced is direct .
the notation dbp( ,  ) = 1 means that the base pair distance ( 30 ) between  ,  is 1 ; i.e.  ,  differ by one base pair .
moreover , since dbp( , b ) = dbp( , b ) + 1 , each iteration in the while loop ensures advancement by one base pair to the target structure b. it follows that the while loop involves dbp(a , b ) iterations .
greedy method to determine direct pathway between a and b , as described by voss et al .
secondary structures a , b can be considered to be sets of base pairs , so the requirement that   a  b means that every base pair of  belongs to either a or b. this condition ensures that the pathway produced is direct .
the notation dbp( ,  ) = 1 means that the base pair distance ( 30 ) between  ,  is 1 ; i.e.  ,  differ by one base pair .
moreover , since dbp( , b ) = dbp( , b ) + 1 , each iteration in the while loop ensures advancement by one base pair to the target structure b. it follows that the while loop involves dbp(a , b ) iterations .
however , we can benefit from a more randomized approach in which we randomly add or remove a valid base pair that yields a structure that is among the k lowest energy structures .
since the result is clearly dependent on a parameter k we can iterate the same approach for several values of k and return the route with the lowest energy barrier . 

figure 9.semi-greedy method to determine direct pathway between a and b. the only difference between the greedy and semi - greedy method is that the latter randomly selects one of the k lowest energy neighbors ( step 5 ) rather than the minimum energy neighbor .
benchmarking indicates that the semi - greedy method generally outperforms the greedy method when determining low energy pathways between conformers of a riboswitch .
semi - greedy method to determine direct pathway between a and b. the only difference between the greedy and semi - greedy method is that the latter randomly selects one of the k lowest energy neighbors ( step 5 ) rather than the minimum energy neighbor .
benchmarking indicates that the semi - greedy method generally outperforms the greedy method when determining low energy pathways between conformers of a riboswitch .
indirect routes present more opportunities and challenges , since the space of possible routes increases considerably .
also , a purely greedy approach is not possible since the algorithm would not be able to escape from cycles . indeed , suppose that the structure a is the minimum free energy structure for the given rna sequence ; then the first step would add or remove a base pair , yielding a structure that is no longer the minimum free energy structure . in the next step , the added ( respectively removed ) base pair
would then be removed ( respectively added ) , in order to return to the minimum free energy structure .
for that reason , it makes sense to exclude certain moves at certain times during the search .
tabu search ( 20 ) is a well - studied combinatorial optimization method that entails a greedy strategy where a list of recently taken moves is placed temporarily on a 
tabu list , and can not be applied until removed from the tabu list . in figure 10 , we present pseudocode for a tabu semi - greedy algorithm ( rnatabupath ) to find nearly optimal , possibly indirect pathways between designated secondary structures a and b. the algorithm starts with the initial structure . at each successive step in the execution of the algorithm , we choose to add or remove that base pair resulting in the lowest energy ( greedy ) , after which the base pair is placed in the tabu list , hence can not be added or removed for a certain number of steps .
figure 10.tabu semi - greedy algorithm to compute near - optimal folding pathway between two designated structures a , b for a given rna sequence . in line 11 , we assume that  is obtained from  without using a base pair in the tabu list . the tabu list contains base pairs that were recently added or removed from an intermediate structure . when added to the tabu list , a base pair is given a time stamp .
fitness f of a structure is defined by f = e + w  bp , where e is energy of current structures , w is weight defined in line 2 and bp is  incremental  distance toward the target , i.e. 1 .
tabu semi - greedy algorithm to compute near - optimal folding pathway between two designated structures a , b for a given rna sequence . in line 11 , we assume that  is obtained from  without using a base pair in the tabu list .
the tabu list contains base pairs that were recently added or removed from an intermediate structure . when added to the tabu list , a base pair
fitness f of a structure is defined by f = e + w  bp , where e is energy of current structures , w is weight defined in line 2 and bp is  incremental  distance toward the target , i.e. 1 . as in every optimization algorithm we need to define the fitness function , f. the fitness function is a measure of quality of each state . in the case at hand ,
a state is a secondary structure  , and the fitness function must account for the free energy e( ) as well as the distance dbp( , b ) from  to the target structure b. hence , the fitness f( ) of secondary structure  is defined by 

 where w represents a weight that regulates the importance of reaching the target structure versus choosing a low energy structure .
a low weight has the potential of driving the algorithm to structures that are too far away from the target , b , while a higher weight can quickly converge to the target structure at the expense of including higher energy intermediate structures in the path produced .
an intermediate value for weight w will tend to cause the algorithm to behave in a manner similar to that of the greedy algorithm for direct pathways . in order to avoid the latter , we have developed a  weight oscillation strategy that can be explained in the following steps : start with a given initial weight w0.increase the value of w when the distance to the target has not been improved for a number of iterations and restart from the structure found to be closest to the target.decrease the weight when the distance to the target is improved.if the weight reaches a certain value wmax , increase the value of w0 and restart the search ( with w = w0 ) . 
  start with a given initial weight w0 .
increase the value of w when the distance to the target has not been improved for a number of iterations and restart from the structure found to be closest to the target . decrease the weight when
if the weight reaches a certain value wmax , increase the value of w0 and restart the search ( with w = w0 ) .
our tabu strategy starts with the initial structure a , and in each step either adds or removes the base pair that minimizes the fitness function f. the base pair that has just been added or removed will be kept in a tabu list for a certain number of steps during which time it can not be added or removed to any structure in the pathway being constructed .
the fitness function f is adaptive , since it depends on the weight oscillation scheme .
additionally , the algorithm introduces an aspiration criterion for which a base pair can be changed ( even if it is tabu ) when the resulting structure reduces the best distance to the target found so far , provided that its free energy does not exceed that of the maximum energy of a structure in the pathway constructed so far .
additionally , we introduce two stochastic aspects to the tabu algorithm : the time a base pair remains on the tabu list , and the way to break ties when choosing the best base pair .
consequently , we can start with a given value and iterate the algorithm using different values while maintaining the best pathway so far found .
note that we assume that in line 10 of the pseudocode of figure 10 , we assume that  is obtained from  without using a base pair in the tabu list unless the aspiration criterion just mentioned has been applied , and that the tabu list is updated .
, we found that by adding a semi - greedy component to tabu search , the resulting algorithm was substantially improved .
( 19 ) , described in figure 8 , is improved by adding a semi - greedy component for the search .
clearly , one could apply monte carlo and simulated annealing strategies to sample low energy folding pathways , as well as envision a genetic algorithm , that permits the crossover between folding pathways having a common source a and target b. nevertheless , the tabu semi - global approach of rnatabupath appears to be a very fast method to quickly determine near - optimal folding pathways .
the web site for rnapathfinder includes additional tools to determine the frequency of occurrence of secondary startures in ( say ) 1000 low energy folding pathways , and to determine the similarity between two pathways . in this section ,
we survey several known heuristics for determining folding pathways between two secondary structures , as well as present our novel semi - greedy and rnatabupath methods .
in this section , we present summary results on folding pathways and energy barriers computed for each of the algorithms : greedy ( 19 ) , semi - greedy , rnatabupath , morgan  higgs direct ( 6 ) , morgan  higgs indirect ( 6 ) , findpath ( 18 ) and barriers ( 2,15 ) . due to the stochastic nature of the semi - greedy method and rnatabupath
, we report the best results found over 1000 runs . in the case of rnatabupath ,
fitness of the current structure is defined by f = e + w  bp , where e is the free energy of the current structure , and bp is the 
rnatabupath allows the user to input parameters wmin , wmax that confine the weight w  [ wmin , wmax ] .
reported values were for wmin = 1 and wmax = 7 , which are the default values on the web server .
table 1 presents the  energy barrier in the pathway a=0 , 1 ,  , n = b between low energy structures a and b of known conformational switches , where energy barrier is defined to be max{e(i )  e(a ) : i = 1 ,  , n}. structures a and b are two metastable states of five riboswitches , guanine riboswitch from bacillus subtilis ( rb1 ) , adenine riboswitch from v. vulnificus ( rb2 ) , s - adenylmethionine riboswitch from thermoanaerobacter tecongensis ( rb3 ) , thymine pyrophosphate riboswitch from t. tecongensis ( rb4 ) , and xpt - pbux riboswitch from b. subtilis ( rb5 ) , whose metastable secondary structures were found by in - line probing experiments of various groups .
table 1 also contains results for some conformational switches found on the parnass web site , http://bibiserv.techfak.uni-bielefeld.de/parnass/examples.html ; however , since the metastable structures for the latter conformational switches have not been experimentally determined , we ran the software rnabor , which determines for each integer value of  , the minimum free energy structure mfe( ) and partition function z( ) over all -neighbors of the minimum free energy structure . here a structure  is said to be a -neighbor of structure  if base pair distance between  ,  is . [ see ( 25 ) for details on rnabor . ] for the conformational switches taken from the parnass web site , we defined a to be the minimum free energy structure and b to be that structure which is the minimum free energy structure over all -neighbors of a , where 10 <  and the output of rnabor indicated a second peak at the value . 
 table 1.algorithm benchmarks for computing folding pathways between two low energy secondary structuresinstancegreedysemi - greedyrnatabupathgreedymhindirectmhfindpathbarriersrb132.8025.2424.0426.2423.9924.04rb214.649.207.2510.0010.008.20*rb324.8022.7017.9028.4020.0022.40rb416.9016.9016.9016.9016.9016.90rb533.3025.6724.5426.7426.7424.54hok36.3733.7029.6636.3036.3028.5sl14.0914.0912.9018.2016.2013.0011.80attenuator11.509.008.6012.6014.708.708.30s157.107.106.609.709.707.106.60s - box leader7.105.305.2010.209.305.20*thim leader21.4416.6714.8420.5731.0016.13ms28.306.606.6011.7011.706.60*hdv25.5021.7017.0023.5322.5017.4dsra8.308.308.2014.6010.778.308.00ribd leader13.8411.709.5018.1116.9010.71amv10.006.405.8015.610.45.80*alpha operon6.506.506.509.906.506.50*hiv-1 leader14.2813.4911.3017.9018.509.30greedy refers to our implementation of voss et al . ( 19 ) , where a direct path is constructed by choosing the lowest energy base pair to remove or add at each step ;  semi - greedy refers to to our modification of voss et al . ( 19 ) , where a direct path is constructed by choosing one of the k lowest energy base pairs to remove or add at each step ; rnatabupath refers to our semi - greedy tabu search method described in the text ; 
indirectmh refers to our implementation of morgan  higgs method ( 6 ) to produce a possibly indirect path ; 
18 ) with look - ahead parameter k = 10 ; barriers refers to the exact method of flamm et al .
in each case , near - optimal low energy pathways between two low energy secondary structures of five different riboswitches : guanine riboswitch from b. subtilis ( rb1 ) , adenine riboswitch from v. vulnificus ( rb2 ) , s - adenylmethionine riboswitch from t. tecongensis ( rb3 ) , thymine pyrophosphate riboswitch from t. tecongensis ( rb4 ) and xpt - pbux riboswitch from b. subtilis ( rb5 ) .
secondary structures for rb1rb5 were experimentally determined ; see wakeman et al . ( 16 ) for rb1rb4 and mandal et al .
sequences of additional conformational switches were taken from the parnass web site http://bibiserv.techfak.uni-bielefeld.de/parnass/examples.html , courtesy of the giegerich lab .
for the latter , the two low energy structures were taken to be the minimum free energy structure a and the structure b determined by rnabor ( 25 ) to be the minimum free energy structure over all structures having base pair distance k with a , where 10  k and a second peak was found at position k in the output of rnabor .
energy barrier in the pathway a = 0 , 1 ,  , n = b from a to b is here defined to be max{e(i )  e(a ) } : i = 1 ,  , n } , where free energy is measured in kcal / mol , as computed by rnaeval from vienna rna package .
notation used in last column given as follows : means barriers could not converge ; * means that either structure a or b is not locally optimal , hence barriers could not be directly applied .
however , one could apply barriers in the following manner , as suggested by an anonymous referee .
given non - locally optimal structures a , b , one can first determine locally optimal structures  ,  that , respectively , contain a , b , then apply barriers to find an optimal path between  , . this will yield a near - optimal path from a to b. boldface numbers indicate the minimal barrier energy found by one of the heuristic algorithms , while underlined numbers indicate the minimum barrier energy found by the exact algorithm , barriers .
( 19 ) , where a direct path is constructed by choosing the lowest energy base pair to remove or add at each step ;  semi - greedy refers to to our modification of voss et al .
( 19 ) , where a direct path is constructed by choosing one of the k lowest energy base pairs to remove or add at each step ; rnatabupath refers to our semi - greedy tabu search method described in the text ; 
indirectmh refers to our implementation of morgan  higgs method ( 6 ) to produce a possibly indirect path ; 
18 ) with look - ahead parameter k = 10 ; barriers refers to the exact method of flamm et al .
near - optimal low energy pathways between two low energy secondary structures of five different riboswitches : guanine riboswitch from b. subtilis ( rb1 ) , adenine riboswitch from v. vulnificus ( rb2 ) , s - adenylmethionine riboswitch from t. tecongensis ( rb3 ) , thymine pyrophosphate riboswitch from t. tecongensis ( rb4 ) and xpt - pbux riboswitch from b. subtilis ( rb5 ) .
sequences of additional conformational switches were taken from the parnass web site http://bibiserv.techfak.uni-bielefeld.de/parnass/examples.html , courtesy of the giegerich lab .
for the latter , the two low energy structures were taken to be the minimum free energy structure a and the structure b determined by rnabor ( 25 ) to be the minimum free energy structure over all structures having base pair distance k with a , where 10  k and a second peak was found at position k in the output of rnabor .
energy barrier in the pathway a = 0 , 1 ,  , n = b from a to b is here defined to be max{e(i )  e(a ) } : i = 1 ,  , n } , where free energy is measured in kcal / mol , as computed by rnaeval from vienna rna package .
notation used in last column given as follows : means barriers could not converge ; * means that either structure a or b is not locally optimal , hence barriers could not be directly applied .
however , one could apply barriers in the following manner , as suggested by an anonymous referee . given
non - locally optimal structures a , b , one can first determine locally optimal structures  ,  that , respectively , contain a , b , then apply barriers to find an optimal path between  , . this will yield a near - optimal path from a to b. boldface numbers indicate the minimal barrier energy found by one of the heuristic algorithms , while underlined numbers indicate the minimum barrier energy found by the exact algorithm , barriers .
for technical reasons having to do with computation of the partition function , the treatment of dangles in rnabor is identical to that of vienna rna package rnafold with option d2 . in some instances ,
the metastable structure we chose using rnabor was no longer locally optimal under the d1 treatment of dangle , which latter is used in all the algorithms appearing in table 1 . in particular , we should mention that one must explicitly use d1 option with rnasubopt , to ensure that rnafold , rnaeval and barriers all use the same treatment of dangles .
due to the energy model differences ( d2 versus d1 ) in using rnabor to choose one of the metastable structures , barriers could not be used in some instances  rb2 , s - box leader , ms2 , amv and alpha operon .
we see that the greedy approach is simple , but yields considerably poorer results than other methods tested . however , a small change such as a semi - greedy component yields great improvements .
tabu search for indirect routes outperforms both greedy and semi - greedy approaches ( data for the tabu greedy method is not shown ) . in the semi - greedy algorithm and rnatabupath , we experimented with different choices of the value k , where randomly one of the best k neighbors is chosen .
after computational experiments over the range of lengths typical for conformational switches , we fixed the value k = 8 for semi - greedy algorithm and k = 5 for rnatabupath .
the initial weight w0 in rnatabupath ranges from 1 to 7 , the default setting for the web server , although best results for this range depend on the input sequence . in general , w0  [ 4 , 7 ] works better for larger sequences .
higgs direct and indirect algorithms did not perform well in all but one instance ; morgan  higgs indirect algorithm curiously outperformed all algorithms for rb1 .
in general , findpath is a very fast algorithm that produces excellent quality direct pathways , with barrier energies often equal or close to those of rnatabupath . in the case of hok - rna and hiv-1 leader
, findpath outperformed all other approaches . in our benchmarking , we set the look - ahead of findpath to be 10 ; often increasing the look - ahead to 100 did not change the results . however , in the 396 nt hok - rna , findpath improved dramatically with increased look - ahead k : barrier energy of 28.5 for k = 10 , 28.17 for k = 20 , 23.5 for k = 100 , 22.7 for k = 200 , 21.4 for k = 500 and k = 1000 .
figures 2 and 3 demonstrate cases where a well - chosen indirect pathway necessarily has lower barrier energy than that of any direct pathway . applying this principle to 304 examples derived from pseudoknotted structures in pseudobase ( 21 )
, we found that in roughly half the examples , rnatabupath and findpath produced the same barrier energy , while in all other instances , rnatabupath produced a lower barrier energy barrier than did findpath ; indeed , the maximum difference in barrier energy was 6.51 , while the average was 1.93 kcal / mol with standard deviation of 1.45 .
figure 11 depicts a folding pathway computed by rnatabupath between the two meta - stable secondary structures of the adenine riboswitch from v. vulnificus ( rb2 ) ( 16 ) .
figure 12 depicts the free energy of intermediate structures in this pathway as a function of step number . 

figure 11.comparison of best direct pathway ( above ) and best indirect pathway ( below ) , as found by rnatabupath between the two metastable secondary structures of the adenine riboswitch from v. vulnificus ( rb2 ) , as reviewed in wakeman et al .
figure 12.graph of free energy of intermediate structure as a function of step number or index in the rnatabupath folding pathway between two metastable secondary structures of the adenine riboswitch from v. vulnificus ( rb2 ) , corresponding to data from figure 11 .
dotted lines depict a similar graph for a folding pathway computed by the semi - global method .
comparison of best direct pathway ( above ) and best indirect pathway ( below ) , as found by rnatabupath between the two metastable secondary structures of the adenine riboswitch from v. vulnificus ( rb2 ) , as reviewed in wakeman et al .
graph of free energy of intermediate structure as a function of step number or index in the rnatabupath folding pathway between two metastable secondary structures of the adenine riboswitch from v. vulnificus ( rb2 ) , corresponding to data from figure 11 .
dotted lines depict a similar graph for a folding pathway computed by the semi - global method . one useful application of rnatabupath is to provide an energy upper bound for subsequent application of barriers , an observation pointed out by an anonymous referee . specifically , given rna sequence s and two metastable structures a , b , let e0 denote the minimum free energy of s and let e(a ) denote the free energy of structure a. if e is the barrier energy computed by rnatabupath ( or another method ) for a folding pathway from a to b , then barriers with bound e + ( e(a )  e0 ) will compute an optimal pathway , provided it converges . the barrier energies obtained by barriers in table 1
since barriers is the only exact algorithm , when it converges , a provably optimal pathway is produced . in the cases of rb1 , rb3 ,
rb4 , rb5 , thim leader , ribd leader and hiv-1 leader , barriers did not converge , even when started with the energy bound obtained by rnatabupath .
molecular folding pathways are low energy routes taken along an energy surface . as previously noticed by morgan and higgs ( 6 ) ,
this is clear from the toy example presented in figures 2 and 3 . in other data ,
we see how the creation of a base pair in a region with no base pairs leads to the stabilization of other secondary structures along the folding pathway .
since barriers is an exact algorithm , it should be used whenever possible ; i.e. one should first apply findpath or rnatabupath to obtain an energy upper bound for subsequent application of barriers . in other cases , findpath and rnatabupath
if large type - h pseudoknots appear in the structure obtained by adjoining two metastable structures , then rnatabupath is likely to be the best algorithm , since indirect pathways will have lower barrier energy in this case . to assist those interested in computing near - optimal folding pathways , we have created the web server rnapathfinder , located at http://bioinformatics.bc.edu/clotelab/rnapathfinder .
in addition to supporting rnatabupath computations , source code can be downloaded for several algorithms discussed in this article .
fundacion caja madrid ( to i.d . ) ; national science foundation ( grants dbi-0543506 and dms-0817971 to p.c . and w.a.l . ) ; rna ontology consortium ( to i.d .
; deutscher akademischer austauschdienst ( to p.c . for funding a visit to martin vingron 's ; group in the max planck institute of molecular genetics ) ; digiteo foundation ( to p.c . ) . funding for open access charge : national science foundation ( grant dbi-0543506 ) .
 given an rna sequence and two designated secondary structures a , b , we describe a new algorithm that computes a nearly optimal folding pathway from a to b. the algorithm , rnatabupath , employs a tabu semi - greedy heuristic , known to be an effective search strategy in combinatorial optimization . folding pathways , sometimes called routes or trajectories , 
 are computed by rnatabupath in a fraction of the time required by the barriers program of vienna rna package . 
 we benchmark rnatabupath with other algorithms to compute low energy folding pathways between experimentally known structures of several conformational switches . 
 the rnapathfinder web server , source code for algorithms to compute and analyze pathways and supplementary data are available at http://bioinformatics.bc.edu/clotelab/rnapathfinder . 