The energy-spectrum of bicompatible sequences

Genotype-phenotype maps provide a meaningful filtration of sequence space and RNA secondary structures are particular such phenotypes. Compatible sequences, which satisfy the base-pairing constraints of a given RNA structure, play an important role in the context of neutral evolution. Sequences that are simultaneously compatible with two given structures (bicompatible sequences), are beacons in phenotypic transitions, induced by erroneously replicating populations of RNA sequences. RNA riboswitches, which are capable of expressing two distinct secondary structures without changing the underlying sequence, are one example of bicompatible sequences in living organisms. We present a full loop energy model Boltzmann sampler of bicompatible sequences for pairs of structures. The sequence sampler employs a dynamic programming routine whose time complexity is polynomial when assuming the maximum number of exposed vertices, κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document}, is a constant. The parameter κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} depends on the two structures and can be very large. We introduce a novel topological framework encapsulating the relations between loops that sheds light on the understanding of κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document}. Based on this framework, we give an algorithm to sample sequences with minimum κ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} on a particular topologically classified case as well as giving hints to the solution in the other cases. As a result, we utilize our sequence sampler to study some established riboswitches. Our analysis of riboswitch sequences shows that a pair of structures needs to satisfy key properties in order to facilitate phenotypic transitions and that pairs of random structures are unlikely to do so. Our analysis observes a distinct signature of riboswitch sequences, suggesting a new criterion for identifying native sequences and sequences subjected to evolutionary pressure. Our free software is available at: https://github.com/FenixHuang667/Bifold.


Background
Bicompatible sequences in evolution RNA evolution has been studied extensively in the framework of theoretical evolutionary optimization, center staging the genotype-phenotype mapping from RNA sequences to their structures [1][2][3][4][5][6][7]. RNA secondary structures are particular such phenotypes, contact structures that can be represented as diagrams, where vertices are nucleotides and arcs are base pairs drawn in the upper half-plane. If the arcs are not crossing, the represented RNA secondary structure is pseudoknot-free [8,9]. These pseudoknot-free structures correspond to tree structures based on the lengths of contiguous subsequences. As a result, finding a structure with minimum free energy (mfe) to a given sequence can be computed efficiently using dynamic programming (DP) routines [10,11]. Diagrams containing crossing arcs represent RNA pseudoknots [12,13]. These pseudoknot structures can not be decomposed into trees in a straightforward Page 2 of 18 Huang et al. Algorithms Mol Biol (2021) 16:7 fashion. Recently, pseudoknots have been studied from a topological perspective and their structural complexity was characterized by the topological genus [14][15][16][17]. In our paper, we are interested in RNA riboswitches, a small segment of an mRNA molecule that regulates mRNA translation and can express two alternative secondary structures. As riboswitches are small ( < 500 nts), they are not likely to exhibit pseudoknots [13,18,19]. In the following we restrict our analysis to pseudoknot-free secondary structures.
In [4], the authors realized that genotype-phenotype mappings provide a natural filtration of the sequence space by means of considering sequences "equivalent" if they fold into the same mfe-secondary structure. This perspective naturally leads to consider the induced subgraphs of preimages of the folding map in sequence space, i.e., the neutral networks of RNA secondary structures [5]. Neutral network are graphs, consisting of all sequences that fold into a distinguished structure, with edges connecting two such sequences if they differ in exactly one nucleotide. These networks provide a framework to quantify well-known evolutionary theories, such as Motoo Kimura's neutral theory of evolution.
A plethora of work has been done on the diffusionlike process of sequences searching for an optimal structure, ranging from simulation-based studies [20] to the mathematical analysis of the cluster-size distribution depending on the structure of the neutral net [5]. These studies have shown that connectivity and density of neutral networks are of central importance for the understanding of how sequences evolve.
One prominent phenomenon is that of spontaneous, rapid transitions of evolving populations of RNA sequences from one structure to another-even in absence of fitness advantages [6,21]. Despite of the Intersection Theorem [5], guaranteeing the existence of bicompatible sequences for any two RNA secondary structures, transitions between neutral networks are only observed for certain structure pairs. In [21], the authors showed that in the course of a phenotypic transition an evolving population of RNA sequences tunnels through bicompatible sequences. These sequences constitute de facto a gateway between different phenotypes.
Bicompatible sequences play furthermore a prominent role in the analysis of RNA riboswitches [22]. Riboswitch sequences express two distinct structures, each of which appearing in a specific biophysical contexts, see Fig. 1. Both structures are typically thermodynamically suboptimal and exhibit a large base pair distance [23]. Specific mechanisms are observed, most prominently that of the existence of a switching sequence, a contiguous subsequence that engages for each respective structure in a unique fashion. The two structures are mutually exclusive, since bases of the switching sequence pair downstream in one and upstream in the other configuration [22].  [24] and its switching sequence (blue), involved two respective helices

Finding bicompatible sequences for two structures
The studies of bicompatible sequences of a pair of secondary structures of the same length, or in more general cases where a sequence satisfies multiple structure constraints, are motivated by the computational design of RNA sequences [25][26][27][28][29]. Early developments such as [11,25] design a sequence by simulation approaches. These entail considering an objective function that involves the energy contribution of multiple target structures on a common sequence. They then replicate a sequence by introducing a single random mutation, and a sequence survives if it improves the performance of the designated objective function. Later developments such as [26,27] generate sequences satisfying multiple structure constraints based on a multi-objective genetic algorithm.
Recently, new approaches such as [30][31][32][33], design sequences on a single secondary structure using a Boltzmann sequence sampler based on the "dual" partition function of with respect to McCaskill's partition function of secondary structures for a fixed sequence [34]. These approaches consider a sequence and a structure as a pair, and the pair can be mapped via Here, η is the energy function of a sequence-structure pair (σ , S) , Q n 4 and S n are the collections of all sequences and secondary structures of length n respectively, K is a constant(the Boltzmann constant), and T is a temperature parameter. The Eq. 1 resembles the notion of a scalar product of vector spaces: V × V → K , whence the notion of dual partition function. The partition function for a fixed structure induces a probability space with Boltzmann distribution where a sequence has a higher probability if it is energetically more stable with respect to the given structure.
For multiple structures, [35,36] made use of a combination of graph coloring theory and heuristic local optimization in order to find sequences whose energy landscapes are dominated by the prescribed conformations. Later development in [28] presents an algorithm to sample sequences with multiple structure constraints with uniform probability.
In [29], the authors developed a Boltzmann sampler to generate sequences for two secondary structures of the same length such that the total energy of the two sequence-structure pairs is minimized. Their sampler is based on a simplified energy model focusing on the energy contribution of helices, while the energy contributions from hairpin loops, interior loops, and multi-loops are ignored. The key challenge is to decompose the loops contained in the two structures in a tree structure such that a DP-routine can be employed to (1) compute their partition function. In their approach, a hyper-graph model is used to describe the intersection relationships among the loops. Assume a tree decomposition of the hyper-graph is given with tree-width w, then the time complexity of computing the partition function is O(4 w+1 n) , thus fixed-parameter tractable (FPT). However, in the general case, finding a tree decomposition of a hyper-graph with minimum tree-width is NPhard. Bodlaender's famed result [37] shows that checking whether a tree decomposition with tree-width ≤ k exists and constructing a tree decomposition with tree-width k, if such k exists, can be performed in linear time. However, it remains unknown how big the gap between such k and the minimum tree-width is. As the sampler in [29] considers only simple loops formed by helices, the induced hyper-graph is usually simple, whence an optimal tree decomposition can be approximated. Though the authors in [29] claim their approach can be generalized to a full-loop energy model [38], the hyper-graph will become more involved, and an optimal tree decomposition is arguably no longer easy to obtain.

Boltzmann sampling
In order to understand how the time complexity is affected when considering a full-loop energy model [38], we express the loop intersection relationship by means of a novel topological framework [39]. The framework considers each loop as a vertex (a 0-simplex), and d loops having nontrivial intersection as a(d − 1)-simplex. Thus, the simplicial complex obtained by gluing all d-simplices, for all d ≥ 0 , gives rise to a topological space. The simplicial complex is called a loop nerve, see Fig. 2. Specifically, the 0-simplices representing the loops are equivalent to the hyper-edges in the hyper-graph model discussed in [29], while the d-simplices for d > 0 capture additional information that is not present in the hyper-graph model. We observe that the time complexity of computing the partition function is closely related to the structure of the topological space, in particular whether or not it contains d-dimensional holes. In topology, d-dimensional holes are elements of the dth homology group and the latter constitute a key signature of the space. The easiest way to understand the significance of this is to consider a disc and a punctured disc. While the former can continuously be contracted into a point, the latter does not allow for such a contraction. We illustrate the geometric interpretation of the 0th-, 1st-, and 2nd-homology group in Fig. 3. Homological signatures play a central role in the context of differentiating topological spaces, since any homeomorphism of topological spaces induces an isomorphism of homologies. For instance, the characterization of closed surfaces as a sphere, a finite connected sums of projective planes, or a finite connected sums of tori uses the fact that the homologies of spheres, projective planes and tori are distinct.
In [39] it is shown that the loop nerve induced by a pair of secondary structures has an associated topological quotient space of a very particular type: it is obtained by gluing tetrahedra, spheres and ribbons. It turns out that the tetrahedra can be processed via homotopies and the resulting space may be depicted as an apple tree, the spheres representing the apples and the ribbons the branches. In the absence of spheres, we shall obtain an optimal algorithm to compute the partition function of a pair of secondary structures, i.e., we can obtain in this case a tree decomposition of the hyper-graph in [29] with minimal tree-width. It is the spheres, that are responsible for the computational complexity and an optimal polynomial time algorithm is at present unknown.
The topological framework introduced does, at present, not solve the NP-hardness rooted in the tree decomposition with minimum tree-width, it does however present a fundamentally different approach to the study of computational complexity problems associated with computing mfe-sequences to two or more secondary structures. One might speculate that simplicial complexes specifically the recently developed framework of weighted complexes (a generalization in which simplices are endowed with a certain weight) may prove useful for larger classes of computational complexity problems.
The framework presented here is tailored to express intrinsic relations between loops of different structures and reveals a detailed blueprint of how they organize, namely as topological wedge-sums of spheres.
It is well-known [40,41] that even identical genotypes can lead to many different phenotypes in response to environmental changes, and riboswitches are natural sequences that can easily access a variety of phenotypes. Phenotypic accessibility is closely connected to bicompatible sequences and we shall employ the sequence sampler for two secondary structures in a full-loop energy model [38] to investigate properties of RNA riboswitch sequences on multiple levels. First, we compare the energy spectra of sample sequences under the constraint of being compatible and bicompatible, respectively. We show that for the two alternative structures of a riboswitch, the energy spectra of compatible sequences with either one structure is similar to the sequences that are bicompatible, i.e. compatible with both structures. This implies that riboswitch sequences can switch between their assumed structures easily, without affecting the free energy since there exist multiple, energy favorable bicompatible sequences. In contrast, random structure pairs do not facilitate such phenotypic switches, i.e. native riboswitches are distinguished sequences.
Second, we analyze how sequence-structure pairs rank within the partition function by comparing riboswitch sequences with random sequences. The rank analysis allows us to conclude that the relative rank of native riboswitch sequences is distinctively higher than the relative rank of the sampled sequences. Accordingly, native sequence exhibit an optimized thermodynamic stability with respect to the pair of structures.
Thirdly, we investigate a certain adaptability index, measuring the capability of a structure R to change into the structure S. We find that the adaptability of the alternative structures of riboswitches are distinctively different from that of random structure pairs. This indicates that the two alternative structures of a riboswitch are more evolutionary accessible by sequences than random structure pairs.

Bistructures
We present an RNA secondary structure as an arc diagram, a graph whose vertices are drawn on a horizontal line and the Watson-Crick as well as Wobble base pairs are drawn as arcs in the upper half-plane [42][43][44], see Fig. 4B. The vertices are labeled by V = {1, 2, . . . , n} from left to right, representing the nucleotides. The linear order of the vertices indicates the direction of the RNA strand from 5 ′ -end to 3 ′ -end. Here we consider only the canonical Watson-Crick and Wobble base pairs in an RNA secondary structure. As a result, for any pair of nucleotides, there can be at most one such canonical base pair, each vertex can be only incident to one arc.
An arc, (i, j), represents the base pair between the ith and jth nucleotides. Two arcs (i, j) and (r, s) are called crossing if and only if i < r < j < s holds. An RNA structure is called a secondary structure, if it does not contain any crossing arcs. Furthermore, the arcs of a secondary structure can be endowed with the partial order: (r, s) ≺ (i, j) if and only if i < r < s < j . We shall introduce two "formal" vertices associated with positions 0 and (n + 1) , respectively and add the formal arc (0, n + 1) , referred to as the rainbow. An interval, [i, j], is the set of vertices {i, i + 1, . . . , j − 1, j}.  (1,17). E A bistructure B is a collection of loops {L 1 , . . . , L 9 } . X = {L 3 , L 7 , L 9 } (blue) is an irreducible substructure of B with its complement X = {L 1 , L 2 , L 4 , L 5 , L 6 , L 8 } . We mark the exposed vertices E X = V X ∩ V X in red. The closure of X is given by X = {L 2 , L 3 , L 4 , L 6 , L 7 , L 8 , L 9 } In a loop-based energy model [38,45], arcs and unpaired vertices are organized in loops contributing to the energy. A loop, L, is a subset of vertices, represented as a disjoint union of S-intervals, are arcs (including the rainbow arc (0, n + 1) ) and where any other interval-vertices are unpaired, see Fig. 4D. It can be represented by a maximal arc (a 1 , b k ) with respect to the partial order ≺ . Given a loop, this maximal arc is unique, whence a loop can be represented by L (a 1 ,b k ) . In particular, the rainbow arc, (0, n + 1) , represents an exterior loop, that is not nested in any arc in the arc diagram, see Fig. 4C. Furthermore, each non-rainbow arc appears in exactly two loops, being maximal for exactly one of them. Loops correspond to the boundary components of the secondary structure viewed as a fatgraph [46]. In the following, we shall identify loops with their sets of vertices.
Given two secondary structures, R and S, having the same vertex set V = {1, . . . , n} , we draw the vertices on a horizontal line, the arcs of R in the upper and the arcs of S in the lower half-plane. We refer to this arc diagram as a bistructure, B(R, S). The idea of considering a secondary structure pair of given length has been studied in [47,48]. Here we shall distinguish the R-arcs from the S-arcs even though they might have the exact same endpoints. For example an arc (i, j) in R is denoted by (i, j) R and an arc (i, j) in S is denoted by (i, j) S . In a loop-based model, the R-loops and the S-loops are distinct since their represented R-arcs from the S-arcs are distinct. Hence, a bistructure B(R, S) can be considered as the set of loops is not necessarily empty, since paired and unpaired vertices can be contained in the intersection of the B ′ -and B ′ -loops. Furthermore, for a given substructure X = {L 1 , · · · , L k } , we define the boundary of X by X C = {L ∈ X|∃L i ∈ X, L ∩ L i � = ∅} , i.e., X C is the set of all loops in the complement of X that have nontrivial intersection with X. We call X = X ∪ X C the closure of X. A substructure is called reducible if the loop set can be bi-partitioned into two sets of loops The intersection E B ′ = V B ′ ∩ V B ′ is called the set of exposed vertices of B ′ . The exposed vertices are key elements in computing the partition function of a bistructure, since the vertices are contained in multiple loops and their nucleotide information needs to be remembered until the energies of the loops containing the exposed vertices are calculated.

Partition function and Boltzmann sampler
We first recall the notion of a partition function for sequences that are compatible to a single structure R [34].
Here C n (R) denotes the set of R-compatible sequences while η(σ , R) is the energy of the sequence-structure pair (σ , R) . Lastly, K is the Boltzmann constant and T the temperature. In Turner's model [38,45], η(σ , R) = L∈R η(σ , L) , where L is a loop contained in the secondary structure R. The energy of a loop L is a function of its type and of the nucleotides associated to the arcs and the unpaired bases it contains. In practice, the energy computation takes into account a maximum of two specific arcs and four unpaired vertices, as well as the number of arcs and the number of unpaired bases.
For a bistructure B(R, S) and a sequence σ , we set η(σ , B(R, S)) = 1 2 (η(σ , R) + η(σ , S)) . Then we define the partition function of sequences bicompatible to R and S by where C n (R, S) denotes the set of bicompatible sequences to both R and S.
A decomposition of B is a block sequential loop removal of the bistructure. Let us first illustrate the computation of Q(R, S) when a specific decomposition is given. Suppose X = {L 1 , . . . , L k } is a substructure of B(R, S) with vertex set V, and exposed vertex set E X .
Then we can compute the energy η(σ X , X) since the nucleotide information of the vertices contained in V is specified.
Here, � ℓ is the collection of RNA sequences of length ℓ.
By definition, if X is an irreducible substructure, then removing a loop L from X produces a set of irreducible substructures X 1 , . . . , X k . We investigate how the exposed vertex set evolves with a loop removal. To this end let x ∈ E X be an exposed vertex. If x ∈ L , then either KT . (a) ∃L ′ ∈ X, L ′ = L such that x ∈ L ′ , or, (b) at least one such L ′ loop exists. In the first case (a), we have x is no longer exposed, while in the second case (b), we have x ∈ E X i for some 1 ≤ i ≤ k . Finally, if x / ∈ L to begin with, then we have x ∈ E X i for some 1 ≤ i ≤ k after removing L form X.
Let τ X denote a fixed subsequence over E X , τ X i a subsequence over E X i , 1 ≤ i ≤ k , and σ L a subsequence over the loop L. We consider all possible subsequences Then, the partition function Q(X, τ X ) can be computed recursively by For a given decomposition, the terms Q(X i , τ X i ) , for 1 ≤ i ≤ k , can be computed in parallel. We illustrate the recursion 2 in Fig. 5.
When Q(R, S) is computed, we can Boltzmann sample RNA sequences following the classical stochastic backtracking method introduced by [49], which is of linear time complexity. Given an irreducible substructure X that is decomposed into a loop L and a set of irreducible substructures X 1 , . . . , X k . Assume the nucleotides in X i are sampled, then with a fixed subsequence τ X over the exposed vertex set E X , the subsequence (σ v ) v , v ∈ L \ ∪ i E X i is sampled with probability Multiplying all inside probabilities of each iteration, we conclude that a sequence is sampled with probability P(σ ) = e − η(σ ,B(R,S)) KT /Q(R, S).

The topology of a bistructure
The partition function Q(R, S) is computed recursively based on substructures, where a loop is removed from the substructure to calculate its energy. The removal yields a collection of substructures having fewer loops.
In this recursion, a nucleotide τ i at position i needs to be stored until the energy of all loops containing τ i are calculated. For a fixed decomposition D = {X 0 , . . . , X k } , the number of nucleotides that need to be stored at each step of the recursion 2 is |L ∪ k i=1 E X i | . We denote the maximum number of |L ∪ k i=1 E X i | in all recursion steps by κ D (B) . For a fixed decomposition D, κ D (B) is a constant.
Assume D is a decomposition of B(R, S), we can implement a dynamic programming (DP) routine to compute Q(R, S) recursively. The time complexity of the algorithm It is impossible to consider all decomposition since there are exponentially many of them. The algorithm introduced in [29] computed the partition function following an analogous DP-routine as Eq. 2. The key question is then to find a "smart" decomposition of the bistructure. The authors in [29] develop a hyper-graph model to interpret the overlapping relationships among all loops. In their hyper-graph model, a labeled vertex represents a nucleotide at a fixed position, and a hyperedge represents a loop. Then a tree decomposition of the hyper-graph induces a hierarchy tree structure for the loops. Following the tree decomposition, one can derive a removal order of loops in the bistructure, and by construction the tree-width w = κ D (B) − 1.
Finding a tree decomposition with minimum treewidth for a general hyper-graph is NP-hard. However, for the simple energy model in [29], the hyper-graph may be simple enough such that a tree decomposition with minimum tree-width can be found by approximation algorithms [37]. It is not clear whether this is still feasible when passing to a full-loop energy model. Fig. 5 Illustration of the recursion for the partition function. Given a bistructure X with exposed vertices v 1 and v 2 , we decomposed X into a loop L(v 1 , u 1 , u 2 , u 3 ) and a substructure X ′ (v 1 , u 1 , u 2 , u 3 , v 2 ) . Here we assume v 1 , u 1 , u 2 , u 3 , and v 2 are contained in loops of X ′ , hence v 1 and v 2 remain exposed and u 1 , u 2 , and u 3 are newly exposed vertices Loop intersections are studied in [39] via a simplicial complex, where a loop in a bistructure B(R, S) is represented by an abstract 0-simplex, i.e., a vertex. This line of work goes beyond the hyper-graph approach in that intersections of multiple loops can be consistently expressed and unlocks powerful concepts from algebraic topology. If d loops have nonempty mutual intersection, they are represented by a (d − 1)-simplex. The collection of all d ≥ 0-simplices forms a simplicial complex, called a loop nerve, see Fig. 2. A loop removal is tantamount to deleting the corresponding 0-simplex as well as all higher dimensional simplices that contain it. We shall show that understanding the structure of the topological space provides insight into designing an optimal decomposition.
We first give an overview of how to design the decomposition of a bistructure via the topological framework. In [39], the authors show that the induced loop nerve of a bistructure B(R, S) has very specific properties. The topological space is uniquely classified by the rank of its second homology group r 2 (B) , which counts the number of 2-dimensional holes. As mentioned before, the geometric realization is comprised of ribbons glued to filled tetrahedra and spheres. Each sphere has a combinatorial interpretation within the bistructure as a crossing component and their number equals the rank of the second homology group. We shall show that the challenge of the decomposition problem stems from the spheres, as the ribbons and tetrahedra are organized in a tree-like fashion. The global tree-like structure induces a tree decomposition naturally, while the sphere can be resolved locally. To resolve the spheres we can map the problem to a known NP-problem such as, for instance, the traveling salesman problem (TSP). This allows us to solve the spheres via approximation algorithms [50]. We illustrate this idea in Fig. 6.

Topological framework
In the following we discuss the results in [39].  correspond to hyper-edges in [29]. In the loop nerve the collection of d-simplices, encapsulates the information of loop intersections not articulated in the hyper-graph model. The loop nerve of a bistructure B contains no d-simplex for d > 3 [39] and there are only two nontrivial homology groups of T(B), both being free and abelian: H 0 (T (B)) ∼ = Z and H 2 (T (B)) ∼ = ⊕ r k Z . In view of connectivity, the rank of H 2 (T (B)) , r 2 (B) , is the only determinant. An R-arc (i, j) and an S-arc (r, s) are called crossing if i < r < j < s holds. We shall proceed by discussing overlaps and crossing components.
An overlap is a degree four vertex in its arc diagram. An overlap corresponds to a 3-simplex in K 3 (X) in the loop nerve. Assume x is an overlap being the endpoint of the arcs p 1 ∈ R and p 2 ∈ S . We split x into two adjacent vertices x 1 and x 2 , where x 1 carries the endpoint of p 1 and x 2 the endpoint of p 2 . This is done such that after the split p 1 does not cross p 2 , see Fig. 8 We next convince ourselves that the split does not "really" affect the induced topological space, where "really" means "upto homotopy". Let x be an overlap. This vertex is contained in four loops and is the endpoint of two arcs p 1 and p 2 . Let L 1 and L 2 be the loops in R that contain p 1 , where p 1 is the maximal arc of L 2 . Furthermore let L 3 and L 4 be the loops in S that contain p 2 , where p 2 is the maximal arc of L 4 . Clearly, ∩ 4 i=1 L i = {x} . Splitting x into x 1 and x 2 [39] is tantamount to removing an edge of the corresponding filled tetrahedron as well as its interior. We end up with two triangles that are still glued along the opposite edge from the edge we removed, see Fig. 8. This splitting does not change r 2 (B) , whence we can restrict ourselves to the non-overlap case.
In the arc diagram of a bistructure we adopt the notion of crossing component as in [39]. We illustrate the * -graph of a bistructure without overlaps and crossing arcs in Fig. 9. Note that, by Lemma 3, if B has no crossings, the induced topological space T(B) is a "ribbon tree". Namely, each ribbon is obtained by gluing a sequence of triangles along their edges such that each triangle has at most two edges glued to other triangles. These ribbons are then glued together along some of the edges of their constituent triangles such that no closed bands appear. Now we are in position to describe the structure of the topological space T(B). If an irreducible substructure X is induced by a crossing component, then the induced topological space is "sphere"-like. Otherwise if X is noncrossing, the induced topological space is "ribbon tree"-like. T(B) is a ribbon tree modulo edge contraction of spheres [39]. Finally, we have the combinatorial interpretation of r 2 (B) and given a bistructure B(R, S) with r crossing components. then r 2 (B) = r.

Scheduling
We next discuss how to design a decomposition based on the properties of the loop nerve. The global treelike structure induces a tree decomposition naturally, while the sphere will be resolved locally. We first consider the case where B contains no crossing arc. In this case, we extend the partial order ≺ for a bistructure by the following: for any two arcs (i, j), (r, s) ∈ B we say (i, j) ≺ B (r, s) if and only if i < r < s < j . Then, we show in the SM that for an irreducible substructure X ⊆ B , X contains a unique maximal arc with respect to ≺ B .
We decompose X by removing the loop L m , where m is the maximal arc of X. The loop removal produces a set of irreducible substructure X 1 , . . . , X k . Repeating this loop removal for any produced irreducible substructures gives a unique loop removal order D 0 (B) . We show A bistructure B with overlaps can be mapped to a bistructure B ′ without overlaps by the above splitting of overlapping vertices. The decomposition D 0 on B ′ induces a natural decomposition D on B by the one-toone correspondence between the B-arcs and the B ′ -arcs. We illustrate in Fig. 10 how to derive a decomposition of a hyper-graph with minimum tree-width for the case where B(R, S) is a bistructure without crossing arcs or overlaps.
We next discuss how to resolve spheres. Recall that the NP-hardness of the decomposition problem stems from the spheres. In this case, we consider mapping the problem to a known NP-problem as, for instance, the traveling salesman problem (TSP). To this end we remove a set of loops from X with a minimum number of exposed vertices, such that X has no crossing arcs. The remaining noncrossing substructure can be decomposed using the optimal algorithm presented before, see Fig. 6. This allows to solve the problem via approximation algorithms of the TSP [50]. The approximation approach is the subject of future work and beyond the scope of this paper, for the analysis presented below, we employ a greedy approach to resolve the spheres.

Structural adaptability
In this section, we introduce three measures for quantifying the structural adaptability using the bicompatible sequence sampler. Given two structures R and S, we first sample sequences with single-compatible and bicompatible constraints respectively, and compare their energy spectrum of sampled sequences paired with R and S. In case of bicompatibility not affecting the energy spectrum significantly, we can conclude that switching from R to S (and vice versa) is feasible. Secondly, we investigate the energy ranking of (σ , R) and (σ , S) , within the partition function of σ . This shows how stable R and S are in the Boltzmann ensemble of a sampled sequence. Finally, we introduce an index, called the adaptability, measuring the capability of a structure R to transform into S. The adaptability is obtained by comparing the proportion of the partition function for bicompatible sequences w.r.t. the partition function of single-compatible sequences. In case of all sequences being bicompatible, the index equals 1, while in the case of sequences not being bicompatible, the index equals 0.

Energy spectra
Let us begin by introducing the spectrum over a partition function. Let Q(X) be a partition function of sequences compatible with X, where X is a secondary or bistructure and Q(X) = σ e − η(σ ,X) KT . To simplify notation, we shall write Q instead of Q(X), if we do not need to emphasize the context of the underlying structure X. Naturally, Q induces the discrete probability space (� n , P Q ) , where P Q (σ ) = e − η(σ ,X) KT /Q . We consider a real-valued random variable f : (� n , P Q ) −→ R and refer to the induced measure P f on R , P f (r) = {σ |f (σ )=r} P Q (σ ) , as the f-spectrum over Q.
For practical purposes, an f-spectrum, P f (r) , cannot be computed directly, since we have to consider all σ ∈ � n and potentially infinitely many r ∈ R . To approximate the f-spectrum we first discretize by means of a monotone increasing multiset (a s ) , where = a s − a s−1 , setting P f (a s ) = a s−1 <r≤a s P f (r) . We employ a Boltzmann sampler to generate sequences of the probability space (Q n 4 , P Q ) and approximate P f (a s ) by P f (a s ) ≈ 1 m | {σ | a s−1 < f (σ ) ≤ a s }| , where σ is a sequence sampled from the partition function Q and m denotes the sample size. Here we set m = 10 4 .
We proceed by introducing some particular choices for the pair (f, Q), which we shall denote by f Q : We call P f R Q (r) the R-spectrum of Q and P f S Q (r) the S-spectrum of Q.

Ranking
Next we investigate how stable R and S are in the Boltzmann ensemble of Q(σ ) , where σ ∈ (� n , P Q(R,S) ) is a sequence sampled via Q(R, S). We compare the energies, η(σ , R) and η(σ , S) to η(σ , M(σ )) , where M(σ ) denotes the mfe-structure of σ . Then we consider the ratios For a fixed sequence σ , the ratios reflect the gap between the energies obtained when considered with R (or with S) and the mfe-structure.

Adaptability
We discuss the energy-spectrum over a partition function, Q as an induced measure of a random variable, f. By construction we normalize, when working with the probability measure P Q (σ ) , the value of Q. As a result, the absolute values of the different partition functions, for instance, when comparing Q(R) and Q(R)| S is not a factor.
Comparing a plethora of riboswitches, as well as sequences of various random structure pairs, we end up with relating the partition functions of sequences over an entire spectrum of lengths. The free energy of any sequence is however the sum of loop energies, and each loop has a unique maximal arc. In Jin and Reidys [52] it is shown that the number of arcs in random structures satisfies a central limit theorem, whence its mean scales linearly with n. This implies that the number of loops grows linearly with n, which in turn suggests that the free energy of a sequence grows linearly with n.
Accordingly, we consider the scaled partition function and set We call w R and w S the densities of the structure R and S respectively. The adaptability w R is a real number in [0, 1] measuring the proportion of the partition function of R composed by bicompatible sequences relative to that . composed by sequences compatible only with R. The closer this adaptability is to 1, the more energy favorable bicompatible sequences there are, suggesting that the structure R can change into the structure S more easily. Note that by construction w R and w S are asymmetric, namely, the transitions from R to S and those from S to R are not necessarily equal.

Results
In this section we focus on using the RNA sequence sampler to study the sequence-structure relations of RNA riboswitch sequences. The time complexity of the algorithm presented here is O(4 κ(B) n) , where κ(B) is constant depending only on the bistructure B, and n is the length of the two input sequences. Currently, we can only deal reliably in case of κ(B) ≤ 20 . In the following, we select six riboswitch sequences from the literature, across different classes, organisms and switching mechanisms, exhibiting a sequence length below 250nts, and κ(B) ≤ 20 , see Table 1.

Energy spectra
We are now in position to study the R-and S-spectra of Q-Boltzmann sampled, bicompatible sequences of specific structure pairs, or equivalently the R-and S-spectra of sequences compatible to bistructures using the measure of "Energy spectra" section. It will be interesting to compare these with the spectra of the compatible sequences of each respective secondary structures, R and S and to provide a comparative analysis of the R-and S-spectra of native structures with that of random structure pairs. Let R, S be two secondary structures, to begin the comparative analysis of compatible and bicompatible sequences we compare Q(R) and Q(S), given by with the partition functions Q(R) and Q(S) are recursively computed and their Boltzmann samplers are introduced in [32,33]. On an abstract level, Q(R)| S and Q(S)| R were a priori available by means of rejection Boltzmann samplers for Q(R) and Q(S). However, sampling such sequences is impractical as the probability of randomly encountering a bicompatible sequence is too low. As a first application, our framework developed in "Materials and methods" section we observe that Q| S (R) and Q| R (S) can be computed by replacing the energy function η(σ , B(R, S) by η(σ , R) and η(σ , S) , respectively.
Comparing the R-and S-spectra of Q(R) with Q(R)| S and Q(S) with Q(S)| R , respectively, allows us to draw conclusions about the difficulty of the process of modifying a R-compatible sequence into an R, S-bicompatible sequence, while maintaining the energy with respect to R and the analogue statement for S.
In Fig. 11d we compare f R Q(R) with f R Q(R)| S , i.e. the R-energy spectra over Q(R) and Q(R)| S (LHS) and f S Q(S) with f S Q(S)| R , i.e. the S-energy spectra over Q(S) and Q(S)| R (RHS). We remark that these are pairwise comparisons of measures over distinctively different, nested probability spaces.
We show that the R-energy spectra over Q(R) and Q(R)| S and S-energy spectra over Q(S) and Q(S)| R pairwise coincide. This means, that modifying a R-compatible sequence into a R, S-bicompatible sequence can be done without affecting the free energy with respect to R and vice versa for S. Moreover, this finding holds for all native, as well as random structure pairs we analyzed. The results for the other riboswitch as well as random structures pairs are shown in the SM. The next step is to relate bicompatible sequences, Boltzmann sampled from Q(R, S) to those Boltzmann sampled from Q(R). In the context of the Q(R)| S versus  16:7 Q(R) analysis, we now factor in the energy with respect to S. While R-energy levels can be maintained while satisfying the S base-pairing conditions, it turns out to be much more intricate to derive bicompatible sequences that are well suited for both R and S at the same time. Thus we consider and compare the energy spectra f R Q(R,S) and f R Q(R) as well as f S Q(R,S) and f S Q(S) for a variety of riboswitch sequences and their respective native structure pairs, see Table 1. We display the energy spectrum of the riboswitch add in Fig. 12(top). The energy spectra and detailed analysis of the other riboswitch sequences is presented in the SM.
In order to put our results into context, we present an analysis of the spectra f R Q(R,S) and f R Q(R) as well as f S Q(R,S) and f S Q(S) for structure pairs that are not riboswitches. We consider two random RNA sequences and fold them to their mfe-structures, where each structure is thermodynamically stable to not only one of the sequences, but to both. For practical reasons, we consider sequences of length 150, which is close to the average length of natural riboswitches. We show a Q(R, S) = σ ∈C n (R,S) e −1/2·(η(σ ,R)+η(σ ,S)) KT representative result for in Fig. 12(bottom). Our observations are remarkably robust: the energy spectra of random structure pairs are literally identical and we provide a detailed analysis of additional spectra in the SM. Figure 12(top left) shows that the R-spectrum, f R Q(R,S) , is practically identical to the R-spectrum, f R Q(R) . The S-spectra behave completely analogous, there is no significant difference between the S-spectrum f S Q(R,S) and f S Q(S) , see Fig. 12(top right). In the SM we provide additional data and the spectra of the riboswitch sequences listed in Table 1. The phenomenon holds robustly for all riboswitch sequences we analyzed.
For random structure pairs, however, the picture changes: the R-spectrum, f R Q(R,S) , is shifted distinctively to the left of the R-spectrum, f R Q(R) , see Fig. 12(bottom left). The same holds for the S-spectra: f S Q(R,S) is shifted distinctively to the left of the S-spectrum f S Q(S) , see Fig. 12(bottom right). More data for random structure pairs are presented in the SM. The different patterns of the R-spectrum and S-spectrum for riboswitches in contrast to random structure pairs reveals a strong signal as the two structures corresponding to a riboswitch are strongly correlated. As such, the signal can be used to identify riboswitch features.

Ranking
We compute the (r R , r S ) for the riboswitches mht and lys according to the test introduced in "Ranking". We display our result in Fig. 13. Figure 13 displays the ratios (r R , r S ) for the riboswitch mgt (LHS) and lys (RHS). The figure is obtained  For each riboswitch we Boltzmann sample 10 3 sequences and compute (r R , r S ) for the sampled sequences (blue). We contrast this with (r R , r S ) of for the respective, native riboswitch sequence (red) based on the Boltzmann sampling of 10 3 sequences from (Q(R, S)) (blue) and displays in addition the ratios of the native sequences of mgt and lys (red). We find for mgt that all ratios satisfy r R > 70% and r S > 62% and furthermore the respective (coordinatewise) means are (85%, 74%) . For lys we have r R > 64% and r S > 72% with a mean of (78%, 85%) . Accordingly, R and S are suboptimal structures in the Boltzmann sampled sequences. Furthermore, we find that the ratio of ratios, r R /r S is almost constant within the set of sampled sequences. For both, mgt and lys alike, (r R , r S ) of the native sequence is distinctively higher than the ratio pairs obtained from the sampled sequences. This indicates that the native sequence exhibits an evolved thermodynamic stability with respect to the pair of structures (R, S).

Density
We compute (w R , w S ) for all riboswitch structures pairs in Table 1 according to measure introduced in "Adaptability" section, and augment the analysis by inspecting 50 random structure pairs as control set. The random structure pairs are obtained using the method described in "Energy spectra" section. We consider random structure pairs having length of 150, which is close to the average sequence length of riboswitches. We display in Fig. 14 the pairs (w R , w S ). Figure 14 shows that the pairs (w R , w S ) of riboswitch sequences are distinctively different from random structure pairs and appear in the top right corner, while the (w R , w S ) of random structure pairs are shifted towards the bottom left corner. The closeness of ratio pairs displayed in Fig. 14 to the top right corner represents how likely a sequence, Boltzmann sampled from Q (R) , is contained in Q (R)| S . This reflects how dense the bicompatible sequences sampled from Q (R) are within the compatible sequence Boltzmann sampled from Q (R) . Figure 14 shows that this adaptability is significantly higher for native riboswitch sequences, compared to sequences Boltzmann sampled from Q (R)| S for random structure pairs (R, S).

Discussion
The time complexity of computing the partition function for a given pair of secondary structures is determined by the decomposition of the bistructure and the computation of the partition function. Hammer et al. [29] shows that when assuming a tree decomposition is given with a constant tree-width, the partition function can be computed in polynomial time. However, obtaining such a tree decomposition with minimum tree-width remains an NP-hard problem. The work in [29] focuses on the former, presenting an algorithm that computes the partition function on a simplified loop-based energy model. Although in general a hyper-graph can be decomposed into a tree-like structure having tree-width k, the parameter k is not understood in a systematic way.  Table 1 (blue). Furthermore: (w R , w S ) of 50 random structure pairs of length 150, as control set (red)