For a given sncRNA sequence (of length *L*), there is a thermodynamic ensemble (Ω) that contains the different suboptimal structures, each with a given free energy (*G*
_{
i
}) [50]. Thus, the partition function reads
, and the free energy of the ensemble is *G* = -*kT* ln(*Z*). Then, the probability that the sncRNA folds into the structure *i* is given by
. We assumed *T* = 37°C, then *kT* = 0.616 Kcal/mol. In this work, instead of comparing the MFE structures to analyze two different sequences, we compared the two ensembles of structures corresponding to the sequences. We introduced the base-pair distance between two structures (*d*
_{
BP
}), which is more accurate than the Hamming distance, to evaluate the difference between two structures [26]. The base-pair distance (*d*
_{
BP
}) between the different structures of Ω (*S*
_{
i
} denotes structure *i*), referred as intrinsic distance, is given by
(doubly probabilistically averaged). This magnitude accounts for the structural variability within Ω of a given sequence, and then allows defining plasticity (*P*). Lower values of *d*
_{
0
} indicate that Ω is dominated by the MFE structure, while higher values correspond to ensembles with more suboptimal structures within a given energy gap. More plastic is a sequence when it presents more structural fluctuations at the equilibrium. Therefore, we defined plasticity as
. This magnitude can then be used to distinguish very stable RNAs.

To compute *R*
_{
m
} we need to compare different mutant sequences. The average distance between structural ensembles after one single point mutation (*d*
_{
1
}) follows
(where Ω_{1} is the ensemble of mutants and Π_{
j
} is calculated using the partition function of Ω_{1}, denoted by *Z*
_{
1
}). Since *d*
_{
1
} only accounts for one mutant, we need to average several calculations. Here 〈•〉 indicates average for perturbations and Δ• standard deviation. Hence, 〈*d*
_{1}〉 is the average structural distance after 1 single point mutation (*L* calculations). Then, we defined mutational robustness as
. As in the definition of *P*, we rescaled by *L*/2 to have an absolute value, since 〈*d*
_{
i
}〉 scales with *L* and because the number of base-pairs of a structure is bounded by *L*/2 (i.e., *d*
_{
BP
} between certain structure and the unfolded state is bounded by *L*/2). Certainly, the lower the distance, closer to 1 (maximum) should be the robustness. Here we considered that robustness follows a linear trend with the relative structural distance, although quadratic expressions could also be employed [13]. Analogously, we calculated the distance between structural ensembles after 2 single point mutations (*d*
_{
2
}), and the distance between ensembles after one environmental perturbation (*d*
_{
e
}), which was simulated as a random variation over the value of all the energetic parameters that define the model for base-pair interactions. For that, all parameters determining the energies for base pairing and stacking are perturbed simultaneously [22]. More in detail, to perform environmental perturbations, we took variations up to 20% of the nominal values following normal random distributions, i.e., being *β*
_{
0
} the nominal value of an energetic parameter, the perturbed value reads *β* = (1+0.2*ξ*)*β*
_{0}, where *ξ* ~ Ν(0,1). Hence, 〈*d*
_{2}〉 is the average structural distance after 2 single point mutations (10 *L* calculations), and 〈*d*
_{
e
}〉 is the average structural distance after an environmental perturbation (1,000 calculations). Then, we defined environmental robustness as
. We further defined epistasis as
, which measures the interference between mutations. *E* > 0 means antagonistic epistasis (i.e., 〈*d*
_{2}〉 < 2〈*d*
_{1}〉, resulting in compensatory effects), while *E* < 0 synergistic epistasis (i.e., 〈*d*
_{2}〉 > 2〈*d*
_{1}〉, resulting in enhancement effects) [15].

In addition, for each sncRNA we computed its degree of functionality (*V*), given by
(probabilistically averaged), where *V*
_{
i
} is the number of times that a motif involving consecutively three free nucleotides and three bound nucleotides (in the 5' sense or in the 3') appears in the structure *i* of the ensemble. Two overlapping motifs were counted as a single event. While *V*
_{
i
} is a magnitude that corresponds to one structure, *V* corresponds to a sequence. We call this magnitude functionality because it quantifies the number of different mechanisms for potential interactions with further RNA molecules [2, 42]. In other words, the degree of functionality accounts for the number of regions that may provide accessibility for RNA-RNA interactions. Moreover, *V*
_{
i
} is roughly proportional to the number of hairpins of the structure, and that metric of functionality also accounts for the complexity of the molecule.

Structural robustness was tested for significance by comparing it to a distribution of robustness values generated from a large set of artificially originated sequences. Artificial sequences shared the property of yielding the same MFE structures as the real sequences. For each sncRNA, we generated 69 random sequences, resulting in a population of 5,451 elements. Results were primary maintained when using smaller random populations. We constructed three different random samples. Sample I was obtained by iteratively solving the corresponding inverse folding problems using different initial sequences [10] with Vienna RNA package (default energetic parameters, dangles = 2, MFE objective) [23, 51]. Notably, this allows sharing the MFE structure, but the thermodynamic ensembles may differ. Subsequently, sample II was obtained by combining inverse folding and neutral evolution, introducing mutations that do not change the MFE structure [12], thereby minimizing the bias introduced by the optimization method itself. For that, we performed *L* mutations. This process, nevertheless still produces biased sequences because mutations would be accumulated in regions with unpaired nucleotides (e.g., loops or tails). By definition, mutations affecting paired nucleotides are not neutral, with the exception of G-U/G-C paired regions. To counterbalance this bias, we constructed a sample III by which sequences were subjected to a neutral evolution process accounting for potential compensatory mutations (also *L* mutations). This process was based, in the case of paired nucleotides, on mutating the complementary nucleotides as well. Following this procedure, the simulated neutral evolution process accounts for both neutral single-point mutations and neutral base-pair mutations. This allowed enlarging considerably the sequence space and avoid more efficiently the bias produced by inverse folding methods.