BACKGROUND
in  <cit> , borenstein and ruppin define neutrality of an rna sequence a=a
<dig> …,a
n by η=1−〈d〉n, where in this section 〈d〉 denotes the average, taken over all 3n single-point mutants of a, of the base pair distance d
bp between the minimum free energy  structure s
<dig> of a and the mfe structures of single-point mutants of a. an rna sequence a is then defined to be robust if η is greater than the average neutrality of  <dig> control sequences generated by the program rnainverse  <cit> , which fold into the same target structure s
<dig>  the main finding of  <cit>  is that precursor micrornas  exhibit a significantly higher level of mutational robustness than random rna sequences having the same structure.
to control for sequence composition bias in their computational study, the authors selected sequences from the output of rnainverse, whose dinucleotide composition was similar to that of wild type pre-mirna .
since the filtering step required enormous run time and computational resources, the authors restricted their attention to a small set of  <dig> micrornas, generating only  <dig> control sequences per microrna.
borenstein and ruppin conclude that robustness of precursor micrornas is not the byproduct of a base composition bias or of thermodynamic stability.
subsequently rodrigo et al.
<cit>  undertook a similar analysis for bacterial small rnas, also using the program rnainverse, albeit using somewhat different definitions – precise definitions are given in “formal definitions of robustness” section.
the main finding of  <cit>  was that bacterial sncrnas are not significantly robust when compared with  <dig> sequences having the same structure, as computed by rnainverse; however, bacterial sncrnas tend to be significantly plastic, in the sense that the ensemble of low energy structures is structurally diverse.
unlike the case of precursor micrornas  <cit> , rodrigo et al.
did not control for sequence compositional bias.
this raises the question of whether the control sequences analyzed in  <cit>  are representative or to what extent features shared by sequences output by the program rnainverse are artifacts of the program used.
indeed, the number of rna sequences that fold into a given target structure can be astronomically large.
over a few weeks, before we elected to terminate the execution, our state-of-the-art inverse folding software rnaifold  <cit>  generated  <dig> , <dig> many 52-nt sequences that fold exactly into the mfe secondary structure s
<dig> of hiv- <dig> ribosomal frameshift stimulating signal from the gag-pol overlap region af <dig> /1631- <dig>  and which additionally code 17-mer peptides in the gag and pol reading frames having amino acids that appear in gag/pol peptides found in the los alamos hiv- <dig> database  <cit> .
the number of  <dig> nt rna sequences that fold into target s
<dig> without additionally imposing the constraint of coding particular peptides in overlapping gag/pol reading frames is certain to dwarf the previous number.
moreover, the number of sequences that fold into the mfe structure of an animal precursor microrna  or into the mfe structure of bacterial sncrna  is certain to be even more daunting.
different inverse folding algorithms have adopted different strategies to generate sequences that fold into a user-specified target secondary structure s
<dig>  for instance, rnainverse  <cit>  performs an adaptive walk, in one step of which a nucleotide in the current sequence is mutated and subsequently accepted if the base pair distance between the minimum free energy  structure of the mutated sequence and the target structure s
<dig> is reduced.
nupack design  <cit>  selects a candidate mutation position with probability proportional to its contribution to the ensemble defect ∈s
<dig> and s
0=i indicates i is unpaired in s
0).
rnaifold cp-design  <cit>  uses constraint programming to systematically explore the search tree of all inverse folding solutions in an order determined by certain heuristics.
accordingly, one cannot claim that the collection of sequences generated by any particular inverse folding algorithm is representative of the astronomically large space of all inverse folding solutions – indeed, each inverse folding algorithm has an inherent but unknown bias.
in this paper, we describe the algorithm rnadualpf, which generates sequences which have low free energy with respect to a user-specified target structure s
<dig> – i.e.
the inherent bias of rnadualpf is known, unlike the situation of other inverse folding algorithms.
we show that rnadualpf is extremely fast software for generating sequences that approximately fold into s
0; moreover, in a postprocessing step, one can filter the output of rnadualpf to select sequences that exactly fold into s
<dig>  rnadualpf additionally allows the user to specify iupac codes to constrain certain nucleotide positions as well as to control the gc-content of all generated sequences.
sampling is performed in a manner distinct but somewhat analogous to that by which sfold  <cit>  and rnasubopt -p  <cit>  sample representative secondary structures from the boltzmann ensemble of all structures of a given sequence.
using rnadualpf, we perform a pilot study that is similar, though not identical, to that of  <cit>  for two classes of rna:  <dig> c. elegans precursor microrna from mirbase  <cit>  and the bacterial small noncoding rnas previously analyzed in  <cit> .
finally, it should be noted that, although rnadualpf was developed entirely independently of the work of reinharz et al.
<cit> , one can view our c-program as an extension of python program incarnation  <cit>  to the full turner energy model, where additionally gc-content is rigorously handled.
this point will be discussed further in the conclusion.
formal definitions of robustness
let a=a
<dig> …,a
n denote an arbitrary rna sequence, where ai∈n={a, u, g, c}, a secondary structure s of a is a set of base pairs  satisfying the following conditions:  if ∈s then a
i,a
j constitute a watson-crick or gu wobble pair, i.e.
ij∈ℬ which is the set {au,ua,gc,cg,gu,ug}.
if ∈s then i+θ<j, where θ= <dig> .
if ∈s and ∈s, then either i<k<ℓ<j or k<i<j<ℓ or i<j<k<ℓ or k<ℓ<i<j.
the collection of all secondary structures of the rna sequence a is denoted 𝕊𝕊, and the free energy  <cit>  of s is denoted by e, or simply by e provided that the sequence a is clear from context.
the boltzmann probability
p=p
a for structure s of a is defined by exp/rt)/z, where the partition function z=z=∑s∈𝕊𝕊exp/rt).
given two secondary structures s,t of a, the base pair distance
d
bp between s and t is defined to be the size of the symmetric difference of s,t, i.e.
|s−t|+|t−s|.
in  <cit> , rodrigo et al.
define intrinsic distance
<dig> d0=∑s,tp·p·dbp
i.e.
intrinsic distance is another name for ensemble diversity earlier defined in  <cit> , and computed by vienna rna package  <cit> .
plasticity is defined in  <cit>  to be normalized ensemble diversity; i.e.
<dig> p=d0n/ <dig>
obtained by dividing ensemble diversity by  the maximum possible number n/ <dig> of base pairs in a structure of a. given two rna sequences a=a
<dig> …,a
n and b=b
<dig> …,b
n of the same length n, rodrigo et al.
define d
<dig> to be the expected base pair distance between structures of a and structures of b minus the ensemble diversity of a, i.e.
<dig> d1=∑s∈𝕊𝕊∑t∈𝕊𝕊pa·pb·dbp−d <dig>
since d
<dig> is not symmetric, this measure is not a metric.
in contrast, ensemble distance as described in  <cit>  is a valid metric, defined by the following:
<dig> dv=∑s∈𝕊𝕊∑t∈𝕊𝕊pa·pb·dbp−d0+d02=∑i<j−pi,j) <dig>
in  <cit> , rodrigo et al.
define the mutational robustness
<dig> rm=1−〈d1〉n/ <dig>
where 〈d
1〉 denotes the average value of d
<dig> taken over all single point mutants a’ of a. since d
<dig> is not a true metric, we replace it by the metric d
v in our computation of mutational robustness.
clearly both notions are closely related.
implementation
in  <cit> , mccaskill described a cubic time algorithm to compute the partition function
<dig> z=z=∑s∈𝕊𝕊exp/rt)
for an rna sequence a=a
<dig> …,a
n, where the sum is taken over all secondary structures 𝕊𝕊 of a, e denotes the free energy for the structure s of a with respect to the turner energy parameters  <cit> , r denotes the universal gas constant and t is absolute temperature.
subsequently ding and lawrence  <cit>  described how to use the partition function together with a simple backtracking strategy to sample secondary structures of a from the boltzmann ensemble of low energy structures.
if s
<dig> is a given secondary structure of length n, we define the dual partition function
<dig> z∗=z∗=∑a∈𝔸𝔸exp/rt)
where the sum is taken over all rna sequences a=a
<dig> …,a
n of length n that are compatible with structure s
<dig>  i.e.
a
i,a
j constitute a watson-crick or wobble pair for each base pair ∈s
<dig>  the set of all rna sequences that are compatible with s
<dig> is denoted by 𝔸𝔸.
note that if a sequence a is not compatible with the target structure s
<dig>  then the energy e is infinite, so the corresponding boltzmann factor exp/rt) is zero and the sum in eq.
could have been written over all sequences of the same length as s
<dig>  here we describe the efficient software rnadualpf to compute the dual partition function
z
∗ and to sample from the low energy ensemble of sequences that are compatible with a given secondary structure s
<dig>
dual partition function
if s is a secondary structure on sequence a=a
<dig> …,a
n, then the length of s, denoted by ℓ, is equal to n, while the size of s, denoted by |s|, is the number of base pairs belonging to s. similarly, if secondary structure s is restricted to the interval , where 1≤i≤j≤n, then the length of the restriction of s to , denoted by ℓ, is equal to j−i+ <dig>  while the size of the restriction of s to , denoted by |s|, is the number of base pairs  of s that satisfy i≤x<y≤j.
given an rna sequence a=a
<dig> …,a
n, the mccaskill algorithm  <cit>  computes the partition function z defined in eq.
.
when a is clear from context, z is usually denoted by z.
given a target secondary structure s
<dig>  we describe below an algorithm to compute the dual partition function
z
∗, defined as the sum of all boltzmann factors exp), where the sum is taken over all rna sequences a∈𝔸𝔸.
unlike the mccaskill algorithm, which requires time that is cubic in the length of a, the algorithm presented below requires time that is  linear <dig> in the length of s
<dig>  our algorithm is motivated by the initialization step of the algorithm info-rna  <cit> , in which a sequence is determined, for which the free energy with respect to target structure s
<dig> is a minimum – i.e.
info-rna determines argmina
e.
the algorithm specification requires the notation z
∗, which denotes the sum
<dig> z∗=∑a,ai=x,aj=yexp−ea,s0/rt
of boltzmann factors for sequences a=a
i,…,a
j for which a
i=x,a
j=y, and for the restriction s
<dig>  defined by
<dig> s0=∈s0:i≤x<y≤j.
the function z
∗ is defined for all base pairs ∈s
0; these values will be stored in an array, whose rows index base pairs of s
<dig>  and whose columns are indexed by the six canonical base pairs gc, cg, au, ua, gu, ug .
once z
∗ has been computed for all base pairs that are visible, i.e.
for which there is no base pair  for which x<i<j<y, we can compute the full partition function z
∗.
au
cg
gc
ua
gu
ug
z
s
0
the first column indicates the base pair index which dictates the order in which the dual partition function is computed for different loops closed by the base pair , where we the index of base pair  is defined to be the rank of  in the total ordering defined in eq.
.
columns i and j indicate the opening and closing positions of each base pair.
type indicates the type of element in the secondary structure closed by each base pair, where r. bulge stands for right bulge, stack for stacking base pair, and int.
loop for interior loop.
the dual partition function
z
∗ of the substructure closed by base pair  appears in the rightmost column, while the partition function z
∗ for each of the six canonical base pairs is given in columns 5- <dig>  note that for base pair  <dig>  sequence constraints depicted in fig.
<dig> force i and j to be instantiated respectively to g and c, hence the dual partition function z
∗ is zero for any base pair different than gc.
the last column of the last row of the table shows the total dual partition function z
∗ for the target structure s
0
following  <cit> , we define a total ordering on base pairs  belonging to the target structure s
<dig> that satisfy the following precedence rule for any two base pairs ,.
<dig> ≺⇔x<i<j<yori<j<x<y
from this ordering, we assign a base pair index to each base pair , which is defined to be the rank of  in the total ordering.
the following definitions correspond to the turner nearest neighbor energy model  <cit> , which is an additive loop model where a loop closed by external base pair  is designated as a k-loop, if the loop contains k base pairs interior to .
therefore, hairpin loops are 0-loops; base pair stacks, bulge loops and internal loops are 1-loops; and multiloops are k-loops for k≥ <dig> -way junctions), where the additional count is due to the outer component adjacent to   <cit> .
since au-base pairs that close a loop are energetically unfavorable, in the turner energy model, there is an au-penalty we now define:
<dig> eau= <dig> ifis the outermost pair in a stem ofs <dig> having au,ua,gu,ug0otherwise.
this au-penalty is applied only if  is a base pair adjacent to a triloop, a bulge, an internal loop or a multiloop, or if it is the outermost base pair of an external loop in target structure s
<dig>  and  is instantiated by one of the pairs au, ua, gu, ug.
when base-paired positions i,j are clear from the context, we write e
au.
here, we assume that in parsing the input target structure, a list bpcloseelorml has been created of those base pairs , which close either an external loop or a multiloop.
let i be the indicator function, it follows that if  closes an external loop or multiloop, then exp−i·eaurt is the boltzmann factor for a special au-penalty, otherwise this factor equals  <dig>  for clarity in the notation, this factor is denoted by ert).
note that this term is distinct from the factor exprt) applied to base pairs adjacent to a triloop, a bulge or an internal loop, which does not depend on the indicator function.
hairpins
let  close a hairpin in s
<dig>  the hairpin free energy term h, arising solely from entropic considerations, is defined by
<dig> h=hairpineifj−i−1≤30hairpine+ <dig> rtlnj−i−130otherwise
where hairpine designates the hairpin free energy obtained from table look-up, when j−i−1≤ <dig>
triloop
let triloop
x,y denote the collection of special triloops, xabcy, having an energy bonus triloope.
<dig> z∗=e−eauirt·exp−h+eaurt×43−|triloopx,y|+∑abc∈triloopx,yexp−triloopert
tetraloop
let tetraloop
x,y denote the collection of special tetraloops, xabcdy, having an energy bonus tetraloope.
similarly, given nucleotides n <dig> n2∈n, tetraloop
x,y denotes the collection of special tetraloops of the form xn
1
abn
2
y. define z
∗ by
<dig> z∗=e−eauirt·exp−hrt×∑n <dig> n2∈nexp−mismatchrt×42−|tetraloopx,y|+∑ab∈tetraloopx,yexp−tetraloopert
hexaloop let hexaloop
x,y denote the collection of special hexaloops, xabcdefy, having an energy bonus hexaloope.
similarly, given nucleotides n
<dig> n
<dig>  hexaloop
x,y denotes the collection of special hexaloops of the form xn
1
abcdn
2
y. define z
∗ by
<dig> z∗=e−eauirt·exp−hrt×∑n <dig> n2∈nexp−mismatchx,y,n <dig> n2rt×44−|hexaloopx,yn <dig> n2|+∑ab∈hexaloopx,yn <dig> n2exp−hexaloopexn1abcdn2yrt
hairpin size exceeds four and is different than six
define z
∗ by
<dig> z∗=e−eauirt·exp−hrt×∑n <dig> n2∈nexp−mismatchx,y,n <dig> n2rt·4j−i− <dig>
stacked base pairs, bulges and internal loops
here, we consider the case of a 1-loop, which comprises the case of stacked base pairs, bulges and internal loops.
the following cases correspond to each possibility.
stacked base pair
in this case,  stacks on the base pair , and the partition function z
∗ has been computed.
let stack denote the free energy of base stack 5′−xu−3′3′−yv−5′ obtained by table look-up.
<dig> z∗=e−eauirt·∑uv∈ℬexp−stackrt×z∗i+ <dig> j− <dig> u,v
bulge loop
in this case,  closes a bulge in s
<dig>  since bulge size may exceed the values in table look-up, we define the free energy for a bulge of size r by
<dig> bulge=bulgeeifr≤30bulgee+ <dig> rtlnr30otherwise.
if  closes a left bulge of size r in s
<dig>  then the bulge is closed by base pair  involving nucleotide pair u,v, and
<dig> z∗=e−eauirt·∑uv∈ℬexp−eaui,j,x,yrt×exp−bulgert·4r·z∗i+r+ <dig> j− <dig> u,v
while if  closes a right bulge in s
<dig>  then the bulge is closed by base pair  involving nucleotide pair u,v, and
<dig> z∗=e−eauirt·∑uv∈ℬexp−eaurt×exp−bulgert·4r·z∗i+ <dig> j−r− <dig> u,v
internal loop
in this case,  closes an internal loop in s
<dig>  whose left  portion is of size r
<dig> [resp.
r
2].
since internal loop size r=r
1+r
<dig> may exceed the values in table look-up, we define the free energy for an internal loop of size r by
<dig> internal=internaleifr≤30internale+ <dig> rtlnr30otherwise.
the closing base pair  of the internal loop of size r=r
1+r
<dig> may involve the nucleotides uv∈ℬ, while the unpaired  nucleotides in positions i+ <dig> j− <dig> i+r
<dig> j−r
<dig> may involve a,b,c,d∈n.
in addition, there is an energy penalty for non symmetric internal loops, min, where the value of the constants asym and maxasym are given in the turner energy model.
thus
z∗=e−eauirt·exp−minasym·|r1−r2|,maxasymrt×∑uv∈ℬ∑a,b,c,d∈nexp−eaui,j,x,yrt×exp−internalr1+r2rt·4r1+r2−4×exp−mismatch+mismatchrt
<dig> ×z∗i+r1+ <dig> j−r2− <dig> u,v
external loop
despite the fact that, by following the total order on base pairs defined in eq.
, the dual partition function of multiloops is always computed before the dual partition function of the external loop, the computation of the dual partition function of multiloops will be easier to understand if the dual partition function of the external loop is defined in advance.
in order to improve speed, some implementations of rna thermodynamics-based algorithms ignore the contribution of dangling positions, which corresponds to vienna rna package -d <dig> flag.
rnadualpf also includes this option, which dramatically increases the speed of the algorithm.
the reason behind this difference of performance is clear from the following definitions.
suppose that h= constitutes the list of k external base pairs of s
<dig>  where i
1<j
1<i
2<j
2< ⋯ < i
k < j
k. for each , with 1≤r≤k, and for each choice of base pair gc, cg, au, ua, gu, ug, the value z
∗ has been previously computed and stored by dynamic programming, as well as the sum z
∗.
when the contribution of dangles is ignored, the dual partition function of an external loop with ℓ nucleotide positions external to every base pair is defined by
<dig> z∗=4ℓ·∏r=1kz∗
where ℓ=n−∑r= <dig> …,k and n is the length of the target structure s
<dig>
the default treatment of dangles in rnadualpf described below corresponds to vienna rna package -d <dig> flag, where both flanking positions of each external base pair contribute to the free energy.
let d=[a
<dig> b
<dig> …,a
k,b
k]⊆[i
1− <dig> j
1+ <dig> ⋯,i
k,j
k] be a list of those nucleotide positions that are adjacent to the k external base pairs ,…,.
the ordered multiset [a
<dig> b
<dig> …,a
k,b
k] can be considered as a collection of constraints, so that  if a
2=i
2− <dig>  and a
2=j
1+ <dig>  then a
2=b
<dig> and any nucleotide value that is assigned to b
<dig> must simultaneously be assigned to a
<dig>  moreover, there can also be an overlap between the list of base paired positions in h [i
<dig> j
<dig> …,i
k,j
k] and the multiset d=[a
<dig> b
<dig> …,a
k,b
k].
if  j
1=i
2− <dig>  then b
1=i
<dig> and a
2=j
<dig>  therefore, in the computation we have to account for these constraints.
let m denote the number of unpaired positions in d, without repetitions, and define a
r,b
r as the nucleotides instantiated respectively at a
r,b
r. the energy term for a 5′-dangle  on base pair  with nucleotides u,v is denoted by e
d <dig> [resp e
d3] where the dangle position x− <dig>  is assigned nucleotide w. with the notation just described, we have
<dig> z∗=∑〈u <dig> v <dig> …,uk,vk〉∈ℬk∑a <dig> b <dig> …,ak,bk∈n2k4ℓ−m×∏r=1kz∗ir,jr;ur,vr×exp−ed5ir,jr,ar;ur,vr,ar+ed3ir,jr,br;ur,vr,brrt
depending on the target structure s
<dig>  it can happen that the second sum of eq.
must be restricted to range over strictly less than 42k many rna sequences.
this is explained as follows.
if i
1= <dig> [resp.
j
r=n] then there is no position for a 5′  dangle, and hence the nucleotide sequences considered in the second summation would have length strictly less than 2k.
moreover, certain 5′ dangled positions could be identical to 3′ dangle positions, which arises for instance when j
k+2=i
k+1; alternatively, certain dangled positions could be identical with base-paired positions, which arises for instance when j
k+1=i
k+ <dig>  in such situations, instantiations of the 3′-dangle on  and the 5′-dangle on  are not independent, thus leading to a restriction of the range of the second summation in eq.
.
a similar restriction is implicitly assumed in the treatment of external loops in this section and of multiloops in the next section.
the algorithm performance can be improved by dividing the external loop into groups of components having interdependently constrained dangling positions, as just explained.
define two base pairs , as adjacent if x<y<x
′<y
′ and x
′−y≤ <dig> – i.e.
dangling positions of the base pairs , are constrained.
let g denote a maximal collection of adjacent base pairs belonging to h=, together with their associated dangle positions in d=[i
1− <dig> j
1+ <dig> …,i
k− <dig> j
k+1].
it is important to note that h∪d is thus partitioned into a collection of g disjoint groups g=.
therefore, we can divide an external loop of k helices into a collection groups g of size g≤k, and p unpaired positions that are external to every base pair of s
<dig> and not adjacent to any base pair.
for a group g with h base pairs, let h= denote the list of base pairs in g, and let d=[α
<dig> β
<dig> …,α
h,β
h]⊆[κ
1− <dig> λ
1+ <dig> ⋯,κ
h− <dig> λ
h+1] denote their associated dangle positions.
if u
r,v
r,a
r,b
r denote the nucleotides instantiated at the base pair r= and its respective dangling positions α
r,β
r respectively, then the dual partition function of g is the following.
<dig> z∗=∑〈,…,〉∈ℬh∑a <dig> b <dig> …,ah,bh∈n2h×∏r=1hz∗κr,λr;ur,vr×exp−ed5κr,λr,αr;ur,vr,ar+ed3κr,λr,βr;ur,vrbrrt
where the range of the second summation can be constrained by the overlap among positions in d and between positions in d and h, as explained for eq.
.
finally, since there are no shared dangling positions between groups, the dual partition function of an external loop is defined by
<dig> z∗=4p·∏r=1gz∗.
multiloop
suppose that  closes a multiloop in s
<dig>  which is a k-loop, or -way junction, for k> <dig>  where there are ℓ unpaired bases in the multiloop.
suppose that the k components of the multiloop are closed by the base pairs ,…, with the property that i<i
1<j
1<i
2<j
2<⋯<i
k<j
k<j.
assume that for all nucleotide choices in ℬ for each of the k base pairs of the multiloop , for 1≤r≤k, the value z
∗ has previously been computed and stored by dynamic programming, as well as the sum z
∗.
the computation of the dual partition function is similar to that of the external loop.
however, in this case we have to add the contribution of the base pair closing the multiloop , the au-penalties applied to this base pair, and the energetic penalty of a multiloop a+b·+c·ℓ, where the values of the constants a, b and c are given in the turner energy model.
then, the dual partition function of a multiloop without accounting for dangling positions is
<dig> z∗=e−eauirt·exp−a+b·+cℓrt·4ℓ×exp−eaurt·∑〈u <dig> v <dig> …,uk,vk〉∈ℬk×∏r=1kz∗ir,jr;ur,vr
the notation we use to define the dual partition function of multiloops with dangling positions is similar to that described for external loops.
however, some modifications are required in the previously given definitions, since we have to take into account the flanking positions of the base pair  closing the multiloop.
let h= be the collection of k base pairs closing one of the k components of the multiloop, and the base pair  closing the multiloop, and define the multiset d=[a
<dig> b
<dig> …,a
k+ <dig> b
k+1]⊆[i
1− <dig> j
1+ <dig> ⋯,i
k− <dig> j
k+ <dig> i+ <dig> j−1] of nucleotide positions adjacent to the base pairs in h. due to the possible overlap with the base pair closing the multiloop and its flanking positions, there are additional constraints in the ordered multiset [a
<dig> b
<dig> …,a
k+ <dig> b
k+1], so that  if a
1=i
1− <dig>  and i
1=i+ <dig>  then a
1=a
k+ <dig> and any nucleotide value that is assigned to a
<dig> must simultaneously be assigned to a
k+ <dig>  moreover, there can also be an overlap between the list of base paired positions [i
<dig> j
<dig> …,i
k,j
k,i,j] and the multiset [a
<dig> b
<dig> …,a
k+ <dig> b
k+1].
if  i=i
1− <dig>  then a
k+1=i
<dig> and a
1=i.
let m denote the number of unpaired positions in d, without repetitions.
then, the dual partition function of a multiloop with dangling positions is defined as follows.
<dig> z∗=e−eauirt·∑〈u <dig> v <dig> …,uk,vk〉∈ℬk∑a <dig> b <dig> …,ak+ <dig> bk+1∈n2exp−a+b·+cℓrt·4ℓ−m·exp−eaui,j,x,yrt×∏r=1kz∗ir,jr;ur,vr×exprt×exp−ed3j,i,ak+1;y,x,ak+1+ed5j,i,bk+1;y,x,bk+1rt
as explained for eq.
, it can happen that the second summation must be restricted to range over strictly less than 42k many rna sequences.
a decomposition similar that for external loops can be performed to improve the performance in the computation of the dual partition function of a multiloop.
in a multiloop, in addition to the adjacency definition given for external loops, we consider the base pair  that closes the multiloop as adjacent to a base pair  that closes a component of the multiloop, where i<x<y<j, if either x≤i+ <dig> or y≥j− <dig>  then, let g denote a maximal collection of adjacent base pairs belonging to h=, together with their associated dangle positions in d=[i
1− <dig> j
1+ <dig> …,i
k− <dig> j
k+ <dig> i+ <dig> j−1].
this decomposition produces a collection g of g disjoint groups g
<dig> …,g
g, one of which, designated the closing group
g
c contains the closing base pair  of the multiloop, and g− <dig> of which, designated as non-closing groups
g
nc, do not contain the base pair .
non-closing groups have the same composition as those defined for external loops – i.e.
a collection of h base pairs h= and a set of dangling positions d=[α
<dig> β
<dig> …,α
h,β
h]⊆[κ
1− <dig> λ
1+ <dig> ⋯,κ
h− <dig> λ
h+1].
therefore, we can compute the dual partition function
z of a non-closing group as described in eq.
.
in addition, the collection of non-closing groups of size g− <dig> of a multiloop of k components is denoted by gnc, where 0≤≤k.
therefore, a multiloop of k components and ℓ unpaired positions can be decomposed into one closing group g
c, a collection of non-closing groups gnc, and p unpaired positions that are not adjacent to any base pair, with 0≤p≤ℓ.
in a non-closing group, the collection of base pairs of size h+ <dig> is denoted by h=, where the base pair  closing the multiloop is at the last position.
the ordered multiset of adjacent positions is denoted by d=[α
<dig> β
<dig> …,α
h+ <dig> β
h+1]⊆[κ
1− <dig> λ
1+ <dig> ⋯,κ
h− <dig> λ
h+ <dig> i+ <dig> j−1], where the positions adjacent to i and j are at the last positions are respectively denoted by α
h+ <dig> β
h+ <dig>  a graphical example of a closing group and a non-closing group is shown in fig.
1
e, where the positions of a non-closing group with  <dig> base pair are highlighted in green and the positions of the closing group are highlighted in red and blue, and where the base pair  that closes the multiloop is depicted in red.
fig.
<dig> sampling dependency examples in rnadualpf for different structural elements:  stacked base pair,  hairpin,  1× <dig> internal loop,  3× <dig> internal loop and  multiloop.
base pair  to be sampled is highlighted in red, positions whose energy contribution is dependent on the instantiation of  are highlighted in blue, and positions that are mutually dependent, but independent of the instantiation of , are highlighted in green.
unpaired positions where the nucleotide choice has no effect in the free energy of the structure are indicated in black
for a closing group
g
c with h+ <dig> base pairs in h= and their flanking positions d=[α
<dig> β
<dig> …,α
h+ <dig> β
h+1]⊆[κ
1− <dig> λ
1+ <dig> ⋯,κ
h− <dig> λ
h+ <dig> i+ <dig> j−1], let x,y denote the nucleotides assigned to the closing base pair of the multiloop , and let u
r,v
r,a
r,b
r denote the nucleotides assigned respectively to the base pair r= and its flanking positions α
r,β
r. then, the the dual partition function
z
∗ of the closing group is defined by
<dig> e−eauirt·∑〈u <dig> v <dig> …,uk,vh〉∈ℬh∑a <dig> b <dig> …,ah+ <dig> bh+1∈n2×exp−eaui,j,x,yrt·∏r=1hz∗κr,λr;ur,vr×exp−ed5κr,λr,αr;ur,vr,ar+ed3κr,λr,βr;ur,vr,brrt×exp−ed3j,i,αh+1;y,x,ah+1+ed5j,i,βh+1;y,x,bh+1rt
in the same way as in eq.
, the values of the second summation are constrained to the possible choices among overlapping positions.
then, the dual partition function
z
∗ of the multiloop with k components and ℓ unpaired positions, where p of which are not adjacent to any base pair, is defined by
<dig> z∗=exp−a+b·+cℓrt·4p×z∗gc;x,y·∏gnc∈gncz∗gnc
sampling
once the dual partition function
z
∗ and its subcases z
∗ for each base pair  have been computed, it is possible to perform a boltzmann weighted sampling of positions i and j. for example, given the target structure with sequence constraints depicted in fig.
<dig>  rnadualpf computes the dual partition function table shown in table  <dig>  the dual partition function of the substructure enclosed by the base pair  is z
∗, and the dual partition function of the substructure enclosed by the base pair  where i,j are currently instantiated by the nucleotides x,y is denoted by is z
∗.
therefore, the boltzmann probability of x,y at positions i,j in the substructure enclosed by the base pair  is z
∗/z
∗ and can be sampled using the roulette wheel method.
fig.
<dig> target structure with sequence constraints used as input of rnadualpf to compute the dual partition function values shown in table  <dig>  sequence constraints are highlighted in red
due to the turner energy model, it is necessary to determine nucleotide positions whose instantiation influences the energy  of other positions, and subsequently all mutually dependent positions must be instantiated simultaneously.
figure  <dig> illustrates the mutual dependencies that must be considered when sampling different types of elements, where the base pair  to be sampled is highlighted in red, positions whose sampling probability is dependent on the instantiation of  are highlighted in blue, and positions that are mutually dependent, but independent of the instantiation of , are highlighted in green.
since the dynamic programming algorithm for the dual partition function proceeds from inner to outer base pairs, using the total ordering ≺ in eq.
, the sampling order of base pairs proceeds from outer to inner positions, i.e.
from largest base pair index to smallest.
in order to account for mutual dependencies in the sampling step, we define the function sample for each base pair  in s
<dig>  where k indicates the base pair index defined from eq.
, t indicates the type of structural element closed by base pair  in the target rna secondary structure, as shown in table  <dig>  and x,y are the instantiated nucleotides at positions .
due to the mutual dependencies, sampling a base pair with base pair indexk closing an m-loop, for m> <dig>  forces the instantiation of all inner closing base pairs of the m-loop, and the base pair index of each such inner base pair is strictly less than k. for this reason, except in the case of external loops, the outermost base pair  has been always instantiated before sample is called, and therefore the instantiation x,y is given as a parameter of the sampling function.
the boltzmann probability of each possible instantiation of mutually dependent positions can be computed on the fly in the backward step.
however, in order to improve the speed of the algorithm, in the forward step rnadualpf stores  the conditional dual partition function values of instantiations of interdependent positions.
these tables are used by the sampling function, since each value corresponds to the dual partition function conditional on a specific instantiation of the positions to be sampled by sample.
since the sampling procedure depends on the type t of element, we describe the function sample for each type of element – hairpin, stacked base pair, internal loop , multiloop and external loop as depicted in fig.
<dig>  for each of these cases, the values are stored in the conditional dual partition function table associated with the closing base pair .
hairpins
when hairpin size exceeds three , since the base pair  has been previously instantiated, flanking positions i+ <dig> j− <dig> are sampled first.
given the current assignment x,y, the boltzmann probability of sampling respectively the nucleotides u,v at the flanking positions i+ <dig> j− <dig> is
<dig> pi+1=u,j−1=v|i=x,j=y=z∗i,j,i+ <dig> j−1;x,y,u,vz∗i,j;x,y
therefore, in the forward step rnadualpf stores in a table the conditional dual partition function of each possible instantiation {x,y,u,v} of the base pair  and its flanking positions i+ <dig> j− <dig> respectively, defined by
<dig> z∗i,j,i+ <dig> j−1;x,y,u,v=e−eauirt·exp−hrt·exp−mismatchrt·4j−i− <dig>
then, remaining unpaired positions are uniformly sampled, since the nucleotide choice does not change the final free energy.
triloops, tetraloops and hexaloops are exceptions to this rule, since there are special loops that contribute to or penalize the free energy.
in those cases, we have to account for the special loops, as defined in “hairpins” section.
although it could seem to be a waste of space to store a different conditional dual partition function table for each base pair , even for two different hairpins of the same size in the target structure, one should note that rnadualpf allows sequence constraints, and thus z
∗ could possibly differ from z
∗ when  and  close hairpins of the same size.
stacking base pairs
as depicted in fig.
1
b, sampling probability of a base pair with base pair index
k− <dig> is dependent on the value sampled at the adjacent stacking base pair with base pair indexk.
therefore, sample samples the base pair  using the conditional probability given the instantiation of base pair  by x,y, defined as follows:
<dig> pi+1=u,j−1=v|i=x,j=y=z∗i,j,i+ <dig> j−1;x,y,u,vz∗i,j;x,y
the conditional dual partition function values stored in the forward step correspond to each instantiation {x,y,u,v} of the base pairs ,, denoted by
<dig> z∗i,j,i+ <dig> j−1;x,y,u,v=e−eauirt·exp−stackrt×z∗i+ <dig> j− <dig> u,v
internal loops
the energy contribution of internal loops in the turner energy model depends on the flanking unpaired positions of both the inner and outer closing base pairs, hence the sampling probability of the inner base pair cannot be separated from the adjacent unpaired positions.
moreover, for specific sizes of internal loop , the inner and outer closing base pairs share flanking positions.
in these cases, all the unpaired positions and the outer base pair must be sampled at the same time, since the energy contribution of each combination of base pairs and flanking positions is different.
in the 1× <dig> internal loop depicted in fig.
1
c, if the outer base pair  is instantiated by x,y, then let p denote the probability of sampling the nucleotides u,v,a,b,c respectively at positions k,l,n
<dig> n
<dig> n
<dig>  where  is the inner closing base pair, n
<dig> is the flanking position at i+ <dig> shared by the base paired positions i and k, and n
<dig> and n
<dig>  in the following equation, let p abbreviate the conditional probability just defined.
then
<dig> pu,v,a,b,c|x,y=z∗i,j,k,l,n <dig> n <dig> n3;x,y,u,v,a,b,cz∗i,j;x,y
rnadualpf computes and stores the conditional dual partition function of each possible instantiation {x,y,u,v,a,b,c} respectively at positions i,j,k,l,n
<dig> n
<dig> n
<dig>  where the value z
∗ is defined by
<dig> e−eauirt·exp−min−|,maxasym)rt×4j−l−3·exp−eaurt×exp−internalk−i+j−l−2rt×exp−mismatch+mismatchrt×z∗
for internal loops of sizes  similar conditional dual partition function tables are computed following the definitions in “stacked base pairs, bulges and internal loops” section.
other internal loops: when there are no shared flanking positions between the two base pairs that close an internal loop, as depicted in fig.
1
d, the energy contribution of innermost base pair and its respective flanking positions is independent of those of the outermost base pair.
in this case, rnadualpf samples first the flanking positions i+ <dig> j− <dig> of the outermost base pair , whose sampling probability is solely dependent on the instantiated nucleotides x,y at positions i,j.
is not necessary to store any conditional dual partition function for sampling these positions, since the probability of sampling the values a,b at the flanking positions i+ <dig> j− <dig>  given the assignment x,y is defined by
<dig> pi+1=a,j−1=b|i=x,j=y=exp−mismatchrt∑c,d∈nexp−mismatchrt
where mismatch penalties are obtained from table look-up.
finally, the innermost base pair  and its flanking positions k− <dig> l+ <dig> are sampled together.
in this case, we need to store an additional value z
∗, which is given by
<dig> z∗k− <dig> l+1=∑uv∈ℬ∑c,d∈nexp−mismatchrt·z∗
then, following the same notation, the probability of sampling the nucleotides v,u,d,c respectively at positions k,l,k− <dig> l+ <dig> is
<dig> pk=v,l=u,k−1=d,l+1=c=z∗k,l,k− <dig> l+1;v,u,d,cz∗k− <dig> l+ <dig>
therefore, the conditional dual partition function of each possible instantiation {v,u,d,c} stored in the corresponding table is defined as
<dig> z∗k,l,k− <dig> l+1;v,u,d,c=exp−mismatchrt·z∗
finally, since the remaining unpaired position does not contribute to the free energy, it is uniformly sampled.
multiloops and external loops
as explained in “external loop” section, if dangling positions are not included in the computation, sampling an external base pair or the closing base pair  of a multiloop from z
∗ is trivial.
on the other hand, by including dangling positions in the sampling, there is a dramatic increase in the space complexity of rnadualpf, albeit the space used is only a constant factor larger.
however, the decompositions into groups described in “external loop” and “external loop” sections allow to sample the positions of each group independently.
the example shown in fig.
1
e depicts a multiloop with two groups: a non-closing group
g
nc highlighted in green, and a closing group
g
c highlighted in red and blue, where the closing base pair of the multiloop  is marked in red.
in a non-closing group
g
nc all base pairs in h and dangling positions in d must be sampled together.
therefore, the conditional dual partition function of each possible instantiation of nucleotides at the h closing pairs in h and their adjacent positions in d is stored.
let u={u <dig> v <dig> …,uh,vh} denote an instantiation of the h base pairs in h=[κ
<dig> λ
<dig> …,κ
h,λ
h], and let w={a <dig> b <dig> …,ah,bh} denote an instantiation of the h flanking positions in d=[α
<dig> β
<dig> …,α
h,β
h] in the non-closing group
g
nc.
then, the probability of sampling u,w is
<dig> p=u,d=w)=z∗g,hgnc,dgnc;u,wz∗
therefore, the conditional dual partition function of each instantiation u,w at h,d, stored in the table of the group, is defined by
<dig> z∗g,hgnc,dgnc;u,w=∏r=1hz∗κr,λr;ur,vr×exp−ed5κr,λr,αr;ur,vr,ar+ed3κr,λr,βr;ur,vr,brrt
recall that the base pairs in h are adjacent.
therefore, due the constraints given by the overlapping positions within d, and between d and h, explained in “external loop” section, the number of possible instantiations u,w of h,d is ≤.
in a similar way, sampling from the closing group
g
c closed by the base pair , with h+ <dig> base pairs in h and their corresponding flanking positions in d requires us to store the conditional dual partition function of each instantiation of nucleotides {x,y,u,w} respectively at i,j,h,d, where u={u <dig> v <dig> …,uh,vh} denotes an instantiation of the h first base pairs  in h, w={a <dig> b <dig> …,ah+ <dig> bh+1} denotes an instantiation of the 2· flanking positions in d=[α
<dig> β
<dig> …,α
h+ <dig> β
h+1], and x,y denotes an instantiation of .
the probability of the instantiation u,w, given the nucleotides x,y is
<dig> ph=u,d=w|i=x,j=y=z∗gc,i,j,h,d;x,y,u,wz∗gc;x,y
then, the values stored in the table of the closing group correspond to the conditional dual partition function of each instantiation x,y,u,w are given by z∗,d;x,y,u,w), which is defined by the following expression:
<dig> e−eauirt·exp−eaurt·∏r=1hz∗κr,λr;ur,vr·exp−ed5κr,λr,αr;ur,vr,ar+ed3κr,λr,βr;ur,vr,brrt·exp−ed3j,i,αh+1;y,x,ah+1+ed5j,i,βh+1;y,x,bh+1rt
as a final remark, we would like to recall that all the conditional dual partition function values are computed and stored in the forward step at the same time as the dual partition function.
therefore, despite the consequent increase of space complexity in the algorithm, the computation of the values required for correct sampling does not involve a greater time complexity.
scaling
the sequence partition function z
∗ grows much faster than the usual structure partition function z, and so scaling must be used in the implementation.
let c> <dig> be a user-defined constant.
by a slight modification of the previous recursions, we actually compute z†=z∗cj−i+ <dig>  and hence z†=z∗cn, where n is the length of s
<dig>  for instance, the analogue of eq.
is
<dig> z†=z∗i,j;x,ycj−i+1=e−eauirt·exp−hrtcj−i+1×∑n <dig> n2∈nexp−mismatchrt·4j−i− <dig>
and the analogue of eq.
is
<dig> z†=z∗cj−i+1=e−eauirt·12·∑uv∈ℬexp−stackrt×z†i+ <dig> j− <dig> u,v
this modification does not affect properties of sequences sampled from the low energy ensemble, since the same scaling factor appears in both the numerator and denominator of all conditional probabilities.
for instance, the analogue of eq.
is
<dig> pi+1=u,j−1=v|i=x,j=y=z∗i,j,i+ <dig> j−1;x,y,u,vz∗i,j;x,y=z†i,j,i+ <dig> j−1;x,y,u,vz†
controlling gc-content
the gc-content of an rna sequence a=s
<dig> …,s
n is the number of nucleotides that are either g or c. instead of computing z
∗ and z
∗, we can compute z
∗ and z
∗, defined to be the corresponding partition dual partition functions, restricted to sequences having gc-content of α. note well that gc-content α includes the closing nucleotides x and y respectively located at positions i and j; i.e.
<dig> z∗i,j;x,y;α=∑ai,…,aj,gc=αai=x,aj=y,ai+ <dig> …,aj−1∈nexp−eai,…,aj;s0/rt
where s
<dig> denotes the restriction of target structure s
<dig> to the interval .
we describe two particular subcases, to provide the idea of how modifications need to be undertaken.
triloop
note that the number of rna sequences of length m having gc-content of α is mα·2α·2m−α=mα·2m≤4m, since α selected positions must be either g or c, yielding the term 2α, while the remaining m−α positions must be either a or u, yielding the term 2m−α.
assume that γ=|{x,y}∩{g,c}|=β.
then
<dig> z∗i,j;x,y;α=e−eauirt·exp−h+eaurt×j−i−1α−β·2j−i−1−|triloopx,y|+∑abc∈triloopx,yγ=α−βexp−triloopert
multiloop and external loop
assume that  closes a multiloop, which is a -way junction with ℓ unpaired nucleotides.
assume that the ordered multiset of potential dangle positions is d=[a
<dig> b
<dig> …,a
k+ <dig> b
k+1], where a
r=i
r− <dig> and b
r=j
r+ <dig> for r= <dig> …,k, and a
k+1=i and b
k+1=j, and assume that there are m unpaired positions that are not adjacent to a base pair in the multiloop.
if r denotes an rna sequence of arbitrary length, then let the function γ denote the gc-count in r. given an assignment of nucleotide base pairs u
1
v
<dig> …,u
k
v
k to ,…,, where u
r
v
r∈{gc,cg,au,ua,gu,ug}, and given an assignment a
<dig> b
<dig> …,a
k,b
k of dangle nucleotides, where ar,br∈n, for r= <dig> …,k, we let
<dig> γ=γa <dig> …,ak,b <dig> …,bk.
then the dual partition function of a multiloop with a gc-content of α is defined by setting z
∗ equal to the following:
<dig> e−eauirt·∑α1+⋯+αk≤α∑{ur,vr∈ℬ:r= <dig> …,k}∑{a <dig> b <dig> …,ak,bk∈n2k}exp−a+b·+cℓrt·α−∑r=1kαr−γ·2ℓ×exp−eaurt·∏r=1kz∗ir,jr;ur,vr;αr×exp−ed5ir,jr,ar;ur,vr,ar+ed3ir,jr,br;ur,vr,brrt×exp−ed3j,i,ak+1;y,x,ak+1+ed5j,i,bk+1;y,x,bk+1rt
since the modification required in the remaining cases follows similar reasoning as in the treatment of the hairpin and external loop just described, the details for these remaining cases are not given.
an additional challenge of computing the dual partition function with gc-content control is the combinatorial problem of efficiently counting the number n of instantiations of the external loop, consisting of all positions external to every base pair, with gc-content k, where the user can stipulate that certain positions are constrained to contain nucleotides consistent with iupac codes.
to this end, we implemented the combinatorial algorithm defined in supplementary information.
sampling with gc-content
the implementation of sampling with gc-content is performed in a similar manner as described in “sampling” section, with some notable differences.
first, the sampling function is redefined by sample, where k indicates the base pair index in the ordering defined by eq.
for the base pair  that is already instantiated by nucleotide pair xy, and t designates the type of structural element closed by base pair  in the target rna secondary structure, as shown in table  <dig>  the function sample instantiates all positions of the loop having outer closing base pair , including its inner closing base pair and which returns the gc-content of the sampled loop.
moreover, the gc-content of the subsequence a= will be α once the entire sequence a
<dig> …,a
n is sampled.
second, rnadualpf stores a conditional dual partition function table for each base pair  and gc-content  <dig> to j-i- <dig>  the function sample samples from the conditional dual partition function of those sequences which have exactly α gs and cs strictly between the positions i and j, thus guaranteeing a gc-content of α for the subsequence a once the entire sequence a
<dig> …,a
n is sampled.
note that sample samples only the loop closed by the already instantiated outer base pair , and that α is the gc-content of the entire subsequence a=a
i+ <dig> …,a
j− <dig> once the algorithm terminates.
only in the case that base pair  closes a hairpin loop will it generally happen that the gc-content of the loop closed by  is equal to α.
let α be the user-designated gc-content of sequences a=a
<dig> …,a
n to be sampled from a target secondary structure having ℓ base pairs.
the following pseudocode describes how to sample sequences a=a
<dig> …,a
n, whose gc-content exactly equals α. here, an external loop with m components means that there are m exterior base pairs ,…, such that all positions exterior to these base pairs are unpaired; i.e.
each position r∈{ <dig> …,n}−∪c=1m{ic,…,jc} is unpaired in the target structure.
to clarify how the gc-content is sampled in a statistically rigorous manner, suppose that the user has specified the gc-content to be α, and that l is the external loop of the target structure s
<dig> having m components, where the cth component has external closing base pair .
in computing the dual partition function, for all possible choices of non-negative integers α
<dig> …,α
m,β that sum to α and all 6m possible assignments of watson-crick or wobble nucleotide pairs x
<dig> y
<dig> …,x
m,y
m to the base pairs ,…,, the software rnadualpf has computed the sum of ∑c=1mzc∗ plus the boltzmann factor of the external loop with gc-content β. since the dual partition function z
∗ is the sum, taken over all values of α
<dig> …,α
m,β and all watson-crick and wobble pair assignments to the external base pairs, rnadualpf can then use the roulette wheel method to sample values α
<dig> …,α
m,β and x
<dig> y
<dig> …,x
m,y
m in a statistically rigorous manner.
multiloops, and other structural elements, which contain unpaired regions whose sequence does not contribute to the free energy of the structure, are handled in a analogous manner.
RESULTS
robustness and plasticity of c. elegans mirnas and e. coli sncrnas
in  <cit>  borenstein and ruppin used version  <dig>  of the vienna rna package  <cit>  to generate  <dig> rna sequences per wild type precursor microrna  extracted from the database rfam  <dig>   <cit> , with the property that each of the  <dig> control sequences folded into the wild type pre-mirna structure – i.e.
the minimum free energy  structure of each of the  <dig> control sequences was identical to the mfe structure of the wild type pre-mirna.
based on these computational experiments, borenstein and ruppin asserted that the “structure of mirna precursor stem–loops exhibits a significantly high level of mutational robustness in comparison with random rna sequences with similar stem–loop structures”.
noting that the vienna rna package inverse folding program rnainverse does not control for gc-content or other sequence compositional bias, the authors performed a second computational experiment, in which control sequences not only folded into the target wild type structure, but also had similar dinucleotide composition to that of wild type pre-mirna .
since the filtering step required enormous run time and computational resources, the authors restricted their attention to a small set of  <dig> micrornas, and generated only  <dig> control sequences per microrna – note here that rnainverse cannot control for gc-content.
borenstein and ruppin concluded that robustness of precursor micrornas was not the byproduct of a base composition bias or of thermodynamic stability.
subsequently rodrigo et al.
<cit>  undertook a similar analysis for bacterial small noncoding rnas , also using the program rnainverse, albeit using somewhat different definitions – see precise definitions in “formal definitions of robustness” section.
the main finding of  <cit>  was that bacterial sncrnas are not significantly robust when compared with  <dig> sequences having the same structure, as computed by rnainverse; however, the authors found that bacterial sncrnas tend to be significantly plastic, in the sense that the ensemble of low energy structures are structurally diverse.
unlike the case of precursor micrornas  <cit> , rodrigo et al.
did not control for sequence compositional bias.
using rnadualpf, we performed similar computational experiments on  <dig> precursor micrornas of c. elegans from mirbase  <dig>  <cit>  and for the bacterial small noncoding rnas of  <cit> .
below, we discuss each case separately.
for each c. elegans pre-mirna, we used rnadualpf to sample  <dig> sequences with no control over gc-content and  <dig> sequences whose gc-content was identical with that of the wild type pre-mirna.
moreover, each control sequence approximately folded into the mfe wild type pre-mirna structure as computed by vienna rna package  <dig> .
<dig>  <cit> .
table  <dig> shows that length-normalized base pair distance between the mfe structure of the control sequence and that of the pre-mirna is on average  <dig> ± <dig>  for default use of rnadualpf with control over gc-content, and  <dig> ± <dig>  when gc-content of each control sequence is identical to that of the corresponding wild type pre-mirna.
additional measures in table  <dig> show that sequences sampled from rnadualpf  are only modestly more stable thermodynamically,  the ensemble of low energy structures of control sequences deviate slightly more from the target pre-mirna structure, as is the case for wild type pre-mirna sequences, as mesured by ensemble defect  <cit> , expected base pair distance to target  <cit> , expected proportion of native contacts , average positional entropy  <cit> , morgan-higgs structural diversity  <cit> , and vienna structural diversity .
for each of the  <dig> wild type  pre-mirna sequences, rnadualpf sampled sequences, either having exactly the same gc-content as the wt sequence  or with no control over gc-content .
the designation ‘mfe’ indicates that the sampled sequences were subsequently filtered to retain only those, whose minimum free energy structure is identical to the mfe structure of the corresponding wt pre-mirna; otherwise, the designation ‘def’ is used to indicate the default output of rnadualpf, without the subsequent filtering step.
for each wt pre-mirna sequence, rnadualpf generated  <dig> sequences for the default case def , and  <dig> sequences for the non-default case mfe, such that sample mfe structure is identical to wt mfe structure.
various measures were used to compare the properties of rnadualpf sampled sequences to those of wild type sequences: bp dist target: length-normalized average base pair distance d
bp between the mfe structure s
<dig> of sequences sampled by rnadualpf and the target structure s
∗.
energy mfe: length-normalized average free energy e of mfe structure s
<dig>  energy target: length-normalized average free energy e of target s
∗ for the respective sequences.
ensemble defect: length-normalized expected hamming distance to target s
∗  <cit> .
exp bp dist: length-normalized expected base pair distance to target s
∗  <cit> .
prop nat contact: expected proportion of base pairs of target s
∗ that occur in the mfe structure, i.e.
〈|s0∩s∗||s∗|〉.
pos entropy: average positional entropy  <cit> .
gc content: average proportion of positions occupied by g or c. ln dual prob: average natural logarithm of the dual probability exp/rt)/z∗ that sequence a adopts the structure s. ln prob: average natural logarithm of the probability exp/rt)/z that sequence a adopts the structure s. mh str div: length-normalized morgan-higgs structural diversity  <cit> .
vienna str div: length-normalized vienna structural diversity, called ensemble diversity in  <cit> .
values of all measures for default sampled sequences having gc-content within  <dig> % of wild type gc-content  are essential identical to those of exact gc-content control
tables  <dig> and  <dig> display a similar analysis of the collection of bacterial small noncoding rnas of  <cit>  and of rfam  <dig>  database  <cit> .
for the rfam database, we selected one sequence from each of the ≈ <dig> rfam families, with the property that the mfe structure of the sequence most resembled the rfam consensus structure – i.e.
whose mfe structure has smallest base pair distance to the consensus structure.
these tables show similar trends as those displayed in table  <dig>  although values are larger due to increased sequence length of bacterial sncrna and sequences from rfam.
see table  <dig> for an explanation of column headers and various measures.
since bacterial noncoding rna is generally much longer than precursor microrna, no subsequent filtering step was undertaken to ensure that sample sequence mfe structure is identical to that of wild type pre-mirna.
however an additional column is given for sequences required by rnadualpf to have gc-content is within  <dig> % of wt value.
for each rna family from rfam  <dig> , we selected that sequence whose mfe structure had smallest base pair distance to the rfam consensus structure for the family.
these sequences constituted the collection wt.
see table  <dig> for an explanation of column headers and various measures
in agreement with  <cit> , the left panel of fig.
<dig> shows that c. elegans mirna is significantly robust , provided that gc-content is not controlled.
however, in contrast to  <cit> , when gc-content is controlled, we find that c. elegans mirna is significantly non-robust .
to corroborate our findings, for each wild type c. elegans pre-mirna, we performed a second computational experiment, to generate  <dig> sequences with no control over gc-content and  <dig> sequences whose gc-content was identical with that of the wild type pre-mirna.
in contrast to the first experiment, we used rnadualpf to generate sufficiently many sequences to subsequently select  <dig> sequences  and  <dig> sequences , each of whose mfe structure was identical to that of wild type pre-mirna.
the left panel of fig.
<dig> shows that when gc-content is not controlled, c. elegans precursor micrornas are statistically robust , in agreement with the main result of borenstein and ruppin  <cit> .
however, when gc-content of control sequences is identical to that of wild type precursor microrna, we confirm that c. elegans pre-mirna is statistically non-robust .
note that our finding, which is in opposition to results of borenstein and ruppin  <cit> , is based on a larger data set of precursor micrornas, each of which has a larger control set, than in the analysis of  <cit> .
fig.
<dig> z-scores of mutational robustness  and of plasticity  are presented for the bacterial small noncoding rna  collection from  <cit>  and for c. elegans precursor microrna  from mirbase  <dig>  for each wild type pre-mirna  wild type sequence, rnadualpf sampled  <dig>  sequences using the minimum free energy structure of the wild type sequence as target structure.
the gc-content of the sampled sequences was either required to be exactly that of the wild type sequence, or not , as indicated in the legend.
sampled sequences were used to compute the mutational robustness and plasticity, as explained in the main text.
note that c. elegans mirna is significantly robust if gc-content is not controlled, but significantly non-robust if gc-content of rnadualpf samples is identical to that of wild type pre-mirna.
similarly, bacterial sncrnas are not significantly robust if gc-content is not controlled, but significantly non-robust when gc-content is identical to that of wild type sncrna.
for this figure, mutational robustness of rna sequence a=a
<dig> …,a
n is defined by 1−〈dbp〉n, where ensemble distance d
bp between two length n sequences a and b is defined in  <cit> , and the average ensemble distance from all single-point mutants of a is defined by 〈dbp〉=∑bdbp3n where the sum is taken over all single-point mutants b of a. we use this notation of mutational robustness, rather than the notion defined in  <cit> , since the latter notion is not a true metric, as explained in “formal definitions of robustness” section.
the plasticity p=〈dv〉n/2=∑i<jpi,jn is defined in  <cit>  as normalized ensemble diversity, where ensemble diversity  <cit>   d
v is defined by eq.
fig.
<dig> z-scores of mutational robustness  and of plasticity  are presented for c. elegans microrna from rfam  <dig> .
robustness and plasticity were measured as explained in fig.
<dig>  for each wild type  sequence,  <dig> sequences were sampled from rnadualpf both with and without control over gc-content , while  <dig> sequences were generated by rnadualpf, both with and without control over gc-content, to exactly fold into the target structure .
this was achieved by repeatedly sampling sequences in order to obtain a sequence, whose mfe structure was identical to the target structure.
the number of sequences necessary to sample in order to obtain one sequence that folds into the target structure was variable, depending on the target structure and gc-content – in many cases, only  <dig> samples were necessary per selected sequence, in some cases  <dig> samples were necessary.
and in one specific case  <dig> samples were required.
by control over gc-content, we mean that all sampled sequences have identical gc-content to the wild type sequence
turning now to the analysis of bacterial small noncoding rnas, we find that sncrnas are not significantly robust  when gc-content is not controlled, confirming a result from rodrigo et al.
<cit> .
however, when gc-content of the sequences sampled from rnadualpf is required to be identical to that of wild type sncrna, bacterial sncrnas are see to be significantly non-robust .
note here that rodrigo et al.
used rnainverse in their computational experiments, hence could not consider the case with control over gc-content.
the left panels of figs.
3
a and  <dig> summarize our findings that precursor micrornas  are significantly non-robust  with respect to a control set of  <dig>  sequences generated by rnadualpf with identical gc-content to that of the wild type sequence.
finally, in our analysis of plasticity, the right panels of figs.
<dig> and  <dig> show that both c. elegans and bacterial small noncoding rnas exhibit more plasticity when compared with control sequences for which gc-content is not controlled, as well as when compared with control sequences for which gc-content is identical to that of wild type sequences.
structural rna has higher free energy than expected
in figure  <dig> of  <cit> , we showed that the free energy e
<dig> of the minimum free energy  structure s
<dig> of e. coli val-trna , is much higher  than the average free energy 〈e〉 of over four million rnas having the same mfe structure s
<dig> as that of e. coli val-trna.
here, e. coli val-trna rv <dig> was selected, because its mfe structure s
<dig> is identical to the rfam consensus structure for trna family rf <dig>  this preliminary result suggests that naturally occurring transfer rnas may be under selective pressure to be only marginally thermodynamically stable.
since it took a number of days for rnaifold  <cit>  to return over four million solutions of the inverse folding problem for the trna target structure, we now describe how rnadualpf can be used to compute the boltzmann expected free energy of literally all sequences a
<dig> …,a
n with respect to an arbitrary target structure s
<dig>  in this manner, we confirm our preliminary finding concerning e. coli val-trna, and show that the folding energy of structural rna from the rfam database is much higher  than expected.
before presenting results, we need some definitions.
for the turner nearest neighbor energy model  <cit> , the free energy of a secondary structure s of an rna sequence a=a
<dig> …,a
n depends on the  temperature t
<dig>  to indicate this dependence, we write e, where in the sequel, t
<dig> will be designated as table temperature, i.e.
the temperature for which parameters from the turner energy tables are applied.
for an arbitrary, but fixed secondary structure s
<dig> of length n, the dual partition function at temperature t
<dig> is defined by
<dig> z=∑aexp−ea,s <dig> t0/rt
where the sum is taken over all rna sequences a=a
<dig> …,a
n of length n. note that t
<dig> indicates the  temperature at which the energy of a structure s
<dig> and nucleotide sequence a is evaluated using the turner parameters, while all other occurrences of the temperature variable are designated by t, which we call formal temperature.
the distinction between formal and table temperature is made to allow us to use finite difference approximations to derivatives with respect to the formal temperature when when we compute dual expected energy and dual conformational entropy below .
when table temperature t
<dig> equals formal temperature t, and the temperature is clear from the context, we write z
∗; if the target structure s
<dig> is also clear from the context, then we write z
∗.
a similar remark applies to the other thermodynamic functions p
∗,g
∗,〈e
∗〉,s
∗, which we now define.
the dual boltzmann probability
p
∗ is defined by
<dig> p∗a,s <dig> t <dig> t=exp−ea,s <dig> t0z∗s <dig> t <dig> t
the dual ensemble free energy
g
∗ is defined by
<dig> g∗=g∗=gs <dig> t <dig> t=−rtlnz∗s <dig> t <dig> t
where r≈ <dig>  cal/ is the universal gas constant.
the dual expected  energy 〈e
∗〉 is defined by
<dig> 〈e∗s <dig> t <dig> t〉=∑aea,s <dig> t0·pa,s <dig> t <dig> t
straightforward derivations analogous to those in  <cit>  yield the following expressions for dual expected energy 〈e
∗〉 and dual entropy
s
∗:
<dig> 〈e∗s <dig> t <dig> t〉=rt2·∂∂tlnz∗t=t <dig>
<dig> s∗s <dig> t <dig> t=〈ea,s <dig> t <dig> t〉−g∗s <dig> t <dig> tt
programs to compute dual expected energy 〈e
∗〉, dual conformational entropy
s
∗, and dual heat capacity
cp∗ are provided at our web site.
we do not elaborate further on dual entropy or dual heat capacity, since at the present time we have found no compelling applications.
c, and computed 〈e〉.
additionally, for each rfam family, we selected that sequence a
<dig>  whose minimum free energy  structure s
<dig> has smallest base pair distance to the consensus structure s
c. we computed the expected energy 〈e〉, as well as the free energies e and e. figure  <dig> displays box-and-whiskers plots for the fold change 〈e〉e for the consensus structure and the fold change 〈e〉e for the minimum free energy structure.
since the dual boltzmann probability p
∗ is generally larger for sequences a having higher gc-content , rnadualpf computes as well the dual partition function for gc-content k, defined by
<dig> z∗=∑asuch thatgc-content=kexp/rt)
fig.
<dig> analysis of expected free energy 〈e〉 for structures in rfam  <dig>   <cit> .
given a secondary structure s, the expected free energy of all sequences a with respect to s is defined by 〈e〉=∑ae·exp/rt)z∗, where z
∗ is the dual partition function defined in equation .
for each rfam family, we took the family consensus structure s
c, and computed 〈e〉.
additionally, for each rfam family, we selected that sequence a
<dig>  whose minimum free energy  structure s
<dig> has smallest base pair distance to the consensus structure s
c. the expected energy 〈e〉 was computed, as well as the free energies e and e. the fold change 〈e〉e for the consensus structure and the fold change 〈e〉e for the minimum free energy structure were computed.
the box-and-whiskers plots show the mean, 25th and 75th percentile, minimum and maximum values.
as indicated in the legend, these computations were performed either with respect to all sequences or with respect to all sequences having the same  gc-content.
these data clearly indicate that natural rna sequences, whose mfe structures most closely resemble the rfam consensus structures, have higher free energy than expected
in this fashion, we can exactly compute the dual expected energy 〈e
∗〉 of all sequences having gc-content k which approximately fold into target structure s
<dig>  tables  <dig>   <dig> and  <dig> analyze what we mean by approximately folding into the target structure – i.e.
sequences a are preferentially sampled when free energy e is low, hence have large dual boltzmann probability.
rnadualpf, even when exact gc-content is controlled, is faster than inverse folding programs by orders of magnitude, hence providing an effective alternative manner of solving inverse folding.
CONCLUSIONS
in this paper we describe the algorithm and software rnadualpf, which computes the dual partition function
z
∗, defined as the sum of boltzmann factors exp/rt) of all sequences
a with respect to the target structure s
<dig>  using rnadualpf, we efficiently sample rna sequences that  fold into s
<dig>  where additionally the user can specify iupac sequence constraints at certain positions, and whether to include dangles .
moreover, the user can require that all sampled sequences have a precisely specified gc-content, since, optionally, we compute the dual partition function
z
∗ simultaneously for all values k=g+c.
this sampling strategy is complementary to the use of rnaifold  <cit> , since it allows the study of the properties of long rna structures whose number of solutions for the inverse folding problem is astronomically large.
we use rnadualpf to corroborate previous studies  <cit>  using rnainverse  <cit> , by confirming that precursor micrornas are significantly mutationally robust when gc-content is not controlled.
however, in contrast to  <cit> , we find that precursor micrornas are significantly non-robust when gc-content is controlled.
we confirm and extend previous findings  <cit>  that bacterial small noncoding rnas display plasticity  and are not statistically robust, when gc-content is not controlled.
additionally, we obtain the new finding that when when gc-content is controlled, bacterial small noncoding rnas are significantly non-robust, as in the case of precursor micrornas.
one possible reason for the discrepancy between our results and those of  <cit>  could be related with the fact that the energy parameters of vienna rna package  <dig>   differ from those of vienna rna package  <dig> .
<dig> .
another possible reason is that the inverse folding solutions returned by the program rnainverse used in  <cit>  show a different bias than sequences returned by rnadualpf .
as mentioned in the introduction, there is a relation between our c program rnadualpf and the python program incarnation  <cit> , although our work is independent of that of reinharz et al.
<cit> .
incarnation is a weighted sampling algorithm that computes the dual partition function for a simple energy model, which only considers base stacking free energies – unlike rnadualpf, the program incarnation includes no energy contributions for hairpins, bulges, internal loops, multiloops, dangles, or mismatches.
if the user specifies a desired gc-content α, then incarnation does not compute the dual partition function for gc-content, but rather applies an adjustable heuristic so that after a suitable burn-in period, sequences tend to approximately have gc-content α. see table  <dig>  of  <cit>  for benchmarking results on rnadualpf and incarnation, which show conclusively that rnadualpf is not only faster, but its sequences have a higher probability of folding into the target structure, its sequences have a smaller gc-content in default mode, where gc-content is not controlled, etc.
our original motivation in designing rnadualpf was to generate an unbiased sample of near-solutions  to the inverse folding problem.
at present, it seems clear that no program can claim to generate an unbiased sample of inverse folding solutions, since  the solution space so large that this hypothesis cannot be tested by brute force methods, and  different inverse folding algorithms return solution sequences having different properties, as shown in table  <dig> of  <cit> .
nevertheless, in the same manner that structures sampled by the algorithm of ding and lawrence  <cit>  constitute an unbiased, representative set of low energy secondary structures for a given rna sequence, as implemented in sfold  <cit>  and rnasubopt -p  <cit> , the collection of rna sequences sampled by the algorithm rnadualpf constitute an unbiased, representative set of sequences having low energy with respect to a given target structure s
<dig>  although the minimum free energy structure of such sequences may indeed be distinct from s
<dig>  it is likely that the mfe structure and the target s
<dig> be similar, as shown in tables  <dig>   <dig> and  <dig>  moreover, fig.
<dig> presents relative frequency plots that suggest that when gc-content is controlled, the sequences returned by rnadualpf have similar properties to those of wild type sequences:  similar expected base pair distance to the wild type target structure  <cit> ,  similar ensemble defect to the target wild type structure  <cit> ,  similar positional entropy  <cit> ,  similar vienna structural diversity ,  similar morgan-higgs diversity  <cit> ,  similar expected proportion of native contacts .
these graphs were produced by using rnadualpf to sample  <dig>  sequences for each of the  <dig> c. elegans precursor micrornas from the mirbase  <dig> database  <cit> , in each of the following cases:  gc-content identical to that of the rfam sequence,  no control for the gc-content.
figure  <dig> presents additional data, computed in the same manner for c. elegans pre-mirna from mirbase  <dig>  showing that when gc-content is controlled, sequences sampled by rnadualpf satisfy the following:  the average length-normalized base pair distance between the minimum free energy and target structures is ≈ <dig> ,  wild type and rnadualpf sampled sequences have similar free energy with respect to the wild type target structure,  as well as similar minimum free energy,  similar dual probability, and  similar probability to wild type rna sequences.
taken together, this data shows that if gc-content is controlled, then rnadualpf returns sequences whose low energy structures tend to resemble the target structure.
figures  <dig> and  <dig> are similar to figs.
<dig> and  <dig>  except that for each c. elegans pre-mirna,  <dig> sequences were generated by rnadualpf, each of whose mfe structure is identical to the wild type target structure .
taken together, figs.
<dig>   <dig>   <dig> and  <dig> present convincing evidence that rnadualpf generates sequences that  fold into the user-specified target structure, hence supporting our finding that c. elegans precursor micrornas are statistically non-robust, contrary to the finding of  <cit> .
fig.
<dig> for each of the  <dig> c. elegans precursor micrornas from mirbase  <dig> and for each of the following cases , indicated in black, and , indicated in red, rnadualpf sampled  <dig> sequences without any subsequent filtering step.
case  - black lines: all rnadualpf sequences have gc-content exactly equal to that of the rfam sequence .
case  - red lines: rnadualpf was used in default mode, without controlling gc-content .
case  - blue lines: wild type  c. elegans data.
density plots are shown for  the expected base pair distance to target structure s
<dig>  <cit> ,  the ensemble defect to target structure s
<dig>  <cit> ,  the positional entropy  <cit> ,  vienna structural diversity ,  morgan-higgs diversity  <cit> ,  expected proportion of native contacts .
all measures were normalized by sequence length
fig.
<dig> additional measures for the data described in the previous fig.
<dig>  density plots are shown for  the base pair distance between the minimum free energy  structure and the target structure,  the gc-content,  the free energy e of the rna sequences a with respect to the target structure s
<dig>   the free energy e of each sequence with respect to its own minimum free energy  structure,  the log dual probability p∗=∑xexp/rt)z∗, and  the log probability p=∑sexp/rt)z
fig.
<dig> for each of the  <dig> c. elegans precursor micrornas from mirbase  <dig> and for each of the following cases , indicated in black, and , indicated in red, rnadualpf generated  <dig> sequences, whose minimum free energy structure was identical to that of the corresponding wild type pre-mirna .
case  - black lines: all rnadualpf sequences have gc-content exactly equal to that of the rfam sequence .
case  - red lines: rnadualpf was used in default mode, without controlling gc-content .
case  - blue lines: wild type  c. elegans data.
density plots are shown for  the expected base pair distance to target structure s
<dig>  <cit> ,  the ensemble defect to target structure s
<dig>  <cit> ,  the positional entropy  <cit> ,  vienna structural diversity ,  morgan-higgs diversity  <cit> ,  expected proportion of native contacts .
all measures were normalized by sequence length
fig.
<dig> additional measures for the data described in the previous fig.
<dig>  density plots are shown for  the base pair distance between the minimum free energy  structure and the target structure,  the gc-content,  the free energy e of the rna sequences a with respect to the target structure s
<dig>   the free energy e of each sequence with respect to its own minimum free energy  structure,  the log dual probability p∗=∑xexp/rt)z∗, and  the log probability p=∑sexp/rt)z
additionally, we have shown that natural rnas from the rfam  <dig>  database have higher minimum free energy than expected, thus supporting our results in  <cit>  which suggest that functional rnas are under evolutionary pressure to be only marginally thermodynamically stable.
the applications described in this paper demonstrate that rnadualpf is a useful and extremely fast software tool for evolutionary and synthetic biology.
endnote
<dig> when dangling positions are not included in the computation , the algorithm clearly requires linear time.
when dangling positions are included , run time is exponential in the number of components of the largest multilooop; however, in practice the algorithm is extremely fast, and it is possible to modify the algorithm to always run in linear time.
additional file
additional file  <dig> supplementary information.
title is supplementary information for rnadualpf: software to compute the dual partition function with sample applications in molecular evolution theory.
this is a 3-page pdf file containing a tricky derivation for an efficient computation of the number of external loops of size n with gc-content k, where the user can stipulate that certain positions are constrained to contain nucleotides consistent with iupac codes.
abbreviations
mfeminimum free energy
plmvdpeach latent mosaic viroid
sncrnasmall noncoding rna

 in, borenstein and ruppin define neutrality of an rna sequence a=a …,a n by η=1−〈d〉n, where in this section 〈d〉 denotes the average, taken over all 3n single-point mutants of a, of the base pair distance d bp between the minimum free energy structure s of a and the mfe structures of single-point mutants of a. an rna sequence a is then defined to be robust if η is greater than the average neutrality of control sequences generated by the program rnainverse, which fold into the same target structure s the main finding of is that precursor micrornas exhibit a significantly higher level of mutational robustness than random rna sequences having the same structure. to control for sequence composition bias in their computational study, the authors selected sequences from the output of rnainverse, whose dinucleotide composition was similar to that of wild type pre-mirna. since the filtering step required enormous run time and computational resources, the authors restricted their attention to a small set of micrornas, generating only control sequences per microrna. borenstein and ruppin conclude that robustness of precursor micrornas is not the byproduct of a base composition bias or of thermodynamic stability. subsequently rodrigo et al. undertook a similar analysis for bacterial small rnas, also using the program rnainverse, albeit using somewhat different definitions – precise definitions are given in “formal definitions of robustness” section. the main finding of was that bacterial sncrnas are not significantly robust when compared with sequences having the same structure, as computed by rnainverse; however, bacterial sncrnas tend to be significantly plastic, in the sense that the ensemble of low energy structures is structurally diverse. unlike the case of precursor micrornas, rodrigo et al. did not control for sequence compositional bias. this raises the question of whether the control sequences analyzed in are representative or to what extent features shared by sequences output by the program rnainverse are artifacts of the program used. indeed, the number of rna sequences that fold into a given target structure can be astronomically large. over a few weeks, before we elected to terminate the execution, our state-of-the-art inverse folding software rnaifold generated, many 52-nt sequences that fold exactly into the mfe secondary structure s of hiv- ribosomal frameshift stimulating signal from the gag-pol overlap region af /1631- and which additionally code 17-mer peptides in the gag and pol reading frames having amino acids that appear in gag/pol peptides found in the los alamos hiv- database. the number of nt rna sequences that fold into target s without additionally imposing the constraint of coding particular peptides in overlapping gag/pol reading frames is certain to dwarf the previous number. moreover, the number of sequences that fold into the mfe structure of an animal precursor microrna or into the mfe structure of bacterial sncrna is certain to be even more daunting. different inverse folding algorithms have adopted different strategies to generate sequences that fold into a user-specified target secondary structure s for instance, rnainverse performs an adaptive walk, in one step of which a nucleotide in the current sequence is mutated and subsequently accepted if the base pair distance between the minimum free energy structure of the mutated sequence and the target structure s is reduced. nupack design selects a candidate mutation position with probability proportional to its contribution to the ensemble defect ∈s and s 0=i indicates i is unpaired in s 0). rnaifold cp-design uses constraint programming to systematically explore the search tree of all inverse folding solutions in an order determined by certain heuristics. accordingly, one cannot claim that the collection of sequences generated by any particular inverse folding algorithm is representative of the astronomically large space of all inverse folding solutions – indeed, each inverse folding algorithm has an inherent but unknown bias. in this paper, we describe the algorithm rnadualpf, which generates sequences which have low free energy with respect to a user-specified target structure s – i.e. the inherent bias of rnadualpf is known, unlike the situation of other inverse folding algorithms. we show that rnadualpf is extremely fast software for generating sequences that approximately fold into s 0; moreover, in a postprocessing step, one can filter the output of rnadualpf to select sequences that exactly fold into s rnadualpf additionally allows the user to specify iupac codes to constrain certain nucleotide positions as well as to control the gc-content of all generated sequences. sampling is performed in a manner distinct but somewhat analogous to that by which sfold and rnasubopt -p sample representative secondary structures from the boltzmann ensemble of all structures of a given sequence. using rnadualpf, we perform a pilot study that is similar, though not identical, to that of for two classes of rna: c. elegans precursor microrna from mirbase and the bacterial small noncoding rnas previously analyzed in. finally, it should be noted that, although rnadualpf was developed entirely independently of the work of reinharz et al., one can view our c-program as an extension of python program incarnation to the full turner energy model, where additionally gc-content is rigorously handled. this point will be discussed further in the conclusion. formal definitions of robustness let a=a …,a n denote an arbitrary rna sequence, where ai∈n={a, u, g, c}, a secondary structure s of a is a set of base pairs satisfying the following conditions: if ∈s then a i,a j constitute a watson-crick or gu wobble pair, i.e. ij∈ℬ which is the set {au,ua,gc,cg,gu,ug}. if ∈s then i+θ<j, where θ=. if ∈s and ∈s, then either i<k<ℓ<j or k<i<j<ℓ or i<j<k<ℓ or k<ℓ<i<j. the collection of all secondary structures of the rna sequence a is denoted 𝕊𝕊, and the free energy of s is denoted by e, or simply by e provided that the sequence a is clear from context. the boltzmann probability p=p a for structure s of a is defined by exp/rt)/z, where the partition function z=z=∑s∈𝕊𝕊exp/rt). given two secondary structures s,t of a, the base pair distance d bp between s and t is defined to be the size of the symmetric difference of s,t, i.e. |s−t|+|t−s|. in, rodrigo et al. define intrinsic distance d0=∑s,tp·p·dbp i.e. intrinsic distance is another name for ensemble diversity earlier defined in, and computed by vienna rna package. plasticity is defined in to be normalized ensemble diversity; i.e. p=d0n/ obtained by dividing ensemble diversity by the maximum possible number n/ of base pairs in a structure of a. given two rna sequences a=a …,a n and b=b …,b n of the same length n, rodrigo et al. define d to be the expected base pair distance between structures of a and structures of b minus the ensemble diversity of a, i.e. d1=∑s∈𝕊𝕊∑t∈𝕊𝕊pa·pb·dbp−d since d is not symmetric, this measure is not a metric. in contrast, ensemble distance as described in is a valid metric, defined by the following: dv=∑s∈𝕊𝕊∑t∈𝕊𝕊pa·pb·dbp−d0+d02=∑i<j−pi,j) in, rodrigo et al. define the mutational robustness rm=1−〈d1〉n/ where 〈d 1〉 denotes the average value of d taken over all single point mutants a’ of a. since d is not a true metric, we replace it by the metric d v in our computation of mutational robustness. clearly both notions are closely related. implementation in, mccaskill described a cubic time algorithm to compute the partition function z=z=∑s∈𝕊𝕊exp/rt) for an rna sequence a=a …,a n, where the sum is taken over all secondary structures 𝕊𝕊 of a, e denotes the free energy for the structure s of a with respect to the turner energy parameters, r denotes the universal gas constant and t is absolute temperature. subsequently ding and lawrence described how to use the partition function together with a simple backtracking strategy to sample secondary structures of a from the boltzmann ensemble of low energy structures. if s is a given secondary structure of length n, we define the dual partition function z∗=z∗=∑a∈𝔸𝔸exp/rt) where the sum is taken over all rna sequences a=a …,a n of length n that are compatible with structure s i.e. a i,a j constitute a watson-crick or wobble pair for each base pair ∈s the set of all rna sequences that are compatible with s is denoted by 𝔸𝔸. note that if a sequence a is not compatible with the target structure s then the energy e is infinite, so the corresponding boltzmann factor exp/rt) is zero and the sum in eq. could have been written over all sequences of the same length as s here we describe the efficient software rnadualpf to compute the dual partition function z ∗ and to sample from the low energy ensemble of sequences that are compatible with a given secondary structure s dual partition function if s is a secondary structure on sequence a=a …,a n, then the length of s, denoted by ℓ, is equal to n, while the size of s, denoted by |s|, is the number of base pairs belonging to s. similarly, if secondary structure s is restricted to the interval, where 1≤i≤j≤n, then the length of the restriction of s to, denoted by ℓ, is equal to j−i+ while the size of the restriction of s to, denoted by |s|, is the number of base pairs of s that satisfy i≤x<y≤j. given an rna sequence a=a …,a n, the mccaskill algorithm computes the partition function z defined in eq.. when a is clear from context, z is usually denoted by z. given a target secondary structure s we describe below an algorithm to compute the dual partition function z ∗, defined as the sum of all boltzmann factors exp), where the sum is taken over all rna sequences a∈𝔸𝔸. unlike the mccaskill algorithm, which requires time that is cubic in the length of a, the algorithm presented below requires time that is linear in the length of s our algorithm is motivated by the initialization step of the algorithm info-rna, in which a sequence is determined, for which the free energy with respect to target structure s is a minimum – i.e. info-rna determines argmina e. the algorithm specification requires the notation z ∗, which denotes the sum z∗=∑a,ai=x,aj=yexp−ea,s0/rt of boltzmann factors for sequences a=a i,…,a j for which a i=x,a j=y, and for the restriction s defined by s0=∈s0:i≤x<y≤j. the function z ∗ is defined for all base pairs ∈s 0; these values will be stored in an array, whose rows index base pairs of s and whose columns are indexed by the six canonical base pairs gc, cg, au, ua, gu, ug. once z ∗ has been computed for all base pairs that are visible, i.e. for which there is no base pair for which x<i<j<y, we can compute the full partition function z ∗. au cg gc ua gu ug z s 0 the first column indicates the base pair index which dictates the order in which the dual partition function is computed for different loops closed by the base pair, where we the index of base pair is defined to be the rank of in the total ordering defined in eq.. columns i and j indicate the opening and closing positions of each base pair. type indicates the type of element in the secondary structure closed by each base pair, where r. bulge stands for right bulge, stack for stacking base pair, and int. loop for interior loop. the dual partition function z ∗ of the substructure closed by base pair appears in the rightmost column, while the partition function z ∗ for each of the six canonical base pairs is given in columns 5- note that for base pair sequence constraints depicted in fig. force i and j to be instantiated respectively to g and c, hence the dual partition function z ∗ is zero for any base pair different than gc. the last column of the last row of the table shows the total dual partition function z ∗ for the target structure s 0 following, we define a total ordering on base pairs belonging to the target structure s that satisfy the following precedence rule for any two base pairs,. ≺⇔x<i<j<yori<j<x<y from this ordering, we assign a base pair index to each base pair, which is defined to be the rank of in the total ordering. the following definitions correspond to the turner nearest neighbor energy model, which is an additive loop model where a loop closed by external base pair is designated as a k-loop, if the loop contains k base pairs interior to. therefore, hairpin loops are 0-loops; base pair stacks, bulge loops and internal loops are 1-loops; and multiloops are k-loops for k≥ -way junctions), where the additional count is due to the outer component adjacent to. since au-base pairs that close a loop are energetically unfavorable, in the turner energy model, there is an au-penalty we now define: eau= ifis the outermost pair in a stem ofs having au,ua,gu,ug0otherwise. this au-penalty is applied only if is a base pair adjacent to a triloop, a bulge, an internal loop or a multiloop, or if it is the outermost base pair of an external loop in target structure s and is instantiated by one of the pairs au, ua, gu, ug. when base-paired positions i,j are clear from the context, we write e au. here, we assume that in parsing the input target structure, a list bpcloseelorml has been created of those base pairs, which close either an external loop or a multiloop. let i be the indicator function, it follows that if closes an external loop or multiloop, then exp−i·eaurt is the boltzmann factor for a special au-penalty, otherwise this factor equals for clarity in the notation, this factor is denoted by ert). note that this term is distinct from the factor exprt) applied to base pairs adjacent to a triloop, a bulge or an internal loop, which does not depend on the indicator function. hairpins let close a hairpin in s the hairpin free energy term h, arising solely from entropic considerations, is defined by h=hairpineifj−i−1≤30hairpine+ rtlnj−i−130otherwise where hairpine designates the hairpin free energy obtained from table look-up, when j−i−1≤ triloop let triloop x,y denote the collection of special triloops, xabcy, having an energy bonus triloope. z∗=e−eauirt·exp−h+eaurt×43−|triloopx,y|+∑abc∈triloopx,yexp−triloopert tetraloop let tetraloop x,y denote the collection of special tetraloops, xabcdy, having an energy bonus tetraloope. similarly, given nucleotides n n2∈n, tetraloop x,y denotes the collection of special tetraloops of the form xn 1 abn 2 y. define z ∗ by z∗=e−eauirt·exp−hrt×∑n n2∈nexp−mismatchrt×42−|tetraloopx,y|+∑ab∈tetraloopx,yexp−tetraloopert hexaloop let hexaloop x,y denote the collection of special hexaloops, xabcdefy, having an energy bonus hexaloope. similarly, given nucleotides n n hexaloop x,y denotes the collection of special hexaloops of the form xn 1 abcdn 2 y. define z ∗ by z∗=e−eauirt·exp−hrt×∑n n2∈nexp−mismatchx,y,n n2rt×44−|hexaloopx,yn n2|+∑ab∈hexaloopx,yn n2exp−hexaloopexn1abcdn2yrt hairpin size exceeds four and is different than six define z ∗ by z∗=e−eauirt·exp−hrt×∑n n2∈nexp−mismatchx,y,n n2rt·4j−i− stacked base pairs, bulges and internal loops here, we consider the case of a 1-loop, which comprises the case of stacked base pairs, bulges and internal loops. the following cases correspond to each possibility. stacked base pair in this case, stacks on the base pair, and the partition function z ∗ has been computed. let stack denote the free energy of base stack 5′−xu−3′3′−yv−5′ obtained by table look-up. z∗=e−eauirt·∑uv∈ℬexp−stackrt×z∗i+ j− u,v bulge loop in this case, closes a bulge in s since bulge size may exceed the values in table look-up, we define the free energy for a bulge of size r by bulge=bulgeeifr≤30bulgee+ rtlnr30otherwise. if closes a left bulge of size r in s then the bulge is closed by base pair involving nucleotide pair u,v, and z∗=e−eauirt·∑uv∈ℬexp−eaui,j,x,yrt×exp−bulgert·4r·z∗i+r+ j− u,v while if closes a right bulge in s then the bulge is closed by base pair involving nucleotide pair u,v, and z∗=e−eauirt·∑uv∈ℬexp−eaurt×exp−bulgert·4r·z∗i+ j−r− u,v internal loop in this case, closes an internal loop in s whose left portion is of size r [resp. r 2]. since internal loop size r=r 1+r may exceed the values in table look-up, we define the free energy for an internal loop of size r by internal=internaleifr≤30internale+ rtlnr30otherwise. the closing base pair of the internal loop of size r=r 1+r may involve the nucleotides uv∈ℬ, while the unpaired nucleotides in positions i+ j− i+r j−r may involve a,b,c,d∈n. in addition, there is an energy penalty for non symmetric internal loops, min, where the value of the constants asym and maxasym are given in the turner energy model. thus z∗=e−eauirt·exp−minasym·|r1−r2|,maxasymrt×∑uv∈ℬ∑a,b,c,d∈nexp−eaui,j,x,yrt×exp−internalr1+r2rt·4r1+r2−4×exp−mismatch+mismatchrt ×z∗i+r1+ j−r2− u,v external loop despite the fact that, by following the total order on base pairs defined in eq., the dual partition function of multiloops is always computed before the dual partition function of the external loop, the computation of the dual partition function of multiloops will be easier to understand if the dual partition function of the external loop is defined in advance. in order to improve speed, some implementations of rna thermodynamics-based algorithms ignore the contribution of dangling positions, which corresponds to vienna rna package -d flag. rnadualpf also includes this option, which dramatically increases the speed of the algorithm. the reason behind this difference of performance is clear from the following definitions. suppose that h= constitutes the list of k external base pairs of s where i 1<j 1<i 2<j 2< ⋯ < i k < j k. for each, with 1≤r≤k, and for each choice of base pair gc, cg, au, ua, gu, ug, the value z ∗ has been previously computed and stored by dynamic programming, as well as the sum z ∗. when the contribution of dangles is ignored, the dual partition function of an external loop with ℓ nucleotide positions external to every base pair is defined by z∗=4ℓ·∏r=1kz∗ where ℓ=n−∑r= …,k and n is the length of the target structure s the default treatment of dangles in rnadualpf described below corresponds to vienna rna package -d flag, where both flanking positions of each external base pair contribute to the free energy. let d=[a b …,a k,b k]⊆[i 1− j 1+ ⋯,i k,j k] be a list of those nucleotide positions that are adjacent to the k external base pairs,…,. the ordered multiset [a b …,a k,b k] can be considered as a collection of constraints, so that if a 2=i 2− and a 2=j 1+ then a 2=b and any nucleotide value that is assigned to b must simultaneously be assigned to a moreover, there can also be an overlap between the list of base paired positions in h [i j …,i k,j k] and the multiset d=[a b …,a k,b k]. if j 1=i 2− then b 1=i and a 2=j therefore, in the computation we have to account for these constraints. let m denote the number of unpaired positions in d, without repetitions, and define a r,b r as the nucleotides instantiated respectively at a r,b r. the energy term for a 5′-dangle on base pair with nucleotides u,v is denoted by e d [resp e d3] where the dangle position x− is assigned nucleotide w. with the notation just described, we have z∗=∑〈u v …,uk,vk〉∈ℬk∑a b …,ak,bk∈n2k4ℓ−m×∏r=1kz∗ir,jr;ur,vr×exp−ed5ir,jr,ar;ur,vr,ar+ed3ir,jr,br;ur,vr,brrt depending on the target structure s it can happen that the second sum of eq. must be restricted to range over strictly less than 42k many rna sequences. this is explained as follows. if i 1= [resp. j r=n] then there is no position for a 5′ dangle, and hence the nucleotide sequences considered in the second summation would have length strictly less than 2k. moreover, certain 5′ dangled positions could be identical to 3′ dangle positions, which arises for instance when j k+2=i k+1; alternatively, certain dangled positions could be identical with base-paired positions, which arises for instance when j k+1=i k+ in such situations, instantiations of the 3′-dangle on and the 5′-dangle on are not independent, thus leading to a restriction of the range of the second summation in eq.. a similar restriction is implicitly assumed in the treatment of external loops in this section and of multiloops in the next section. the algorithm performance can be improved by dividing the external loop into groups of components having interdependently constrained dangling positions, as just explained. define two base pairs, as adjacent if x<y<x ′<y ′ and x ′−y≤ – i.e. dangling positions of the base pairs, are constrained. let g denote a maximal collection of adjacent base pairs belonging to h=, together with their associated dangle positions in d=[i 1− j 1+ …,i k− j k+1]. it is important to note that h∪d is thus partitioned into a collection of g disjoint groups g=. therefore, we can divide an external loop of k helices into a collection groups g of size g≤k, and p unpaired positions that are external to every base pair of s and not adjacent to any base pair. for a group g with h base pairs, let h= denote the list of base pairs in g, and let d=[α β …,α h,β h]⊆[κ 1− λ 1+ ⋯,κ h− λ h+1] denote their associated dangle positions. if u r,v r,a r,b r denote the nucleotides instantiated at the base pair r= and its respective dangling positions α r,β r respectively, then the dual partition function of g is the following. z∗=∑〈,…,〉∈ℬh∑a b …,ah,bh∈n2h×∏r=1hz∗κr,λr;ur,vr×exp−ed5κr,λr,αr;ur,vr,ar+ed3κr,λr,βr;ur,vrbrrt where the range of the second summation can be constrained by the overlap among positions in d and between positions in d and h, as explained for eq.. finally, since there are no shared dangling positions between groups, the dual partition function of an external loop is defined by z∗=4p·∏r=1gz∗. multiloop suppose that closes a multiloop in s which is a k-loop, or -way junction, for k> where there are ℓ unpaired bases in the multiloop. suppose that the k components of the multiloop are closed by the base pairs,…, with the property that i<i 1<j 1<i 2<j 2<⋯<i k<j k<j. assume that for all nucleotide choices in ℬ for each of the k base pairs of the multiloop, for 1≤r≤k, the value z ∗ has previously been computed and stored by dynamic programming, as well as the sum z ∗. the computation of the dual partition function is similar to that of the external loop. however, in this case we have to add the contribution of the base pair closing the multiloop, the au-penalties applied to this base pair, and the energetic penalty of a multiloop a+b·+c·ℓ, where the values of the constants a, b and c are given in the turner energy model. then, the dual partition function of a multiloop without accounting for dangling positions is z∗=e−eauirt·exp−a+b·+cℓrt·4ℓ×exp−eaurt·∑〈u v …,uk,vk〉∈ℬk×∏r=1kz∗ir,jr;ur,vr the notation we use to define the dual partition function of multiloops with dangling positions is similar to that described for external loops. however, some modifications are required in the previously given definitions, since we have to take into account the flanking positions of the base pair closing the multiloop. let h= be the collection of k base pairs closing one of the k components of the multiloop, and the base pair closing the multiloop, and define the multiset d=[a b …,a k+ b k+1]⊆[i 1− j 1+ ⋯,i k− j k+ i+ j−1] of nucleotide positions adjacent to the base pairs in h. due to the possible overlap with the base pair closing the multiloop and its flanking positions, there are additional constraints in the ordered multiset [a b …,a k+ b k+1], so that if a 1=i 1− and i 1=i+ then a 1=a k+ and any nucleotide value that is assigned to a must simultaneously be assigned to a k+ moreover, there can also be an overlap between the list of base paired positions [i j …,i k,j k,i,j] and the multiset [a b …,a k+ b k+1]. if i=i 1− then a k+1=i and a 1=i. let m denote the number of unpaired positions in d, without repetitions. then, the dual partition function of a multiloop with dangling positions is defined as follows. z∗=e−eauirt·∑〈u v …,uk,vk〉∈ℬk∑a b …,ak+ bk+1∈n2exp−a+b·+cℓrt·4ℓ−m·exp−eaui,j,x,yrt×∏r=1kz∗ir,jr;ur,vr×exprt×exp−ed3j,i,ak+1;y,x,ak+1+ed5j,i,bk+1;y,x,bk+1rt as explained for eq., it can happen that the second summation must be restricted to range over strictly less than 42k many rna sequences. a decomposition similar that for external loops can be performed to improve the performance in the computation of the dual partition function of a multiloop. in a multiloop, in addition to the adjacency definition given for external loops, we consider the base pair that closes the multiloop as adjacent to a base pair that closes a component of the multiloop, where i<x<y<j, if either x≤i+ or y≥j− then, let g denote a maximal collection of adjacent base pairs belonging to h=, together with their associated dangle positions in d=[i 1− j 1+ …,i k− j k+ i+ j−1]. this decomposition produces a collection g of g disjoint groups g …,g g, one of which, designated the closing group g c contains the closing base pair of the multiloop, and g− of which, designated as non-closing groups g nc, do not contain the base pair. non-closing groups have the same composition as those defined for external loops – i.e. a collection of h base pairs h= and a set of dangling positions d=[α β …,α h,β h]⊆[κ 1− λ 1+ ⋯,κ h− λ h+1]. therefore, we can compute the dual partition function z of a non-closing group as described in eq.. in addition, the collection of non-closing groups of size g− of a multiloop of k components is denoted by gnc, where 0≤≤k. therefore, a multiloop of k components and ℓ unpaired positions can be decomposed into one closing group g c, a collection of non-closing groups gnc, and p unpaired positions that are not adjacent to any base pair, with 0≤p≤ℓ. in a non-closing group, the collection of base pairs of size h+ is denoted by h=, where the base pair closing the multiloop is at the last position. the ordered multiset of adjacent positions is denoted by d=[α β …,α h+ β h+1]⊆[κ 1− λ 1+ ⋯,κ h− λ h+ i+ j−1], where the positions adjacent to i and j are at the last positions are respectively denoted by α h+ β h+ a graphical example of a closing group and a non-closing group is shown in fig. 1 e, where the positions of a non-closing group with base pair are highlighted in green and the positions of the closing group are highlighted in red and blue, and where the base pair that closes the multiloop is depicted in red. fig. sampling dependency examples in rnadualpf for different structural elements: stacked base pair, hairpin, 1× internal loop, 3× internal loop and multiloop. base pair to be sampled is highlighted in red, positions whose energy contribution is dependent on the instantiation of are highlighted in blue, and positions that are mutually dependent, but independent of the instantiation of, are highlighted in green. unpaired positions where the nucleotide choice has no effect in the free energy of the structure are indicated in black for a closing group g c with h+ base pairs in h= and their flanking positions d=[α β …,α h+ β h+1]⊆[κ 1− λ 1+ ⋯,κ h− λ h+ i+ j−1], let x,y denote the nucleotides assigned to the closing base pair of the multiloop, and let u r,v r,a r,b r denote the nucleotides assigned respectively to the base pair r= and its flanking positions α r,β r. then, the the dual partition function z ∗ of the closing group is defined by e−eauirt·∑〈u v …,uk,vh〉∈ℬh∑a b …,ah+ bh+1∈n2×exp−eaui,j,x,yrt·∏r=1hz∗κr,λr;ur,vr×exp−ed5κr,λr,αr;ur,vr,ar+ed3κr,λr,βr;ur,vr,brrt×exp−ed3j,i,αh+1;y,x,ah+1+ed5j,i,βh+1;y,x,bh+1rt in the same way as in eq., the values of the second summation are constrained to the possible choices among overlapping positions. then, the dual partition function z ∗ of the multiloop with k components and ℓ unpaired positions, where p of which are not adjacent to any base pair, is defined by z∗=exp−a+b·+cℓrt·4p×z∗gc;x,y·∏gnc∈gncz∗gnc sampling once the dual partition function z ∗ and its subcases z ∗ for each base pair have been computed, it is possible to perform a boltzmann weighted sampling of positions i and j. for example, given the target structure with sequence constraints depicted in fig. rnadualpf computes the dual partition function table shown in table the dual partition function of the substructure enclosed by the base pair is z ∗, and the dual partition function of the substructure enclosed by the base pair where i,j are currently instantiated by the nucleotides x,y is denoted by is z ∗. therefore, the boltzmann probability of x,y at positions i,j in the substructure enclosed by the base pair is z ∗/z ∗ and can be sampled using the roulette wheel method. fig. target structure with sequence constraints used as input of rnadualpf to compute the dual partition function values shown in table sequence constraints are highlighted in red due to the turner energy model, it is necessary to determine nucleotide positions whose instantiation influences the energy of other positions, and subsequently all mutually dependent positions must be instantiated simultaneously. figure illustrates the mutual dependencies that must be considered when sampling different types of elements, where the base pair to be sampled is highlighted in red, positions whose sampling probability is dependent on the instantiation of are highlighted in blue, and positions that are mutually dependent, but independent of the instantiation of, are highlighted in green. since the dynamic programming algorithm for the dual partition function proceeds from inner to outer base pairs, using the total ordering ≺ in eq., the sampling order of base pairs proceeds from outer to inner positions, i.e. from largest base pair index to smallest. in order to account for mutual dependencies in the sampling step, we define the function sample for each base pair in s where k indicates the base pair index defined from eq., t indicates the type of structural element closed by base pair in the target rna secondary structure, as shown in table and x,y are the instantiated nucleotides at positions. due to the mutual dependencies, sampling a base pair with base pair indexk closing an m-loop, for m> forces the instantiation of all inner closing base pairs of the m-loop, and the base pair index of each such inner base pair is strictly less than k. for this reason, except in the case of external loops, the outermost base pair has been always instantiated before sample is called, and therefore the instantiation x,y is given as a parameter of the sampling function. the boltzmann probability of each possible instantiation of mutually dependent positions can be computed on the fly in the backward step. however, in order to improve the speed of the algorithm, in the forward step rnadualpf stores the conditional dual partition function values of instantiations of interdependent positions. these tables are used by the sampling function, since each value corresponds to the dual partition function conditional on a specific instantiation of the positions to be sampled by sample. since the sampling procedure depends on the type t of element, we describe the function sample for each type of element – hairpin, stacked base pair, internal loop, multiloop and external loop as depicted in fig. for each of these cases, the values are stored in the conditional dual partition function table associated with the closing base pair. hairpins when hairpin size exceeds three, since the base pair has been previously instantiated, flanking positions i+ j− are sampled first. given the current assignment x,y, the boltzmann probability of sampling respectively the nucleotides u,v at the flanking positions i+ j− is pi+1=u,j−1=v|i=x,j=y=z∗i,j,i+ j−1;x,y,u,vz∗i,j;x,y therefore, in the forward step rnadualpf stores in a table the conditional dual partition function of each possible instantiation {x,y,u,v} of the base pair and its flanking positions i+ j− respectively, defined by z∗i,j,i+ j−1;x,y,u,v=e−eauirt·exp−hrt·exp−mismatchrt·4j−i− then, remaining unpaired positions are uniformly sampled, since the nucleotide choice does not change the final free energy. triloops, tetraloops and hexaloops are exceptions to this rule, since there are special loops that contribute to or penalize the free energy. in those cases, we have to account for the special loops, as defined in “hairpins” section. although it could seem to be a waste of space to store a different conditional dual partition function table for each base pair, even for two different hairpins of the same size in the target structure, one should note that rnadualpf allows sequence constraints, and thus z ∗ could possibly differ from z ∗ when and close hairpins of the same size. stacking base pairs as depicted in fig. 1 b, sampling probability of a base pair with base pair index k− is dependent on the value sampled at the adjacent stacking base pair with base pair indexk. therefore, sample samples the base pair using the conditional probability given the instantiation of base pair by x,y, defined as follows: pi+1=u,j−1=v|i=x,j=y=z∗i,j,i+ j−1;x,y,u,vz∗i,j;x,y the conditional dual partition function values stored in the forward step correspond to each instantiation {x,y,u,v} of the base pairs,, denoted by z∗i,j,i+ j−1;x,y,u,v=e−eauirt·exp−stackrt×z∗i+ j− u,v internal loops the energy contribution of internal loops in the turner energy model depends on the flanking unpaired positions of both the inner and outer closing base pairs, hence the sampling probability of the inner base pair cannot be separated from the adjacent unpaired positions. moreover, for specific sizes of internal loop, the inner and outer closing base pairs share flanking positions. in these cases, all the unpaired positions and the outer base pair must be sampled at the same time, since the energy contribution of each combination of base pairs and flanking positions is different. in the 1× internal loop depicted in fig. 1 c, if the outer base pair is instantiated by x,y, then let p denote the probability of sampling the nucleotides u,v,a,b,c respectively at positions k,l,n n n where is the inner closing base pair, n is the flanking position at i+ shared by the base paired positions i and k, and n and n in the following equation, let p abbreviate the conditional probability just defined. then pu,v,a,b,c|x,y=z∗i,j,k,l,n n n3;x,y,u,v,a,b,cz∗i,j;x,y rnadualpf computes and stores the conditional dual partition function of each possible instantiation {x,y,u,v,a,b,c} respectively at positions i,j,k,l,n n n where the value z ∗ is defined by e−eauirt·exp−min−|,maxasym)rt×4j−l−3·exp−eaurt×exp−internalk−i+j−l−2rt×exp−mismatch+mismatchrt×z∗ for internal loops of sizes similar conditional dual partition function tables are computed following the definitions in “stacked base pairs, bulges and internal loops” section. other internal loops: when there are no shared flanking positions between the two base pairs that close an internal loop, as depicted in fig. 1 d, the energy contribution of innermost base pair and its respective flanking positions is independent of those of the outermost base pair. in this case, rnadualpf samples first the flanking positions i+ j− of the outermost base pair, whose sampling probability is solely dependent on the instantiated nucleotides x,y at positions i,j. is not necessary to store any conditional dual partition function for sampling these positions, since the probability of sampling the values a,b at the flanking positions i+ j− given the assignment x,y is defined by pi+1=a,j−1=b|i=x,j=y=exp−mismatchrt∑c,d∈nexp−mismatchrt where mismatch penalties are obtained from table look-up. finally, the innermost base pair and its flanking positions k− l+ are sampled together. in this case, we need to store an additional value z ∗, which is given by z∗k− l+1=∑uv∈ℬ∑c,d∈nexp−mismatchrt·z∗ then, following the same notation, the probability of sampling the nucleotides v,u,d,c respectively at positions k,l,k− l+ is pk=v,l=u,k−1=d,l+1=c=z∗k,l,k− l+1;v,u,d,cz∗k− l+ therefore, the conditional dual partition function of each possible instantiation {v,u,d,c} stored in the corresponding table is defined as z∗k,l,k− l+1;v,u,d,c=exp−mismatchrt·z∗ finally, since the remaining unpaired position does not contribute to the free energy, it is uniformly sampled. multiloops and external loops as explained in “external loop” section, if dangling positions are not included in the computation, sampling an external base pair or the closing base pair of a multiloop from z ∗ is trivial. on the other hand, by including dangling positions in the sampling, there is a dramatic increase in the space complexity of rnadualpf, albeit the space used is only a constant factor larger. however, the decompositions into groups described in “external loop” and “external loop” sections allow to sample the positions of each group independently. the example shown in fig. 1 e depicts a multiloop with two groups: a non-closing group g nc highlighted in green, and a closing group g c highlighted in red and blue, where the closing base pair of the multiloop is marked in red. in a non-closing group g nc all base pairs in h and dangling positions in d must be sampled together. therefore, the conditional dual partition function of each possible instantiation of nucleotides at the h closing pairs in h and their adjacent positions in d is stored. let u={u v …,uh,vh} denote an instantiation of the h base pairs in h=[κ λ …,κ h,λ h], and let w={a b …,ah,bh} denote an instantiation of the h flanking positions in d=[α β …,α h,β h] in the non-closing group g nc. then, the probability of sampling u,w is p=u,d=w)=z∗g,hgnc,dgnc;u,wz∗ therefore, the conditional dual partition function of each instantiation u,w at h,d, stored in the table of the group, is defined by z∗g,hgnc,dgnc;u,w=∏r=1hz∗κr,λr;ur,vr×exp−ed5κr,λr,αr;ur,vr,ar+ed3κr,λr,βr;ur,vr,brrt recall that the base pairs in h are adjacent. therefore, due the constraints given by the overlapping positions within d, and between d and h, explained in “external loop” section, the number of possible instantiations u,w of h,d is ≤. in a similar way, sampling from the closing group g c closed by the base pair, with h+ base pairs in h and their corresponding flanking positions in d requires us to store the conditional dual partition function of each instantiation of nucleotides {x,y,u,w} respectively at i,j,h,d, where u={u v …,uh,vh} denotes an instantiation of the h first base pairs in h, w={a b …,ah+ bh+1} denotes an instantiation of the 2· flanking positions in d=[α β …,α h+ β h+1], and x,y denotes an instantiation of. the probability of the instantiation u,w, given the nucleotides x,y is ph=u,d=w|i=x,j=y=z∗gc,i,j,h,d;x,y,u,wz∗gc;x,y then, the values stored in the table of the closing group correspond to the conditional dual partition function of each instantiation x,y,u,w are given by z∗,d;x,y,u,w), which is defined by the following expression: e−eauirt·exp−eaurt·∏r=1hz∗κr,λr;ur,vr·exp−ed5κr,λr,αr;ur,vr,ar+ed3κr,λr,βr;ur,vr,brrt·exp−ed3j,i,αh+1;y,x,ah+1+ed5j,i,βh+1;y,x,bh+1rt as a final remark, we would like to recall that all the conditional dual partition function values are computed and stored in the forward step at the same time as the dual partition function. therefore, despite the consequent increase of space complexity in the algorithm, the computation of the values required for correct sampling does not involve a greater time complexity. scaling the sequence partition function z ∗ grows much faster than the usual structure partition function z, and so scaling must be used in the implementation. let c> be a user-defined constant. by a slight modification of the previous recursions, we actually compute z†=z∗cj−i+ and hence z†=z∗cn, where n is the length of s for instance, the analogue of eq. is z†=z∗i,j;x,ycj−i+1=e−eauirt·exp−hrtcj−i+1×∑n n2∈nexp−mismatchrt·4j−i− and the analogue of eq. is z†=z∗cj−i+1=e−eauirt·12·∑uv∈ℬexp−stackrt×z†i+ j− u,v this modification does not affect properties of sequences sampled from the low energy ensemble, since the same scaling factor appears in both the numerator and denominator of all conditional probabilities. for instance, the analogue of eq. is pi+1=u,j−1=v|i=x,j=y=z∗i,j,i+ j−1;x,y,u,vz∗i,j;x,y=z†i,j,i+ j−1;x,y,u,vz† controlling gc-content the gc-content of an rna sequence a=s …,s n is the number of nucleotides that are either g or c. instead of computing z ∗ and z ∗, we can compute z ∗ and z ∗, defined to be the corresponding partition dual partition functions, restricted to sequences having gc-content of α. note well that gc-content α includes the closing nucleotides x and y respectively located at positions i and j; i.e. z∗i,j;x,y;α=∑ai,…,aj,gc=αai=x,aj=y,ai+ …,aj−1∈nexp−eai,…,aj;s0/rt where s denotes the restriction of target structure s to the interval. we describe two particular subcases, to provide the idea of how modifications need to be undertaken. triloop note that the number of rna sequences of length m having gc-content of α is mα·2α·2m−α=mα·2m≤4m, since α selected positions must be either g or c, yielding the term 2α, while the remaining m−α positions must be either a or u, yielding the term 2m−α. assume that γ=|{x,y}∩{g,c}|=β. then z∗i,j;x,y;α=e−eauirt·exp−h+eaurt×j−i−1α−β·2j−i−1−|triloopx,y|+∑abc∈triloopx,yγ=α−βexp−triloopert multiloop and external loop assume that closes a multiloop, which is a -way junction with ℓ unpaired nucleotides. assume that the ordered multiset of potential dangle positions is d=[a b …,a k+ b k+1], where a r=i r− and b r=j r+ for r= …,k, and a k+1=i and b k+1=j, and assume that there are m unpaired positions that are not adjacent to a base pair in the multiloop. if r denotes an rna sequence of arbitrary length, then let the function γ denote the gc-count in r. given an assignment of nucleotide base pairs u 1 v …,u k v k to,…,, where u r v r∈{gc,cg,au,ua,gu,ug}, and given an assignment a b …,a k,b k of dangle nucleotides, where ar,br∈n, for r= …,k, we let γ=γa …,ak,b …,bk. then the dual partition function of a multiloop with a gc-content of α is defined by setting z ∗ equal to the following: e−eauirt·∑α1+⋯+αk≤α∑{ur,vr∈ℬ:r= …,k}∑{a b …,ak,bk∈n2k}exp−a+b·+cℓrt·α−∑r=1kαr−γ·2ℓ×exp−eaurt·∏r=1kz∗ir,jr;ur,vr;αr×exp−ed5ir,jr,ar;ur,vr,ar+ed3ir,jr,br;ur,vr,brrt×exp−ed3j,i,ak+1;y,x,ak+1+ed5j,i,bk+1;y,x,bk+1rt since the modification required in the remaining cases follows similar reasoning as in the treatment of the hairpin and external loop just described, the details for these remaining cases are not given. an additional challenge of computing the dual partition function with gc-content control is the combinatorial problem of efficiently counting the number n of instantiations of the external loop, consisting of all positions external to every base pair, with gc-content k, where the user can stipulate that certain positions are constrained to contain nucleotides consistent with iupac codes. to this end, we implemented the combinatorial algorithm defined in supplementary information. sampling with gc-content the implementation of sampling with gc-content is performed in a similar manner as described in “sampling” section, with some notable differences. first, the sampling function is redefined by sample, where k indicates the base pair index in the ordering defined by eq. for the base pair that is already instantiated by nucleotide pair xy, and t designates the type of structural element closed by base pair in the target rna secondary structure, as shown in table the function sample instantiates all positions of the loop having outer closing base pair, including its inner closing base pair and which returns the gc-content of the sampled loop. moreover, the gc-content of the subsequence a= will be α once the entire sequence a …,a n is sampled. second, rnadualpf stores a conditional dual partition function table for each base pair and gc-content to j-i- the function sample samples from the conditional dual partition function of those sequences which have exactly α gs and cs strictly between the positions i and j, thus guaranteeing a gc-content of α for the subsequence a once the entire sequence a …,a n is sampled. note that sample samples only the loop closed by the already instantiated outer base pair, and that α is the gc-content of the entire subsequence a=a i+ …,a j− once the algorithm terminates. only in the case that base pair closes a hairpin loop will it generally happen that the gc-content of the loop closed by is equal to α. let α be the user-designated gc-content of sequences a=a …,a n to be sampled from a target secondary structure having ℓ base pairs. the following pseudocode describes how to sample sequences a=a …,a n, whose gc-content exactly equals α. here, an external loop with m components means that there are m exterior base pairs,…, such that all positions exterior to these base pairs are unpaired; i.e. each position r∈{ …,n}−∪c=1m{ic,…,jc} is unpaired in the target structure. to clarify how the gc-content is sampled in a statistically rigorous manner, suppose that the user has specified the gc-content to be α, and that l is the external loop of the target structure s having m components, where the cth component has external closing base pair. in computing the dual partition function, for all possible choices of non-negative integers α …,α m,β that sum to α and all 6m possible assignments of watson-crick or wobble nucleotide pairs x y …,x m,y m to the base pairs,…,, the software rnadualpf has computed the sum of ∑c=1mzc∗ plus the boltzmann factor of the external loop with gc-content β. since the dual partition function z ∗ is the sum, taken over all values of α …,α m,β and all watson-crick and wobble pair assignments to the external base pairs, rnadualpf can then use the roulette wheel method to sample values α …,α m,β and x y …,x m,y m in a statistically rigorous manner. multiloops, and other structural elements, which contain unpaired regions whose sequence does not contribute to the free energy of the structure, are handled in a analogous manner.  robustness and plasticity of c. elegans mirnas and e. coli sncrnas in borenstein and ruppin used version of the vienna rna package to generate rna sequences per wild type precursor microrna extracted from the database rfam, with the property that each of the control sequences folded into the wild type pre-mirna structure – i.e. the minimum free energy structure of each of the control sequences was identical to the mfe structure of the wild type pre-mirna. based on these computational experiments, borenstein and ruppin asserted that the “structure of mirna precursor stem–loops exhibits a significantly high level of mutational robustness in comparison with random rna sequences with similar stem–loop structures”. noting that the vienna rna package inverse folding program rnainverse does not control for gc-content or other sequence compositional bias, the authors performed a second computational experiment, in which control sequences not only folded into the target wild type structure, but also had similar dinucleotide composition to that of wild type pre-mirna. since the filtering step required enormous run time and computational resources, the authors restricted their attention to a small set of micrornas, and generated only control sequences per microrna – note here that rnainverse cannot control for gc-content. borenstein and ruppin concluded that robustness of precursor micrornas was not the byproduct of a base composition bias or of thermodynamic stability. subsequently rodrigo et al. undertook a similar analysis for bacterial small noncoding rnas, also using the program rnainverse, albeit using somewhat different definitions – see precise definitions in “formal definitions of robustness” section. the main finding of was that bacterial sncrnas are not significantly robust when compared with sequences having the same structure, as computed by rnainverse; however, the authors found that bacterial sncrnas tend to be significantly plastic, in the sense that the ensemble of low energy structures are structurally diverse. unlike the case of precursor micrornas, rodrigo et al. did not control for sequence compositional bias. using rnadualpf, we performed similar computational experiments on precursor micrornas of c. elegans from mirbase and for the bacterial small noncoding rnas of. below, we discuss each case separately. for each c. elegans pre-mirna, we used rnadualpf to sample sequences with no control over gc-content and sequences whose gc-content was identical with that of the wild type pre-mirna. moreover, each control sequence approximately folded into the mfe wild type pre-mirna structure as computed by vienna rna package.. table shows that length-normalized base pair distance between the mfe structure of the control sequence and that of the pre-mirna is on average ± for default use of rnadualpf with control over gc-content, and ± when gc-content of each control sequence is identical to that of the corresponding wild type pre-mirna. additional measures in table show that sequences sampled from rnadualpf are only modestly more stable thermodynamically, the ensemble of low energy structures of control sequences deviate slightly more from the target pre-mirna structure, as is the case for wild type pre-mirna sequences, as mesured by ensemble defect, expected base pair distance to target, expected proportion of native contacts, average positional entropy, morgan-higgs structural diversity, and vienna structural diversity. for each of the wild type pre-mirna sequences, rnadualpf sampled sequences, either having exactly the same gc-content as the wt sequence or with no control over gc-content. the designation ‘mfe’ indicates that the sampled sequences were subsequently filtered to retain only those, whose minimum free energy structure is identical to the mfe structure of the corresponding wt pre-mirna; otherwise, the designation ‘def’ is used to indicate the default output of rnadualpf, without the subsequent filtering step. for each wt pre-mirna sequence, rnadualpf generated sequences for the default case def, and sequences for the non-default case mfe, such that sample mfe structure is identical to wt mfe structure. various measures were used to compare the properties of rnadualpf sampled sequences to those of wild type sequences: bp dist target: length-normalized average base pair distance d bp between the mfe structure s of sequences sampled by rnadualpf and the target structure s ∗. energy mfe: length-normalized average free energy e of mfe structure s energy target: length-normalized average free energy e of target s ∗ for the respective sequences. ensemble defect: length-normalized expected hamming distance to target s ∗. exp bp dist: length-normalized expected base pair distance to target s ∗. prop nat contact: expected proportion of base pairs of target s ∗ that occur in the mfe structure, i.e. 〈|s0∩s∗||s∗|〉. pos entropy: average positional entropy. gc content: average proportion of positions occupied by g or c. ln dual prob: average natural logarithm of the dual probability exp/rt)/z∗ that sequence a adopts the structure s. ln prob: average natural logarithm of the probability exp/rt)/z that sequence a adopts the structure s. mh str div: length-normalized morgan-higgs structural diversity. vienna str div: length-normalized vienna structural diversity, called ensemble diversity in. values of all measures for default sampled sequences having gc-content within % of wild type gc-content are essential identical to those of exact gc-content control tables and display a similar analysis of the collection of bacterial small noncoding rnas of and of rfam database. for the rfam database, we selected one sequence from each of the ≈ rfam families, with the property that the mfe structure of the sequence most resembled the rfam consensus structure – i.e. whose mfe structure has smallest base pair distance to the consensus structure. these tables show similar trends as those displayed in table although values are larger due to increased sequence length of bacterial sncrna and sequences from rfam. see table for an explanation of column headers and various measures. since bacterial noncoding rna is generally much longer than precursor microrna, no subsequent filtering step was undertaken to ensure that sample sequence mfe structure is identical to that of wild type pre-mirna. however an additional column is given for sequences required by rnadualpf to have gc-content is within % of wt value. for each rna family from rfam, we selected that sequence whose mfe structure had smallest base pair distance to the rfam consensus structure for the family. these sequences constituted the collection wt. see table for an explanation of column headers and various measures in agreement with, the left panel of fig. shows that c. elegans mirna is significantly robust, provided that gc-content is not controlled. however, in contrast to, when gc-content is controlled, we find that c. elegans mirna is significantly non-robust. to corroborate our findings, for each wild type c. elegans pre-mirna, we performed a second computational experiment, to generate sequences with no control over gc-content and sequences whose gc-content was identical with that of the wild type pre-mirna. in contrast to the first experiment, we used rnadualpf to generate sufficiently many sequences to subsequently select sequences and sequences, each of whose mfe structure was identical to that of wild type pre-mirna. the left panel of fig. shows that when gc-content is not controlled, c. elegans precursor micrornas are statistically robust, in agreement with the main result of borenstein and ruppin. however, when gc-content of control sequences is identical to that of wild type precursor microrna, we confirm that c. elegans pre-mirna is statistically non-robust. note that our finding, which is in opposition to results of borenstein and ruppin, is based on a larger data set of precursor micrornas, each of which has a larger control set, than in the analysis of. fig. z-scores of mutational robustness and of plasticity are presented for the bacterial small noncoding rna collection from and for c. elegans precursor microrna from mirbase for each wild type pre-mirna wild type sequence, rnadualpf sampled sequences using the minimum free energy structure of the wild type sequence as target structure. the gc-content of the sampled sequences was either required to be exactly that of the wild type sequence, or not, as indicated in the legend. sampled sequences were used to compute the mutational robustness and plasticity, as explained in the main text. note that c. elegans mirna is significantly robust if gc-content is not controlled, but significantly non-robust if gc-content of rnadualpf samples is identical to that of wild type pre-mirna. similarly, bacterial sncrnas are not significantly robust if gc-content is not controlled, but significantly non-robust when gc-content is identical to that of wild type sncrna. for this figure, mutational robustness of rna sequence a=a …,a n is defined by 1−〈dbp〉n, where ensemble distance d bp between two length n sequences a and b is defined in, and the average ensemble distance from all single-point mutants of a is defined by 〈dbp〉=∑bdbp3n where the sum is taken over all single-point mutants b of a. we use this notation of mutational robustness, rather than the notion defined in, since the latter notion is not a true metric, as explained in “formal definitions of robustness” section. the plasticity p=〈dv〉n/2=∑i<jpi,jn is defined in as normalized ensemble diversity, where ensemble diversity d v is defined by eq. fig. z-scores of mutational robustness and of plasticity are presented for c. elegans microrna from rfam. robustness and plasticity were measured as explained in fig. for each wild type sequence, sequences were sampled from rnadualpf both with and without control over gc-content, while sequences were generated by rnadualpf, both with and without control over gc-content, to exactly fold into the target structure. this was achieved by repeatedly sampling sequences in order to obtain a sequence, whose mfe structure was identical to the target structure. the number of sequences necessary to sample in order to obtain one sequence that folds into the target structure was variable, depending on the target structure and gc-content – in many cases, only samples were necessary per selected sequence, in some cases samples were necessary. and in one specific case samples were required. by control over gc-content, we mean that all sampled sequences have identical gc-content to the wild type sequence turning now to the analysis of bacterial small noncoding rnas, we find that sncrnas are not significantly robust when gc-content is not controlled, confirming a result from rodrigo et al.. however, when gc-content of the sequences sampled from rnadualpf is required to be identical to that of wild type sncrna, bacterial sncrnas are see to be significantly non-robust. note here that rodrigo et al. used rnainverse in their computational experiments, hence could not consider the case with control over gc-content. the left panels of figs. 3 a and summarize our findings that precursor micrornas are significantly non-robust with respect to a control set of sequences generated by rnadualpf with identical gc-content to that of the wild type sequence. finally, in our analysis of plasticity, the right panels of figs. and show that both c. elegans and bacterial small noncoding rnas exhibit more plasticity when compared with control sequences for which gc-content is not controlled, as well as when compared with control sequences for which gc-content is identical to that of wild type sequences. structural rna has higher free energy than expected in figure of, we showed that the free energy e of the minimum free energy structure s of e. coli val-trna, is much higher than the average free energy 〈e〉 of over four million rnas having the same mfe structure s as that of e. coli val-trna. here, e. coli val-trna rv was selected, because its mfe structure s is identical to the rfam consensus structure for trna family rf this preliminary result suggests that naturally occurring transfer rnas may be under selective pressure to be only marginally thermodynamically stable. since it took a number of days for rnaifold to return over four million solutions of the inverse folding problem for the trna target structure, we now describe how rnadualpf can be used to compute the boltzmann expected free energy of literally all sequences a …,a n with respect to an arbitrary target structure s in this manner, we confirm our preliminary finding concerning e. coli val-trna, and show that the folding energy of structural rna from the rfam database is much higher than expected. before presenting results, we need some definitions. for the turner nearest neighbor energy model, the free energy of a secondary structure s of an rna sequence a=a …,a n depends on the temperature t to indicate this dependence, we write e, where in the sequel, t will be designated as table temperature, i.e. the temperature for which parameters from the turner energy tables are applied. for an arbitrary, but fixed secondary structure s of length n, the dual partition function at temperature t is defined by z=∑aexp−ea,s t0/rt where the sum is taken over all rna sequences a=a …,a n of length n. note that t indicates the temperature at which the energy of a structure s and nucleotide sequence a is evaluated using the turner parameters, while all other occurrences of the temperature variable are designated by t, which we call formal temperature. the distinction between formal and table temperature is made to allow us to use finite difference approximations to derivatives with respect to the formal temperature when when we compute dual expected energy and dual conformational entropy below. when table temperature t equals formal temperature t, and the temperature is clear from the context, we write z ∗; if the target structure s is also clear from the context, then we write z ∗. a similar remark applies to the other thermodynamic functions p ∗,g ∗,〈e ∗〉,s ∗, which we now define. the dual boltzmann probability p ∗ is defined by p∗a,s t t=exp−ea,s t0z∗s t t the dual ensemble free energy g ∗ is defined by g∗=g∗=gs t t=−rtlnz∗s t t where r≈ cal/ is the universal gas constant. the dual expected energy 〈e ∗〉 is defined by 〈e∗s t t〉=∑aea,s t0·pa,s t t straightforward derivations analogous to those in yield the following expressions for dual expected energy 〈e ∗〉 and dual entropy s ∗: 〈e∗s t t〉=rt2·∂∂tlnz∗t=t s∗s t t=〈ea,s t t〉−g∗s t tt programs to compute dual expected energy 〈e ∗〉, dual conformational entropy s ∗, and dual heat capacity cp∗ are provided at our web site. we do not elaborate further on dual entropy or dual heat capacity, since at the present time we have found no compelling applications. c, and computed 〈e〉. additionally, for each rfam family, we selected that sequence a whose minimum free energy structure s has smallest base pair distance to the consensus structure s c. we computed the expected energy 〈e〉, as well as the free energies e and e. figure displays box-and-whiskers plots for the fold change 〈e〉e for the consensus structure and the fold change 〈e〉e for the minimum free energy structure. since the dual boltzmann probability p ∗ is generally larger for sequences a having higher gc-content, rnadualpf computes as well the dual partition function for gc-content k, defined by z∗=∑asuch thatgc-content=kexp/rt) fig. analysis of expected free energy 〈e〉 for structures in rfam. given a secondary structure s, the expected free energy of all sequences a with respect to s is defined by 〈e〉=∑ae·exp/rt)z∗, where z ∗ is the dual partition function defined in equation. for each rfam family, we took the family consensus structure s c, and computed 〈e〉. additionally, for each rfam family, we selected that sequence a whose minimum free energy structure s has smallest base pair distance to the consensus structure s c. the expected energy 〈e〉 was computed, as well as the free energies e and e. the fold change 〈e〉e for the consensus structure and the fold change 〈e〉e for the minimum free energy structure were computed. the box-and-whiskers plots show the mean, 25th and 75th percentile, minimum and maximum values. as indicated in the legend, these computations were performed either with respect to all sequences or with respect to all sequences having the same gc-content. these data clearly indicate that natural rna sequences, whose mfe structures most closely resemble the rfam consensus structures, have higher free energy than expected in this fashion, we can exactly compute the dual expected energy 〈e ∗〉 of all sequences having gc-content k which approximately fold into target structure s tables and analyze what we mean by approximately folding into the target structure – i.e. sequences a are preferentially sampled when free energy e is low, hence have large dual boltzmann probability. rnadualpf, even when exact gc-content is controlled, is faster than inverse folding programs by orders of magnitude, hence providing an effective alternative manner of solving inverse folding.  in this paper we describe the algorithm and software rnadualpf, which computes the dual partition function z ∗, defined as the sum of boltzmann factors exp/rt) of all sequences a with respect to the target structure s using rnadualpf, we efficiently sample rna sequences that fold into s where additionally the user can specify iupac sequence constraints at certain positions, and whether to include dangles. moreover, the user can require that all sampled sequences have a precisely specified gc-content, since, optionally, we compute the dual partition function z ∗ simultaneously for all values k=g+c. this sampling strategy is complementary to the use of rnaifold, since it allows the study of the properties of long rna structures whose number of solutions for the inverse folding problem is astronomically large. we use rnadualpf to corroborate previous studies using rnainverse, by confirming that precursor micrornas are significantly mutationally robust when gc-content is not controlled. however, in contrast to, we find that precursor micrornas are significantly non-robust when gc-content is controlled. we confirm and extend previous findings that bacterial small noncoding rnas display plasticity and are not statistically robust, when gc-content is not controlled. additionally, we obtain the new finding that when when gc-content is controlled, bacterial small noncoding rnas are significantly non-robust, as in the case of precursor micrornas. one possible reason for the discrepancy between our results and those of could be related with the fact that the energy parameters of vienna rna package differ from those of vienna rna package.. another possible reason is that the inverse folding solutions returned by the program rnainverse used in show a different bias than sequences returned by rnadualpf. as mentioned in the introduction, there is a relation between our c program rnadualpf and the python program incarnation, although our work is independent of that of reinharz et al.. incarnation is a weighted sampling algorithm that computes the dual partition function for a simple energy model, which only considers base stacking free energies – unlike rnadualpf, the program incarnation includes no energy contributions for hairpins, bulges, internal loops, multiloops, dangles, or mismatches. if the user specifies a desired gc-content α, then incarnation does not compute the dual partition function for gc-content, but rather applies an adjustable heuristic so that after a suitable burn-in period, sequences tend to approximately have gc-content α. see table of for benchmarking results on rnadualpf and incarnation, which show conclusively that rnadualpf is not only faster, but its sequences have a higher probability of folding into the target structure, its sequences have a smaller gc-content in default mode, where gc-content is not controlled, etc. our original motivation in designing rnadualpf was to generate an unbiased sample of near-solutions to the inverse folding problem. at present, it seems clear that no program can claim to generate an unbiased sample of inverse folding solutions, since the solution space so large that this hypothesis cannot be tested by brute force methods, and different inverse folding algorithms return solution sequences having different properties, as shown in table of. nevertheless, in the same manner that structures sampled by the algorithm of ding and lawrence constitute an unbiased, representative set of low energy secondary structures for a given rna sequence, as implemented in sfold and rnasubopt -p, the collection of rna sequences sampled by the algorithm rnadualpf constitute an unbiased, representative set of sequences having low energy with respect to a given target structure s although the minimum free energy structure of such sequences may indeed be distinct from s it is likely that the mfe structure and the target s be similar, as shown in tables and moreover, fig. presents relative frequency plots that suggest that when gc-content is controlled, the sequences returned by rnadualpf have similar properties to those of wild type sequences: similar expected base pair distance to the wild type target structure, similar ensemble defect to the target wild type structure, similar positional entropy, similar vienna structural diversity, similar morgan-higgs diversity, similar expected proportion of native contacts. these graphs were produced by using rnadualpf to sample sequences for each of the c. elegans precursor micrornas from the mirbase database, in each of the following cases: gc-content identical to that of the rfam sequence, no control for the gc-content. figure presents additional data, computed in the same manner for c. elegans pre-mirna from mirbase showing that when gc-content is controlled, sequences sampled by rnadualpf satisfy the following: the average length-normalized base pair distance between the minimum free energy and target structures is ≈, wild type and rnadualpf sampled sequences have similar free energy with respect to the wild type target structure, as well as similar minimum free energy, similar dual probability, and similar probability to wild type rna sequences. taken together, this data shows that if gc-content is controlled, then rnadualpf returns sequences whose low energy structures tend to resemble the target structure. figures and are similar to figs. and except that for each c. elegans pre-mirna, sequences were generated by rnadualpf, each of whose mfe structure is identical to the wild type target structure. taken together, figs. and present convincing evidence that rnadualpf generates sequences that fold into the user-specified target structure, hence supporting our finding that c. elegans precursor micrornas are statistically non-robust, contrary to the finding of. fig. for each of the c. elegans precursor micrornas from mirbase and for each of the following cases, indicated in black, and, indicated in red, rnadualpf sampled sequences without any subsequent filtering step. case - black lines: all rnadualpf sequences have gc-content exactly equal to that of the rfam sequence. case - red lines: rnadualpf was used in default mode, without controlling gc-content. case - blue lines: wild type c. elegans data. density plots are shown for the expected base pair distance to target structure s, the ensemble defect to target structure s, the positional entropy, vienna structural diversity, morgan-higgs diversity, expected proportion of native contacts. all measures were normalized by sequence length fig. additional measures for the data described in the previous fig. density plots are shown for the base pair distance between the minimum free energy structure and the target structure, the gc-content, the free energy e of the rna sequences a with respect to the target structure s the free energy e of each sequence with respect to its own minimum free energy structure, the log dual probability p∗=∑xexp/rt)z∗, and the log probability p=∑sexp/rt)z fig. for each of the c. elegans precursor micrornas from mirbase and for each of the following cases, indicated in black, and, indicated in red, rnadualpf generated sequences, whose minimum free energy structure was identical to that of the corresponding wild type pre-mirna. case - black lines: all rnadualpf sequences have gc-content exactly equal to that of the rfam sequence. case - red lines: rnadualpf was used in default mode, without controlling gc-content. case - blue lines: wild type c. elegans data. density plots are shown for the expected base pair distance to target structure s, the ensemble defect to target structure s, the positional entropy, vienna structural diversity, morgan-higgs diversity, expected proportion of native contacts. all measures were normalized by sequence length fig. additional measures for the data described in the previous fig. density plots are shown for the base pair distance between the minimum free energy structure and the target structure, the gc-content, the free energy e of the rna sequences a with respect to the target structure s the free energy e of each sequence with respect to its own minimum free energy structure, the log dual probability p∗=∑xexp/rt)z∗, and the log probability p=∑sexp/rt)z additionally, we have shown that natural rnas from the rfam database have higher minimum free energy than expected, thus supporting our results in which suggest that functional rnas are under evolutionary pressure to be only marginally thermodynamically stable. the applications described in this paper demonstrate that rnadualpf is a useful and extremely fast software tool for evolutionary and synthetic biology. endnote when dangling positions are not included in the computation, the algorithm clearly requires linear time. when dangling positions are included, run time is exponential in the number of components of the largest multilooop; however, in practice the algorithm is extremely fast, and it is possible to modify the algorithm to always run in linear time. additional file additional file supplementary information. title is supplementary information for rnadualpf: software to compute the dual partition function with sample applications in molecular evolution theory. this is a 3-page pdf file containing a tricky derivation for an efficient computation of the number of external loops of size n with gc-content k, where the user can stipulate that certain positions are constrained to contain nucleotides consistent with iupac codes. abbreviations mfeminimum free energy plmvdpeach latent mosaic viroid sncrnasmall noncoding rna 
rna inverse folding is the problem of finding one or more sequences that fold into a user-specified target structure s ∗, defined as the sum of boltzmann factors exp/rt) of all rna nucleotide sequences a compatible with target structure s we show that c. elegans precursor microrna is significantly non-robust with respect to mutations, by comparing the robustness of each wild type pre-mirna sequence with sequences of the same gc-content generated by rnadualpf, which approximately fold into the wild type target structure. ∗ and samples sequences having low energy with respect to a target structure, allowing sequence constraints and specified gc-content. we confirm and strengthen earlier findings that precursor micrornas and bacterial small noncoding rnas display plasticity, a measure of structural diversity. where additionally the user can specify iupac sequence constraints at certain positions, and whether to include dangles. using rnadualpf, we efficiently sample rna sequences that approximately fold into s ∗〉, and use it to show that natural rnas from the rfam database have higher minimum free energy than expected, thus suggesting that functional rnas are under evolutionary pressure to be only marginally thermodynamically stable. here we consider the ensemble of all rna sequences that have low free energy with respect to a given target s whose minimum free energy secondary structure is identical to the target s rna secondary structurepartition functionboltzmann ensemblerobustnesshttp://dx.doi.org/ /100000001national science foundationdbi-1262439clote peter issue-copyright-statement© the author 2016 ∗ over all sequences having gc-content k, the user can require that all sampled sequences have a precise, specified gc-content. using different inverse folding software, another group had earlier shown that pre-mirna is mutationally robust, even controlling for compositional bias.
 rna inverse folding is the problem of finding one or more sequences that fold into a user-specified target structure s i.e. whose minimum free energy secondary structure is identical to the target s here we consider the ensemble of all rna sequences that have low free energy with respect to a given target s  we introduce the program rnadualpf, which computes the dual partition function z ∗, defined as the sum of boltzmann factors exp/rt) of all rna nucleotide sequences a compatible with target structure s using rnadualpf, we efficiently sample rna sequences that approximately fold into s where additionally the user can specify iupac sequence constraints at certain positions, and whether to include dangles. moreover, since we also compute the dual partition function z ∗ over all sequences having gc-content k, the user can require that all sampled sequences have a precise, specified gc-content. using z ∗, we compute the dual expected energy 〈e ∗〉, and use it to show that natural rnas from the rfam database have higher minimum free energy than expected, thus suggesting that functional rnas are under evolutionary pressure to be only marginally thermodynamically stable. we show that c. elegans precursor microrna is significantly non-robust with respect to mutations, by comparing the robustness of each wild type pre-mirna sequence with sequences of the same gc-content generated by rnadualpf, which approximately fold into the wild type target structure. we confirm and strengthen earlier findings that precursor micrornas and bacterial small noncoding rnas display plasticity, a measure of structural diversity.  we describe rnadualpf, which rapidly computes the dual partition function z ∗ and samples sequences having low energy with respect to a target structure, allowing sequence constraints and specified gc-content. using different inverse folding software, another group had earlier shown that pre-mirna is mutationally robust, even controlling for compositional bias. our opposite conclusion suggests a cautionary note that computationally based insights into molecular evolution may heavily depend on the software used. c/c++-software for rnadualpf is available at http://bioinformatics.bc.edu/clotelab/rnadualpf. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords rna secondary structurepartition functionboltzmann ensemblerobustnesshttp://dx.doi.org/ /100000001national science foundationdbi-1262439clote peter issue-copyright-statement© the author 2016 