# Designing RNA structures while maximizing stacking interactions

Théo Boury, Laurent Bulteau, Yann Ponty

<figure><img src="image_1.png" /></figure>

## To cite this version:

Théo Boury, Laurent Bulteau, Yann Ponty. Designing RNA structures while maximizing stacking interactions. 2026. hal-05510148v2

HAL Id: hal-05510148

https://hal.science/hal-05510148v2

Preprint submitted on 16 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Designing RNA structures while maximizing stacking interactions

Théo Boury

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

LIX, CNRS, École polytechnique, Institut Polytechnique de Paris, Palaiseau, France

Laurent Bulteau

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

LIGM, Université Gustave Eiffel, Marne-la-Vallée, France

Yann Ponty

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

LIX, CNRS, École polytechnique, Institut Polytechnique de Paris, Palaiseau, France

## Abstract

Inverse folding is a classic problem in RNA bioinformatics, crucial for the design of functional synthetic RNAs, and consists in finding a sequence that uniquely folds into a target secondary structure with respect to energy minimization. In a simple base pairs maximization model (maxBPs), Bonnet et al. showed that a mildly constrained version of inverse folding is NP-hard. By contrast, a linear-time exact algorithm was proposed for maxBPs inverse folding, when restricted to input structures with helices of size 3<sup>+</sup>. However, the maxBPs model artificially induces drastic limitations on the set of designable structures, forbidding the design of well-known RNA families.

In this work, we consider a more realistic energy model accounting for base pairs stacks, a major contributor to RNA stability, and study inverse folding in a stacks maximization model (maxStacks). We propose a linear-time algorithm for maxStacks inverse folding, restricted to helices consisting of at least ⌈log<sub>3.56</sub>(Δ) + 6.2⌋ base pairs, where Δ is the largest degree of a loop in the target structure. Our approach relies on the introduction of the locked property, a concept which in particular induces a unique matching of intervals involved in helices, and represents a sufficient condition for a sequence to represent a solution of maxStacks inverse folding. In particular, our algorithm implies the designability of structures featuring loops of arbitrary degree Δ in the maxStacks model, something that was intrinsically impossible in the maxBPs model beyond Δ > 4.

Moreover, the locked property can be generalized to target structures and competing structure which include crossing base pairs (aka called general pseudoknots). In this challenging setting, we obtain a linear-time algorithm for maxStacks inverse folding, restricted to structures with helix sizes at least ⌈log<sub>3.56</sub>(*m*) + 6.2⌋, where *m* is the number of helices. This result is surprising since validating designs requires solving the RNA folding with general pseudoknots, a problem known to be NP-hard in the maxStacks model.

2012 ACM Subject Classification Applied computing Molecular structural biology

Keywords and phrases RNA structure, String Design, String combinatorics

Digital Object Identifier 10.4230/LIPIcs.CPM.2026.23

<figure><img src="image_9.png" /></figure>

© Théo Boury, Laurent Bulteau and Yann Ponty;
licensed under Creative Commons License CC-BY 4.0

37th Annual Symposium on Combinatorial Pattern Matching (CPM 2026).

Editors: Philip Bille and Nicola Prezza; Article No. 23; pp. 23:1–23:17

<figure><img src="image_10.png" /></figure>

LIPICS

Leibniz International Proceedings in Informatics
Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl Publishing, Germany
1 Introduction

RiboNucleic Acids (RNAs) represent a family of biomolecules equipped with a dual capacity to encode genetic information and perform complex catalytic functions [6]. Such a versatility makes RNA a promising substrate for the fast-paced development of therapeutic tools, including mRNA vaccines, RNA interference or molecular sensors (riboswitches), fueled by a diversity of cellular mechanisms [8, 15]. Such functions often require the adoption (or avoidance) of intricate 3D structures, motivating a rational structural design of RNAs, fueled by the availability of robust predictive methods for RNA folding [9]. Inverse folding is crucial ingredient in the RNA design cookbook, enabling a control of the structure likely adopted by RNA at its thermodynamic equilibrium [9].

In its minimal version [7], called base pair maximization (maxBPs), the RNA folding problem asks, given an input nucleotides sequence $w \in \Sigma^*$, $\Sigma := \{A, C, G, U\}$, for a maximum cardinality subset of, pairwise non-crossing, base pairs in $B := \{\{A, U\}, \{G, C\}, \{G, U\}\}$. RNA folding in the maxBPs model is solvable in $O(n^3)$ using dynamic programming [12]. Given a target structure $T$, RNA inverse folding asks for a sequence whose maxBPs-folding is unique and coincides with a specified target structure [7]. Supplemented with minor additional constraints, inverse folding in the maxBPs model is a NP-hard problem for non-crossing targets [2]. It motivates practical heuristics [5, 13] targeting more complex energy models [17]. However, this complexity result critically focuses on target structures unrealistically intricate, possibly misrepresenting the practical complexity of inverse folding in realistic settings. For instance, restricting the input target $T$ to helices with $3^+$ base pairs enables a $O(n)$ solution for inverse folding in the maxBPs model [3].

The maxBPs model represents a useful simplification, amenable to combinatorial studies [7, 16] and initialization strategies [4], but unfortunately comes with intrinsic limitations. Firstly, it postulates independent energy contributions for base pairs and, as such, disregards more complex structural motifs such as loops [17]. Loops consist of the interaction of multiple base pairs, and are known to majorly contribute to structural stability. Worse, the maxBPs model induces a combinatorial explosion of alternative structures that are equally favorable as the target, leading to an absence of solution for structures having at least one loop involving unpaired positions, or featuring more than 4 helices [7]. This artificially deflates the set of designable structures [16], and does not accurately model the RNA folding process, since such motifs are crucial to many RNAs occurring in nature, including tRNAs, rRNAs...

The limitations of the maxBPs model for inverse folding motivate the consideration of more realistic energy models. Unfortunately, a rigorous algorithmic study of the full Turner (nearest-neighbor) model [17] remains elusive due to its large number of experimentally-measured parameters (several thousands, although not always independent). We argue that the next natural step is to consider a folding model based on maximizing the number of base pair stackings $(i, j)|(i+1:j-1)$ (maxStacks). The folding problem remains solvable in $O(n^3)$ time using dynamic programming, yet maximizing the number of stacks does not immediately reduce to maximizing the number of base pairs. In particular, when crossing interactions (aka general pseudoknots) are allowed, maxBPs folding reduces to the polynomial-time solvable MAX MATCHING, while maxStacks folding is NP-hard [10].

Given a target structure $T$ of length $n$, maxStacks inverse folding asks for a sequence $w^* \in \Sigma^n$, valid for $T (\forall (i,j) \in T: \{w_i, w_j\} \in B)$, such that the number of stacks $(i,j)|(i+1,j-1)$ in $T$ strictly exceeds that of any valid alternative $T' \neq T$ for $w^*$. In this article, we study inverse folding in the maxStacks model, and find that:

■ The scope of maxStacks expands beyond that of the maxBPs model, enabling a design
82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

of structures having arbitrary degrees loops, even in the presence of unpaired positions;

- When the target and space of alternative structures are both restricted to avoid crossing interactions/pseudoknots, we find a $O(n)$ solution for the maxStacks inverse folding, restricted to structures with helices having at least $\lceil \log_{3.56}(\Delta) + 6.2 \rceil$ base pairs, $\Delta$ being the largest degree of a loop;

- When arbitrary crossing interactions/pseudoknots are allowed both in the target and alternatives, we provide a linear-time design algorithm for inverse folding in the maxStacks model, for any structure having helices of length $\lceil \log_{3.56}(m) + 6.2 \rceil$, $m$ being the number of helices in the target. Notably, our algorithm guarantees the validity of the solution despite the NP-hardness of the maxStacks general folding problem.

These results rely on the introduction of a sufficient condition, introduced in Section 4, for a nucleotide sequence to be a solution to the inverse folding problem. Overall, our result shows that the stacking model represents a reasonable tradeoff between accuracy and complexity in the context of design, and that exact strategies can be used to address a variety of relevant use cases. Following a formal definition of inverse folding and connected notions in Section 2, we further motivate our study of maxStacks in Section 3. We then introduce in Section 4 a sufficient condition for RNA sequences to represent solutions to the maxStacks inverse folding, and address the algorithmic search for such sequence in Section 5. Finally, we summarize and discuss the scope of our contribution in Section 6.

101

## 2 Problems statement and definitions

102

103

104

105

106

107

As stated above, an RNA sequence is a string $\omega \in \Sigma^n$, $\Sigma := \{A, C, G, U\}$. An RNA structure $T$ consists of a set of base pairs, i.e. pairs of integers $(i, j) \in [1, n]^2$, $j > i+1$ that are pairwise disjoint. A structure $T$ is compatible with a sequence $\omega$ iff $\forall (i, j) \in S : \{\omega[i], \omega[j]\} \in \mathcal{B}$ with $\mathcal{B} := \{\{A, U\}, \{G, C\}, \{G, U\}\}$. A structure is crossing (aka pseudoknotted), if there exists two base pairs $(i, j)$, $(k, l) \in T$ such that $i < k < j < l$, and non-crossing otherwise. Structures that are free of crossing base pairs are called a secondary structure (2D structure).

108

109

110

As illustrated in Figure 1, a structure $T$ can be uniquely decomposed into a set $\mathcal{H}_T$ of helices, each helix $H = \{(i, j), (i+1, j-1), \dots, (i+h-1, j-h+1)\} \subset T$ being a maximal subset of consecutive base pairs, with $h = |H|$ denoting the size of helix $H$.

111

112

113

114

A stacking base pair, or stack in short, is a pair of consecutive base pairs $(i, j), (i+1, j-1)$ occurring in a structure $T$, and we denote by $S_T$ the set of stacks induced by $T$. Note that an helix of size $h$ consists of $h-1$ stacks, so we have $|S_T| = \sum_{H \in \mathcal{H}_T} |H| - 1$. A structure is *admissible* if each base pair belongs to at least one stack.

115

**Problem 1. (MAXSTACKS) PK-DESIGN**

116

**Input:** Admissible target structure $T$; length $n$.

117

**Output:** RNA sequence $\omega \in \Sigma^n$, compatible with $T$, such that if $T' \neq T$ is an admissible structure compatible with $\omega$ then $|S_T| > |S_{T'}|$; or $\emptyset$ if no such sequence exists.

118

119

**Non-crossing/secondary structure (2D) vs crossing/PseudoKnotted (PK) versions**

120

121

122

123

124

The PK version of the problem defined above does not enforce any particular restriction on structure, so that admissible structures are allowed (but not forced) to be pseudoknotted. We also consider the 2D-DESIGN variant, where *admissible* structures are necessarily non-crossing. Figure 2 shows two sequences, compatible with a target, that admit a competing structure, or represent a solution for the 2D version of the inverse folding problem.

125

126

Note that this is not a mere restriction of the input, since a no-instance in the PK model can become a yes-instance in the 2D model, where any design would only have pseudoknotted
23:4 Designing RNA structures while maximizing stacking interactions

<figure><img src="image_11.png" /><figcaption>■ **Figure 1** (A) Non-crossing (2D) structure, consisting of 5 helices and maximum loop degree 4 with {$H_1, H_2, H_3H_4$}; (B) Associated set of regions; (C) Classic representation of the structure.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>■ **Figure 2** A target structure *T* with 3 helices for a total of 7 stacks, and two sequences $\omega_1, \omega_2$ compatible with *T*. Sequence $\omega_1$ is not a design since there exists an alternative structure $T'$ with the same number of stacks. On the other hand, sequence $\omega_2$ is a solution for the 2D version of PK-DESIGN(*T*).</figcaption></figure>

alternatives. An important difference lies in the problem of deciding whether a given sequence is actually a maxStacks design for some structure, or equivalently building optimal structures for a given sequence. This verification problem can be solved in $O(n^3)$ time in the 2D setting using dynamic programming [7]. On the other hand, verification becomes NP-hard whenever arbitrary pseudoknots are allowed [10]. So the PK-DESIGN problem may not even belong in NP. Here we consider both models for the definition of locked sequences, and our algorithm is written to adapt to either model.

**Further definitions: match, regions and loops**

For any interval I, we write $I_{\text{start}}$ and $I_{\text{end}}$ respectively for its first and last index of I (i.e. $I = [I_{\text{start}}, I_{\text{end}}]$). The complement of x = A (resp U, C, G) is $\bar{x}$ = U (resp. A, G, C). The reverse complement of string $S = x_1...x_\ell$ is $\text{RevComp}(S) = \overline{x_\ell...x_1}$. For an helix $H = \{(i,j),... (i+h-1, j-h+1)\}$, intervals $r = [i, i+h-1]$ and $r' = [j-h+1, j]$ are the regions of $H$, and we write $r' = \bar{r}$ and $r = \overline{r'}$. We write $\mathcal{R}_T$, or simply $\mathcal{R}$ whenever clear, for the overall set of (pairwise disjoint) regions occurring in $T$. Since an helix is fully determined by its two regions, we slightly abuse notations and write $H = (r, \bar{r})$ for the helix with regions $r$ and $\bar{r}$ (in the absence of further precision, we choose $r$ as the leftmost interval among both regions of $H$). See Figure 1 for an illustration of those definitions

Given a target *T* and a sequence *ω*, we say that an interval *I* = [*a*, *a* + *h*] *matches* (resp. *aligns with*) interval [*b* − *h*, *b*] in *T* if for all *k* ∈ [*0*, *h* − 1], (*ω*[*a* + *k*], *ω*[*b* − *k*]) ∈ {{*A*, *U*}, {*G*, *C*}} (resp. ∈ *B*). Finally, we use *ω*[*r*] := (*ω*[*a*], ..., *ω*[*b*]) to denote the substring associated with a region *r* = [*a*, *b*]. With such notations, the compatibility criterion becomes: for each *r* ∈ $\mathcal{R}$, *r*
aligns with $\bar{r}$.

Let $r, r' \in \mathcal{R}_T$. They are *adjacent* if $r'_{\text{start}} = r_{\text{end}} + 1$. Also, $r'$ is the *successor* of $r$ (and $r$ the *predecessor* of $r'$) if $r'$ is the first region to the right of $\bar{r}$. Each region can have at most one successor and predecessor, so $\mathcal{R}_T$ can be partitioned into paths and cycles by this relation. All helices with a region in the same component are called a *loop*, so that each helix $(r, \bar{r})$ belongs to at most two loops (one for $r$, one for $\bar{r}$). The size (number of helices) in a loop is called its *degree*. This definition holds in both PK and 2D setting, but it is mostly used in 2D thanks to a useful property: given an helix $(r, \bar{r})$, we can define an inside $[r_{\text{start}}, \bar{r}_{\text{end}}]$ and an outside $([1, r_{\text{end}}] \cup [\bar{r}_{\text{start}}, n])$, and all helices and loops are included in either side by the non-crossing property.

Due to space constraints, we delegate some proofs of our results to an appendix.

# 3 Why study inverse folding in the maxStacks model?

## Loops of arbitrary degrees can be designed in the maxStacks model

Previous work [7] showed that the maxBPs model fails to produce designs for any target structure with either: i) a loop of degree 5$^+$; or ii) a loop of degree 3$^+$, including 1$^+$ unpaired position [7]. This limitation is largely detached from reality, as such motifs can be found in virtually any RNA structure, including tRNAs, 5s ribosomal RNA...

Contrasting, loops of arbitrarily large degrees can be designed for maxStacks using sufficiently long helices, as illustrated by Figure 3. Such an observation generalizes to any secondary structure, possibly including pseudoknots/crossing interactions as long as sufficiently long helices are considered. However, all helices in the same loop need to be designed with distinct sequences, leading to the following logarithmic lower-bound on the helix size.

**Lemma 1.** Let $\mathcal{L}$ be the 2D structure consisting of a single loop over $\Delta$ helices of size $h$. If $h < \log_6(\Delta)$ then $\mathcal{L}$ is not designable in maxStacks.

## The maxStacks model mimics the Turner model better than maxBPs

We further empirically assess the potential of designs in the maxStacks to produce promising candidates for the Turner energy model [17]. We focus on *undesignable loops*, having either degree 5$^+$, or degree 3$^+$ with 1$^+$ unpaired position, since those are not designable in the maxBPs model. We create secondary structures instances of length $n = 150$ using random uniform generation [3] (min BP distance $\theta = 3$), coupled with a rejection step to generate two sets: *small loops* (0 undesignable loop) and *large loops* (at least one undesignable loop), each consisting of 2000 secondary structures. We then use polynomial-time algorithms [7, 9] to check whether those sequences also represent designs in the maxStacks and the Turner models<fnref n="1" /> (see [3] for details).

For the small loops dataset, we produce solutions for the maxBPs model using uniform random generation of *m*-separated sequences [3]. Note that *m*-separated sequences represent a strict subset of possible solutions in the maxBPs model, which may induce a statistical bias. However, a direct uniform generation of maxBPs solutions is probably unattainable, since finding a single maxBPs design with fixed positions is NP-hard [2], and so is the problem of finding separated sequences, while a uniform random generation algorithm exists

<footnote n="1">¹ Source code is available at https://github.com/theoboury/positioningempiricalmaxstacksdesigns</footnote>
<figure><img src="image_13.png" /><figcaption>Figure 3 Loops of arbitrary degree $\Delta$ can be designed in the maxStacks model using helices of size at least $1 + \lfloor \log \Delta \rfloor$. A solution sequence can be found for a target with loop of degree $\Delta = 5$, and helix size $h = 4$, a structure which cannot be designed in the maxBPs model.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 4 Assessment of the potential of maxStacks solutions as solutions to inverse folding in the Turner energy model. For structures restricted to small loops, maxStacks and maxBPs designs are essentially equivalently promising, featuring ~15% designs in the Turner model. However, for target structures with larger loops, i.e. no solution for maxBPs, 6.5% maxStacks solutions still represent valid designs in the Turner model (vs 0.6% amongst compatible sequences).</figcaption></figure>

for *m*-separated sequences [3]. We compare these maxBPs solutions to those obtained for the maxStacks model, produced by rejection from a set of random uniform compatible sequences (A assigned to unpaired positions). As shown in Figure 4, both maxBPs and maxStacks are equally likely to represent admissible solutions for the Turner model, with 15% of sequences folding uniquely into their respective target [3]. However, 85% of maxBPs designs are maxStacks, while none of the maxStacks designs are solutions to maxBPs. This suggests that, while both models are promising towards inverse folding in the Turner model, maxStacks is more realistic and may ease the resolution of inverse folding.

In the large loops dataset, separated sequences cannot be generated as no solution exists in the maxBPs model. We thus use as baseline a set of sequences compatible with their respective target, and compare those sequences to those obtained for the maxStacks model. As shown in Figure 4, the maxStacks energy model not only enables a designs of structures containing large loops, but also allows produces a substantially larger proportion of designs in the Turner energy model: 6.45% of maxStacks designs directly represent Turner designs, against a mere 0.3% amongst compatible sequences. This motivates the upcoming design of algorithmic solutions for the production of solutions to the maxStacks inverse folding.

4 Sufficient conditions for inverse folding in the maxStacks model

In this section we present the *locked* criterion for sequences and prove that satisfying this criterion is sufficient to yield a design. This criterion is based on the presence of a *key* within each region, with strong constraints on how keys can be matched within the rest of the sequence, thus forbidding “global” rearrangements where whole regions are entirely mapped differently. The non-extension criterion, on the other hand, is crucial to forbid “local” rearrangements that would locally increase an existing helix by a few base pairs. We first give a general criterion that applies to any restriction on the feasible structures in Definition 1, prove its core property in Lemma 2. We then specify the criterion into 2D-locked and PK-locked, leading finally to Theorem 3.
<figure><img src="image_15.png" /><figcaption>Figure 5 Illustration for the proof of Lemma 2. Top: a locked sequence for the target structure T shown with one interval per region, and a darker stripe for the key. Bottom: an alternative structure T'. In order to show that T' has fewer stacks than T, we carefully count helices and base pairs in T', with a case distinction for each region of T depending on its overlap with regions of T'.</figcaption></figure>

**Definition 1 (locked).** Given a sequence $\omega$ of $T$, we say that $\omega$ satisfies the locked property if there exists an interval $\lambda^r$ for each $r \in \mathbb{R}$, named the key of $r$, such that:

1. **Base compatibility:** $\omega$ is compatible with $T$, using only $\{A, U\}$ and $\{G, C\}$ base pairs;

2. **Unpaired A or C:** All unpaired positions are assigned A or C;

3. **Helix keys:** $\forall r \in \mathbb{R}$, each interval $[r_{\text{start}}, \lambda_{\text{start}}^r]$ and $[\lambda_{\text{end}}^r, r_{\text{end}}]$ contains at least one position $i$ with $\omega[i] \in \{G, U\}$;

4. **Keys match keys:** $\forall r \in \mathbb{R}$, if $\lambda^r$ matches any interval $I$, then $\exists r' \in \mathbb{R}, I = \lambda^{r'}$;

5. **No Interval Extension:** $\forall r \in \mathbb{R}$, if an interval $I$ with $\lambda^r \subseteq I$ and $I \not\subseteq r$ aligns with any other interval, then $r$ contains at least one position in a $\{G, U\}$ pair.

**Lemma 2 (Key Matching Lemma).** Let $T$ be a structure and $\omega$ be a locked sequence for $T$. For any $T'$ compatible with $\omega$, if $|S_T| \le |S_{T'}|$, then all positions labeled G or U are paired in $T'$ (without any $\{G, U\}$ pair), and each key $\lambda^r$ of $\omega$ matches $\lambda^{r'}$ in $T'$ with some $r' \in \mathbb{R}_T$.

**Proof.** The overall proof strategy is depicted in Figure 5. We write $\rho = |\mathbb{R}_T|$ for the number of regions of $T$, $\pi_{GU}$ for the number of positions belonging to a pair $\{G, U\}$ in $T$, and $\mu$ (resp. $\mu_{GU}, \mu_{AC}$) for the number of unpaired positions in $T$ (resp unpaired positions with label in $\{G, U\}$ or $\{A, C\}$ in $\omega$), i.e. $\mu = \mu_{GU} + \mu_{AC}$. We use the corresponding primed notations for $T'$ ($\rho', \pi'_{GU}, \mu', \mu'_{GU}$ and $\mu'_{AC}$). We first show the following:

$$|S_T| - |S_{T'}| = \frac{\mu' - \mu + \rho' - \rho}{2} \qquad (1)$$

Indeed, there are $n - \mu$ paired positions in $T$, thus accounting for $\frac{n-\mu}{2}$ base pairs. Each size-h helix contributes to $h-1$ stacks, so the overall number of stacks is the number of base pairs minus the number of helices. There are $\frac{\rho}{2}$ helices in $T$, so $|S_T| = \frac{n-\mu-\rho}{2}$ and similarly $|S_{T'}| = \frac{n-\mu'-\rho'}{2}$, thus yielding Equation (1).

We first compare numbers of unpaired positions with label A or C in $T$ and $T'$. Bases A or C can only be paired with G and U. Writing $n_{GU}$ for the total number of Gs and Us in $\omega$, there are $n_{GU} - \pi_{GU} - \mu_{GU}$ pairs in $\{\{A, U\}, \{C, G\}\}$. Thus, the difference in number of unpaired A and C is: $\mu'_{AC} - \mu_{AC} = (\pi'_{GU} + \mu'_{GU}) - (\pi_{GU} + \mu_{GU})$. By the locked property, $\pi_{GU} = 0$ (Condition 1), $\mu_{GU} = 0$ (Condition 2), so $\mu'_{AC} - \mu_{AC} = \pi'_{GU} + \mu'_{GU}$. Overall, Equation (1) can then be improved into the following:

$$|S_T| - |S_{T'}| = \frac{\mu'_{AC} + \mu'_{GU} - \mu_{AC} - \mu_{GU} + \rho' - \rho}{2} = \frac{\pi'_{GU} + 2\mu'_{GU} + \rho' - \rho}{2} \qquad (2)$$
## 23:8 Designing RNA structures while maximizing stacking interactions

We then focus on regions of $T$ and $T'$. We partition regions of $T$ into the following categories, so that $\rho = \rho_{\simeq} + \rho_{\subset} + \rho_2 + \rho_1 + \rho_0$:

- $\rho_{\simeq} \to \textit{conserved}$ regions $r \in \mathcal{R}_T$ for which there exists $r' \in \mathcal{R}_{T'}$ with $\lambda^r \subseteq r' \subseteq r$

- $\rho_{\subset} \to \textit{extended}$ regions $r \in \mathcal{R}_T$ for which there exists $r' \in \mathcal{R}_{T'}$ with $\lambda^r \subset r'$ and $r' \not\subseteq r$

- $\rho_2 \to \textit{2-split}$ regions $r \in \mathcal{R}_T$ that are neither conserved nor extended, and that intersect at least two distinct regions in $\mathcal{R}_{T'}$.

- $\rho_1 \to \textit{1-split}$ regions $r \in \mathcal{R}_T$ that are neither conserved nor extended, and that intersect exactly one region in $\mathcal{R}_{T'}$

- $\rho_0 \to \textit{lost}$ regions $r \in \mathcal{R}_T$ that do not intersect any region in $\mathcal{R}_{T'}$

Further, a 2-split or extended region is a border if it is the first region in $T$ or if the region directly to its left is conserved, we write $\beta$ for the number of border regions. There are overall $2\rho'$ region extremities in $T'$, with at least 2 extremities in conserved and 2-split regions of $T$, and at least one extremity in 1-split regions. Plus, for a border region $r$, the left extremity of the left-most region of $T'$ intersecting $r$ is not counted in the above (it can only fall in a conserved region by definition of border, and it would count as an additional extremity in such a conserved region).

$$2\rho' \geq 2\rho_{\simeq} + 2\rho_2 + \rho_1 + \beta. \qquad (3)$$

Consider now a 1-split or lost region $r$. By Condition 3, both intervals $[r_{\text{start}}, \lambda^r_{\text{start}}]$ and $[\lambda^r_{\text{end}}, r_{\text{end}}]$ contain a position $i$ with $\omega[i] \in \{\text{G}, \text{U}\}$. If $r$ is lost, then it contains at least two unpaired positions labeled $\{\text{G}, \text{U}\}$. If $r$ is 1-split and intersecting $r' \in \mathcal{R}_{T'}$, then $r'$ cannot intersect both intervals (otherwise it would also contain $\lambda^r$ and $r$ would be conserved or extended), so $r$ features at least one unpaired $\{\text{G}, \text{U}\}$. Thus

$$\mu'_{\text{GU}} \geq 2\rho_0 + \rho_1. \qquad (4)$$

Finally by Condition 5, if $r$ is extended with $r'$, then it features at least one endpoint of a $\{\text{G}, \text{U}\}$-pair in $T'$, i.e.

$$\pi'_{\text{GU}} \geq \rho_{\subset} \qquad (5)$$

Combining the inequalities (3), (4) and (5) with Equation (2), we get:

$$\begin{aligned}
|S_T| - |S_{T'}| &= \frac{1}{2}(\pi'_{\text{GU}} + 2\mu'_{\text{GU}} + \rho' - \rho) \\
&\geq \frac{1}{2}(\pi'_{\text{GU}} + 2\mu'_{\text{GU}} + (\rho_{\simeq} + \rho_2 + \frac{\rho_1}{2} + \frac{\beta}{2}) - (\rho_{\subset} + \rho_{\simeq} + \rho_2 + \rho_1 + \rho_0)) \\
&\geq \frac{1}{2}(\pi'_{\text{GU}} + 2\mu'_{\text{GU}} + \frac{\beta}{2} - (\rho_{\subset} + \frac{1}{2}\rho_1 + \rho_0)) \\
&\geq \frac{1}{2}(\pi'_{\text{GU}} + 2\mu'_{\text{GU}} + \frac{\beta}{2} - (\pi'_{\text{GU}} + \frac{1}{2}\mu'_{\text{GU}})) \\
&\geq \frac{3}{4}\mu'_{\text{GU}} + \frac{\beta}{2}
\end{aligned}
\qquad (6)$$

We now use $|S_{T'}| \geq |S_T|$ so $\mu'_{\text{GU}} = 0$ and $\beta = 0$. This implies $\rho_1 = \rho_0 = 0$ with inequality (4), so all regions of $T$ are either conserved, 2-split or extended. Moreover, $\beta = 0$: there can be no 2-split or extended region (otherwise the leftmost 2-split or extended region would be a border), and all regions of $T$ are conserved. Thus each $\lambda^r$ is included in some region of $T'$ in a one-to-one mapping, and $\rho' = \rho$. Using Equation (2), we obtain $\pi'_{\text{GU}} = 0$. So for any $r \in \mathcal{R}_T$, $T'$ matches (without $\{\text{G}, \text{U}\}$ pairs) $\lambda^r$ with an interval $I$, and by Condition 4, $I = \lambda^{r''}$ for some $r'' \in \mathcal{R}_T$.

<figure><img src="image_16.png" /></figure>
**Definition 2.** A sequence $\omega$ is PK-locked (resp. 2D-locked) for structure $T$ if $\omega$ is locked, with the additional condition **4a** (resp. **4b**) below, that overrides original condition **4**:

4a **Unique Keys:** $\forall r \in \mathcal{R}, \exists \lambda' \neq \lambda^{\bar{r}}$ such that $\lambda^r$ matches $\lambda$.

4b **Loop Unique Keys:** $\forall r \in \mathcal{R}, \forall \lambda' \neq \lambda^{\bar{r}}$ such that $\lambda^r$ matches $\lambda'$ then $\exists r' \in \mathcal{R}, \lambda' = \lambda^{r'}$ and $r, r'$ are not in the same loop.

**Theorem 3.** Given a structure $T$, if a sequence $\omega$ is PK-locked (resp. 2D-locked), then $\omega$ is a solution of PK-DESIGN($T$) (resp. 2D-DESIGN($T$)).

**Proof.** We first show the following claim:

► Claim 4. For any alternative $T'$ with $|S'_T| \ge |S_T|$ and any helix $(r, \bar{r})$, $T'$ matches $\lambda^r$ with $\lambda^{\bar{r}}$.

**Proof.** In the PK case, this is a direct consequence of Lemma 2: any alternative must match keys with their reverse complement in the sequence, but by Condition **4a** each reverse complement appears only once and in the desired region. In the 2D case, we consider a *wrong match* as a pair of keys from different loops $\lambda^r$, $\lambda^{r'}$ with $r$ and $r'$ from distinct helices (from distinct loops by with condition **4b**), with $\lambda^r$ and $\lambda^{r'}$ matched together in $T'$. We consider a *deepest* wrong match, such that no other wrong match occurs in the inside of the helix formed with this match in $T'$. Due to $\lambda^r$, $\lambda^{r'}$ being from different loops in $T$, their inside $I = [\lambda^r_{start}, \lambda^r_{end}]$ contain another key $\lambda^x$, whose match $\lambda^{\bar{x}}$ (in $T$) is not in $I$. By Lemma 2, $\lambda^x$ matches perfectly another key in $I$, thus different from $\lambda^{\bar{x}}$, and $\lambda^x$ also forms a wrong match: a contradiction with $\lambda^r$, $\lambda^{r'}$ chosen as the deepest wrong match in $T'$, which concludes the proof. ◁

It remains to show that the region endpoints in $T'$ are identical to those in $T'$. By condition **5**, the region of $T'$ containing $\lambda^r$ is necessarily included in $r$. If it is strictly included, then at least two bases paired in $T$ become unpaired in $T'$ so $T'$ contains at least one unpaired GU, which is excluded by Lemma 2. Then the region of $T'$ containing $\lambda^r$ is equal to $r$, and $T'$ matches the same regions as $T$, so $T' = T$. ◁

# 5 Algorithm for producing locked sequences

With the locked criterion in hand, we can now build an algorithm producing such sequences, and thus designs. The main result is summarized below.

**► Theorem 5.** In the model with either 2D or PK admissible structures, if $T$ is a structure featuring sufficiently large helices (as per Table 1), then it admits a design in the stacking energy model maxStacks computed in linear time and space by Algorithm 1.

**Proof.** Algorithm 1 builds a locked sequence in linear time and space by Lemmas 7 and 8, and such output is a solution by Theorem 3. ◁

Algorithm 1 produces locked sequences, provided input helices are large enough. The overall strategy is as follows: each key $\lambda^r$ contains a string composed of AA followed by a unique string (either globally in the PK case, or only within local loops in the 2D case). Using distinct keys results in the logarithmic lower bound on the helix size given in Table 1 ($k \ge \log_{3.56}(m) + 6.2$), within a constant factor of the theoretical lower bound provided by Lemma 1. Pattern AA (and its reverse complement UU) only occurs at such positions, so keys can only match keys, as required by conditions **4a** and **4b**. The design of region borders is controlled by the graph $\mathcal{N}$, defined below: each region is given a prefix corresponding to a
23:10 Designing RNA structures while maximizing stacking interactions

<figure><img src="image_17.png" /><figcaption>■ **Figure 6** Execution of Algorithm 1 on a sample target pseudoknotted structure (top). Helices are represented as wide arcs, unpaired positions as bullets. The resulting design is aligned with the target, with prefixes $P_r$ and suffixes $RevComp(P_{\bar{r}})$ in light font, and the rest (including key cores and there reverse complement) in bold font. The result of the key core and prefix picking procedure is given below, where the key cores are distinct and avoid AA and UU, and prefixes are picked in graph $\mathcal{N}$ (see below) according to region adjacencies.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>■ **Figure 7** Graph $\mathcal{N}$ of adjacent region extremities (cf Observation 6). Algorithm 1 assigns the nodes to all region prefixes, the suffix of matching regions being known by complementarity. Enter-nodes (resp. exit-nodes) are marked with incoming (outgoing) arrows. There exist paths from enter- to exit-nodes and cycles of any lengths (all avoiding CU except size-5 cycles).</figcaption></figure>

node in $\mathcal{N}$, and transitions between consecutive regions follow transitions in $\mathcal{N}$. This way, any interval extension beyond the bounds of the region would imply a certain pattern within the graph, and any such pattern has been carefully avoided in the construction of $\mathcal{N}$, as formalized in Observation 6 below.

► **Observation 6.** Let $\mathcal{N}$ be the graph defined in Figure 7. Nodes GC, CG are enter-nodes and GC, UC, CG are exit-nodes. Nodes $\{UC\}$, $\{AG\}$, $\{CG, CU\}$, $\{GC\}$ are respectively called U-, A-, C- and G-nodes. Then $\mathcal{N}$ satisfies the following properties (with the intuitive purpose for each property indicated above).

▷ Force a base in $\{G, U\}$ in each node and each node complement:

1. Each node has one base in $\{G, U\}$ and the other in $\{A, C\}$

▷ Avoid forbidden patterns AA and UU:

2. No node ends with A and only CU ends with U.

3. There is no arc between A- and U-nodes in either direction

▷ Avoid helix extensions with a pair $(x, \bar{x})$, (G, U) or (U, G) from adjacent regions:

4. No node has both in- and out-arcs to x-nodes for any $x \in \{A, U, C, G\}$.

5. No node has an in-arc from a C-node and an out-arc to a U-node.

6. No node has an in-arc from an A-node and an out-arc to a G-node.

▷ Avoid helix extensions with a pair (C, G) using an unpaired C:

7. No enter-node has an out-arc to a G-node

8. No exit-node has an in-arc from a C-node

▷ Assign nodes to any path or cycle, avoiding CU whenever possible:

9. For any $k \ge 1$, there exists a size-k path avoiding CU from an enter- to an exit-node.
■ **Algorithm 1** Design algorithm both for pseudoknotted or 2D sequences

**Input:** RNA structure *T* with "large enough" helices, sequence length *n*

**Output:** Design $\omega$ for *T* that is either PK-locked (PK case) or 2D-locked (2D case)

1: **function** DESIGN(*T*, *n*)

2: $\omega = X^n$

3: $\omega[i] = \mathbb{C}$ for all unpaired position $i$ in $T$

4: Pick nodes $P_r$ in $\mathcal{N}$ for each $r \in \mathcal{R}_T$ such that:

5: – if $\bar{r}_{end} + 1 = r'_{start}$ for some $r' \in \mathcal{R}_T$, then $\mathcal{N}$ has an arc $(P_r, P_{r'})$,

6: – if position $r_{start} - 1$ is unpaired or $< 1$, then $P_r$ is an enter-node

7: – if position $\bar{r}_{end} + 1$ is unpaired or $> n$, then $P_r$ is an exit-node

8: – $P_r \neq \text{CU}$ or $P_{\bar{r}} \neq \text{CU}$

9: Pick a key core $k_H \in \{\mathbb{A}, \mathbb{C}, \mathbb{U}, \mathbb{G}\}^*$ for each helix $H$ in $T$ such that:

10: – all strings $k_H$ have the same length

11: – no string $k_H$ starts with $\mathbb{A}$ or contains $\mathbb{AA}$ or $\mathbb{UU}$

12: – in the PK case: strings $k_H$ are pairwise distinct

13: – in the 2D case: $k_H \neq k_{H'}$ for any two helices $H, H'$ in the same loop.

14: **for** each helix $H = (r, \bar{r})$ of $T$ **do**

15: **if** $P_{\bar{r}} = \text{CU}$ **then**

16: $(r, \bar{r}) = (\bar{r}, r)$

17: $\omega[r] = P_r + \text{AAA}^* + k_H + \text{RevComp}(P_{\bar{r}})$;

18: ▷ AAA\* = fill remaining positions with As, minimum 2

19: $\omega[\bar{r}] = \text{RevComp}(\omega[r])$.

20: **return** $\omega$

Helix size $\ge$ 7 8 9 10 6+k $\lceil \log_{3.56}(m) + 6.2 \rceil$

Design structures with $m \le$ 3 9 35 125 $\frac{\sqrt{17}-1}{4} \left( \frac{\sqrt{17}+3}{2} \right)^k$ ($\approx 0.78 * (3.56)^k$)

■ **Table 1** Upper bound on parameter *m* (overall number of helices in the PK setting or maximum loop size in the 2D setting) as a function of the available helix size in the target structure, obtained with Lemma 11 (in appendix).

10. For any $k \ge 3$, there exists a size-k cycle containing at most one occurrence of CU.

We show the algorithm correctness in two steps: first in Lemma 7 we show that the algorithm does produce some output if helices are large enough. Then in Lemma 8 we show that the output is indeed a locked sequence (PK-locked in the PK case and 2D-locked in the 2D case).

► **Lemma 7.** If each helix in $\mathcal{T}$ has size at least $\lceil \log_{3.56}(m) + 6.2 \rceil$, with $m$ being the overall number of helices in $T$ (PK case) or the maximum loop degree (2D case), then Algorithm 1 does produce some output string $w$ in linear time.

► **Lemma 8.** Output $\omega$ is PK-locked (resp. 2D-locked) in the PK (resp. 2D) case.

**Proof.** In the following we use the notation $H = (r, \bar{r})$ where $r$ is chosen as in the algorithm, i.e. so that $P_{\bar{r}} \neq \text{CU}$. We write $\ell$ for the length of all key cores. With lines 17 and 19 of the algorithm, we have:

$$ \omega[r] = P_r A^i k_H \text{RevComp}(P_{\bar{r}}); \qquad \omega[\bar{r}] = P_{\bar{r}} \text{RevComp}(k_H) U^i \text{RevComp}(P_r). \quad (7) $$

where $i \ge 2$, $P_{\bar{r}} \neq \text{CU}$, strings $P_r$ and $k_r$ satisfy respectively the conditions lines 4-13.
23:12 Designing RNA structures while maximizing stacking interactions

361 ▷ Claim 9. There are no occurrences of AA or UU in $\omega$, other than those included in $A^i$ and
362 $U^i$ for each helix.

363 ▷ Claim 10. If regions $r$ and $r'$ satisfy $P_r = P_{r'}$, then $\{\omega[r_{\text{start}} - 1], \omega[\bar{r}'_{\text{end}} + 1]\}$ is not a
364 valid base pair (or some position is out of interval $[1, n]$)

365 We now prove that the successive conditions are met for the sequence $\omega$ to be locked.
366 Condition 1 is clear by Line 19. For Condition 2, Line 3 sets $C$ to all unpaired positions, and
367 lines 17 and 19 only affect paired positions.

368 Condition 3. Each region starts and ends with a node in $\mathcal{N}$ (or its reverse complement),
369 and by Observation 6.1 each node contains a base in $\{G, U\}$ and another in $\{A, C\}$ so its
370 reverse complement also contains a base in $\{G, U\}$. We thus define, for each region $r \in \mathcal{R}$,
371 $\lambda^r = [r_{\text{start}} + 1, r_{\text{end}} - 1]$. Then, $\omega[\lambda^r] = xA^iKy$ or $\omega[\lambda^r] = \bar{y}\text{RevComp}_kU^i\bar{x}$ for some base
372 $x \neq A$, some $i \ge 2$ and some key core $k$.

373 Condition 4a or 4b. Let $H = (r, \bar{r})$. Assume that $\lambda^r$ can be matched to an interval $I$ (the
374 proof is similar for $\lambda^{\bar{r}}$). Then $\omega[I] = \text{RevComp}(xA^i k_H y) = \bar{y}\text{RevComp}(k_H)U^i\bar{x}$ with $i \ge 2$
375 and $x \neq A$ so $\bar{x} \neq U$. By Claim 9, substring UU necessarily appears in $\omega$ within an interval $U^j$
376 for some key in helix $H' = (r', \bar{r}')$, and since $\bar{k}_H$ does not end with U (condition line 11) and
377 $\bar{x} \neq U$, we necessarily have $i = j$ and $I = \lambda^{\bar{r}'}$. Moreover, $\text{RevComp}(k_H') = \text{RevComp}(k_H)$,
378 so $k_H = k_{H'}$. In the PK setting, we necessarily have $H = H'$ by line 12, yielding condition 4a.
379 In the 2D setting, line 13 implies that $H$ and $H'$ are in distinct loops, thus satisfying
380 Condition 4b.

381 Condition 5. Let $(r, \bar{r})$ be an helix. Assume that an interval $I \supseteq \lambda^r$ aligns with an interval
382 $I'$, without using any $\{G, U\}$ pair within region $r$ (i.e. $I \cap r$ matches the corresponding
383 sub-interval of $I'$), we show that $I \subseteq r$. By Condition 4, $I'$ contains an interval $\lambda^{\bar{r}'}$ for some
384 helix $(r', \bar{r}')$, and both keys have the same size so $|r| = |r'|$. We first show that $I_{\text{end}} \le r_{\text{end}}$:
385 indeed, this is clear if $I_{\text{end}} < r_{\text{end}}$. Otherwise, we have the following relations:

386 $\omega[\bar{r}'_{\text{start}}, \bar{r}'_{\text{start}} + 1] = \text{RevComp}(\omega[r_{\text{end}} - 1, r_{\text{end}}])$ (match with $I$ within region $r$)

387 $\omega[\bar{r}'_{\text{start}}, \bar{r}'_{\text{start}} + 1] = \text{RevComp}(\omega[r'_{\text{end}} - 1, r'_{\text{end}}])$ (by construction of $\omega[r'], \omega[\bar{r}']$)

388 so $\omega[r_{\text{end}} - 1, r_{\text{end}}] = \omega[r'_{\text{end}} - 1, \bar{r}'_{\text{end}}]$ and $P_{\bar{r}} = P_{r'}$. Thus by Claim 10, $\{\omega[\bar{r}'_{\text{start}} - 1], \omega[r_{\text{end}} + 1]\}$ is not a valid base pair, and $I_{\text{end}} \le r_{\text{end}}$. Similarly if $I_{\text{start}} \le r_{\text{start}}$, then $P_r = P_{\bar{r}}$ and by
389 Claim 10, $\{\omega[r_{\text{start}} - 1], \omega[\bar{r}'_{\text{end}} + 1]\}$ is not a valid base pair, so $I_{\text{start}} \ge r_{\text{start}}$. $\triangleleft$

<figure><img src="image_19.png" /></figure>

6 Conclusion and discussion

391 In this work, we have introduced and investigated the inverse folding problem in the stacks
392 maximization (maxStacks) model. We have demonstrated that such a model lifts some
393 of the limitations induced by previous propositions, and in particular enables a design of
394 structures featuring complex loops. For input secondary structures with sufficiently large
395 helices, we proposed a linear-time algorithm that solves maxStacks inverse folding exactly.
396 Our algorithmic design hinges critically on a novel property, called locked, which represents
397 a sufficient condition for a sequence to be a solution to maxStacks inverse folding.

398 Despite unlocking the design in the maxStacks energy model, some desirable targets
399 still remain beyond the reach of our algorithms. Indeed, even if modestly growing with the
400 max degree (resp. number of helices), the minimal helix length requirement reduces the
401 scope of our algorithms, and our algorithms are no longer guaranteed to succeed whenever a
402 single helix is too short. However, we are confident that future extensions can work around
403 this issue. Plausible directions include an exhaustive enumeration of sequences for helices
404
of "small sizes", ensuring an absence of alternative matchings for regions, possibly leading to a Fixed-Parameter Tractable algorithm in the overall length of small helices. Another extension would be to merge helices involved in loops of degree 2, associated with bulges and internal loops.

An interesting extension would be to consider weighted models, with each stack having a different real-valued contribution, depending of its base contents. This would be an important step towards implementing a design algorithm in the full Turner model, where general loops contribute to the structure energy. In such a setting, Wobble pairs ({G, U}) are typically less stable than others, see Table 2 in appendix. Remind that our algorithm does not use any {G, U} and induces more stacks in the target than all other structures (even including {G, U}). It follows that our algorithm already solves the design problem in any weighted stacking models where stacks without {G, U} have greater energetic contribution than those including {G, U}.

In terms of theoretical complexity, this work leaves open the possible NP-hardness of inverse folding in the maxStacks model. The closest related hardness result is from Bonnet et al. [2], who showed NP-hardness for the maxBPs model with pre-defined bases; it remains open whether inverse folding is hard when we lift this restriction under either model. Moreover, the existence of locked sequence for general targets represents an interesting problem, whose complexity remains open. Finally, one could consider the minimal perturbation, e.g. with respect to base pair insertions, such that the modified structure admits a locked sequence. This work already implies that inserting $O(H \log(\Delta))$ (resp. $O(H \log(H))$) base pairs is always sufficient in the 2D (resp. PK) setting, with $H$ the number of helices and $\Delta$ the max degree of a loop. Moreover, this bound is tight up to a multiplicative constant in the worst case, as a direct corollary of Lemma 1. This problem is of specific interest towards applications, as it is known that minor structural variations (e.g. helix extensions) are unlikely to alter biological function, and are often allowed in design methodologies [14].

References

1. Turner 2004 terminal mismatch parameters. URL: https://rna.urmc.rochester.edu/NNDB/turner04/tm.html.

2. Édouard Bonnet, Paweł Rzążewski, and Florian Sikora. Designing RNA secondary structures is hard. *Journal of Computational Biology*, 27(3):302–316, March 2020. doi:10.1089/cmb.2019.0420.

3. Théo Boury, Samuel Gardelle, Laurent Bulteau, and Yann Ponty. RNA inverse folding can be solved in linear time for structures without isolated stacks or base pairs. 20(1):20. doi:10.1186/s13015-025-00278-6.

4. Théo Boury, Leonhard Sidl, Ivo L. Hofacker, Yann Ponty, and Hua-Ting Yao. Old dog, new tricks: Exact seeding strategy improves RNA design performances. In Sriram Sankararaman, editor, *Research in Computational Molecular Biology*, pages 134–152. Springer Nature Switzerland. doi:10.1007/978-3-031-90252-9_9.

5. Alexander Churkin, Matan Drory Retwitzer, Vladimir Reinharz, Yann Ponty, Jérôme Waldispühl, and Danny Barash. Design of RNAs: comparing programs for inverse RNA folding. *Briefings in Bioinformatics*, page bbw120, January 2017. doi:10.1093/bib/bbw120.

6. Raymond F. Gesteland, Thomas R. Cech, and John F. Atkins, editors. *The RNA world*. Number 43 in Cold Spring Harbor monograph series. Cold Spring Harbor Laboratory Press, Cold Spring Harbor, NY, 3. ed. edition, 2006.

7. Jozef Haleš, Alice Héliou, Ján Maňuch, Yann Ponty, and Ladislav Stacho. Combinatorial RNA design: Designability and structure-approximating algorithm in watson–crick and nussinov–jacobson energy models. 79(3):835–856. doi:10.1007/s00453-016-0196-x.
23:14 Designing RNA structures while maximizing stacking interactions

453

454

8 Michelle L. Hastings and Adrian R. Krainer. Rna therapeutics. RNA, 29(4):393-395, March 2023. doi:10.1261/rna.079626.123.

455

456

457

9 I. L. Hofacker, W. Fontana, P. F. Stadler, L. S. Bonhoeffer, M. Tacker, and P. Schuster. Fast folding and comparison of RNA secondary structures. 125(2):167-188. Number: 2. doi:10.1007/BF00818163.

458

459

10 Rune B. Lyngsø. Complexity of Pseudoknot Prediction in Simple Models, pages 919-931. Springer Berlin Heidelberg, 2004. doi:10.1007/978-3-540-27836-8_77.

460

461

462

463

464

11 David H. Mathews, Matthew D. Disney, Jessica L. Childs, Susan J. Schroeder, Michael Zuker, and Douglas H. Turner. Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure. 101(19):7287-7292. Publisher: Proceedings of the National Academy of Sciences. URL: https://www.pnas.org/doi/abs/10.1073/pnas.0401799101, doi:10.1073/pnas.0401799101.

465

466

467

12 R Nussinov and A B Jacobson. Fast algorithm for predicting the secondary structure of single-stranded RNA. Proceedings of the National Academy of Sciences, 77(11):6309-6313, November 1980. doi:10.1073/pnas.77.11.6309.

468

469

470

13 Erik Poppleton, Niklas Urbanek, Taniya Chakraborty, Alessandra Griffo, Luca Monari, and Kerstin Göpfrich. Rna origami: design, simulation and application. RNA Biology, 20(1):510-524, July 2023. doi:10.1080/15476286.2023.2237719.

471

472

473

474

475

14 Matan Drory Retwitzer, Vladimir Reinharz, Alexander Churkin, Yann Ponty, Jérôme Waldispühl, and Danny Barash. incarnafbinv 2.0: a webserver and software with motif control for fragment-based design of rnas. Bioinformatics, 36(9):2920-2922, 01 2020. arXiv:https://academic.oup.com/bioinformatics/article-pdf/36/9/2920/48986446/bioinformatics_36_9_2920.pdf, doi:10.1093/bioinformatics/btaa039.

476

477

15 Siyu Wang, Drew Weissman, and Yizhou Dong. Rna chemistry and therapeutics. Nature Reviews Drug Discovery, 24(11):828-851, July 2025. doi:10.1038/s41573-025-01237-x.

478

479

480

481

16 Hua-Ting Yao, Cédric Chauve, Mireille Regnier, and Yann Ponty. Exponentially few RNA structures are designable. In Proceedings of the 10th ACM International Conference on Bioinformatics, Computational Biology and Health Informatics, BCB '19, pages 289-298. ACM, September 2019. doi:10.1145/3307339.3342163.

482

483

484

17 Jeffrey Zuber, Susan J Schroeder, Hongying Sun, Douglas H Turner, and David H Mathews. Nearest neighbor rules for RNA helix folding thermodynamics: improved end effects. Nucleic Acids Research, 50(9):5251-5262, May 2022. doi:10.1093/nar/gkac261.
A Energy contributions of different base pairs

<table><tr><td>Type of stacks</td><td>{ {G, C}, {G, C}}</td><td>{ {A, U}, {G, C}}</td><td>{ {A, U}, {A, U}}</td></tr><tr><td>Energy (kcal/mol)</td><td>-3,075</td><td>-2,195</td><td>-1,0725</td></tr><tr><td>Type of stack</td><td>{ {G, C}, {G, U}}</td><td>{ {A, U}, {G, U}}</td><td>{ {G, U}, {G, U}}</td></tr><tr><td>Energy (kcal/mol)</td><td>-1,89</td><td>-1,045</td><td>+1</td></tr></table>

Table 2 Average energy of stacks using any two base pairs in the Turner model [11], disregarding uncertainty of measures (values taken from [1]). Averages are taken over all possible configurations of the two base pairs. Note that replacing any {A, U} or {G, C} base pair by {G, U} increases the energy (thus reducing the stability). It overall suggests that the local stability of the stack, and thus the global stability of the structure, is degraded by the presence of {G, U}.

B Proofs for Section 3 (Why study inverse folding in the maxStacks model?)

**Proof of Lemma 1.** Note that $h < \log_6(\Delta) \leftrightarrow 6^h < \Delta$, which means that $\mathcal{L}$ contains more than $6^h$ helices. By contradiction, suppose that there exists a design $\omega$ of $\mathcal{L}$. As there are no more than $6^h$ assignment possibilities for the helices of size $h$, there exists necessarily two helices $H_1 = (r_1, \bar{r}_1)$ and $H_2 = (r_2, \bar{r}_2)$ such that $r_1 = r_2$ and $\bar{r}_1 = \bar{r}_2$. In that case, one could construct $\mathcal{L}'$ containing the same helices as $\mathcal{L}$, with the exception of two helices $H'_1 = (r_2, \bar{r}_1)$ and $H'_2 = (r_1, \bar{r}_2)$ replacing $H_1$ and $H_2$ respectively. Note that $H'_1$ and $H'_2$ are indeed feasible since $r_1 = r_2$ and $\bar{r}_1 = \bar{r}_2$. In addition, $H'_1$ and $H'_2$ feature as much stacks as $H_1$ and $H_2$, making $\mathcal{L}'$ a competitive alternative in maxStacks. ◀

<figure><img src="image_20.png" /></figure>

C Proofs for Section 5 (Algorithm for producing locked sequences)

**► Lemma 11.** A valid key core is a string over {$A, U, C, G$} starting with {$U, G, C$} and avoiding AA and UU. For any $k$ and $m \le \frac{\sqrt{17}-1}{4} \left(\frac{\sqrt{17}+3}{2}\right)^k$ (or $k=1$ and $m \le 3$), there exist at least $m$ distinct valid key cores of length $k$. Conversely, for any $m$, there exist $m$ distinct valid key cores of length $k$ for any $k \ge \log_{3.56}(m) + 0.2$.

**Proof.** In the following we write $r$ and $c$ for the positive roots of the equations $c = 2r + 2$ and $r + 2 = rc$, i.e. $c = \frac{\sqrt{17}+3}{2} \approx 3.56$ and $r = \frac{\sqrt{17}-1}{4} \approx 0.78$. We write $m_x(k)$ for the number of length-$k$ strings starting with $x \in \{A, U, G, C\}$ and avoiding AA and UU. We show by induction that:

$$m_A(k) = m_U(k) \ge rc^{k-1} \quad \text{and} \quad m_C(k) = m_G(k) \ge c^{k-1}$$

Indeed, $m_A(k) = m_U(k)$ and $m_G(k) = m_C(k)$ are clear by symmetry. For $k=1$, $m_A(k) = 1 \ge rc^0$ and $m_G(k) = 1 = c^0$. Also, we have the following recursive relations :

$$ \begin{aligned} m_A(k+1) &= m_U(k) + m_G(k) + m_C(k) & \ge (r+2)c^{k-1} &= rc^k \\ \text{and } m_G(k+1) &= m_A(k) + m_U(k) + m_G(k) + m_C(k) & \ge (2r+2)c^{k-1} &= c^k \end{aligned} $$

Thus, the number of valid key cores is at least $m_U(k) + m_C(k) + m_G(k) = m_A(k+1) = rc^k$. Finally, the logarithmic bound is obtained as follows: for $k \ge \log_{3.56}(m) + 0.2$, using $c \ge 3.56$ and $\log_c(r) \ge -0.2$ we have $k \ge \log_c(m) - \log_c(r)$ so $rc^k \ge m$: there are indeed $m$ distinct valid keys of size $k$. ◀

<figure><img src="image_21.png" /></figure>
**Proof of Lemma 7.** We show that the algorithm can always find nodes $P_r$ satisfying the criteria lines 5 to 8, and key cores satisfying the criteria lines 10 to 13.

Regarding nodes, we consider the graph $\mathcal{G}$ of adjacent successor/predecessors: use one node per region of $\mathcal{R}_T$, and add an arc between $r$ and $r'$ if $\bar{r}_{\text{end}} + 1 = r'_{\text{start}}$ (equivalently, $r'$ is the successor of $r$ and $\bar{r}, r'$ are adjacent). Note that each region has in-degree at most one and out-degree at most one, so it partitions into paths and cycles (subsets of loops). Moreover, if $r$ has out-degree 0 then $\bar{r}_{\text{end}} + 1$ is unpaired (or $> n$), and if $r$ has in-degree 0 then $r_{\text{start}} - 1$ is unpaired (or $< 1$). Thus, lines 5 to 7 are satisfied if each connected component of $\mathcal{G}$ can be mapped to a (not necessarily simple) path or cycle of the length in $\mathcal{N}$, where nodes with in-degree 0 (resp. out-degree 0) are mapped to enter-nodes (resp. exit-nodes) of $\mathcal{N}$. This is always possible by Observations 6.9 and 6.10. Finally, satisfying line 8 is somewhat more technical. We show that nodes can be assigned to regions sequentially such that, at each step, there is at most one region assigned CU whose complement does not have any assigned node. Such a region $r$ with $P_r = \text{CU}$ and $P_{\bar{r}}$ still undefined is hereafter called a *single CU-region*. This is achieved by picking the connected component containing the complement of a single CU-region (if any, otherwise pick any component). If the component is a path, by Observation 6.9 we assign non-CU nodes to all regions: there is no single CU-region after such step. If the component is a cycle, we assign nodes so that at most one node $r$ of the cycle is assigned CU with Observation 6.10 ($r$ is any other region than the complement of the previous CU-region); there is at most one single CU-region at this point.

Regarding key cores, let $\ell = \lceil \log_{3.56}(m) + 0.2 \rceil$. By Lemma 11, key cores can be built with length $\le \ell$, and assigned successively to helices. In the 2D case, note that condition line 13 can easily be satisfied by assigning key cores in a “top-down” ordering of loops: first from the outer loop, then in the intervals inside each helix recursively, so that at most one helix in each loop has an already assigned key core and all other helices can be assigned the $m-1$ available strings. Overall each key core has length $\le \ell$, and lines 17 and 19 design helices with strings of size at least $\ell + 6$, which corresponds to the required bound on the helix size.

The overall algorithm indeed takes linear time: paths and cycles in $\mathcal{G}$ can be enumerated in linear time, and assigned pre-selected paths in $\mathcal{N}$ based on their lengths (the constraint over nodes CU only forces the components to be visited in the right order, without affecting the overall complexity). Then key cores can also be computed in linear time for the total number of helices. Finally, each helix is also designed in linear time of its length. ◀

**Proof of Claim 9.** Such patterns cannot be included in strings $P_r$, $P_{\bar{r}}$ or there reverse complements (nodes or $\mathcal{N}$ are neither AA nor UU), nor in $k_r$, $\bar{k}_r$ (by line 11). They cannot overlap two such parts: $P_r$ does not end with A (Observation 6.2), $k_r$ does not start with A, and $P_{\bar{r}}$ does not end with U and RevComp($P_{\bar{r}}$) does not start with A nor U (using Observation 6.2 and $P_{\bar{r}} \ne \text{CU}$). So patterns AA and UU cannot be included in a single region (beside $A^i$ and $U^i$ in the reverse complement). They cannot use unpaired positions either (since such positions receive a C, as per line 3), and they cannot either use two endpoints of consecutive regions. Indeed, if $\omega[r_{\text{end}}] = \omega[r'_{\text{start}}] = A$ (the proof is similar for U) for any two regions $r, r'$, then $P_{r'}$ is an A-node (since $P_{r'}$ starts with A) and $P_{\bar{r}}$ is a U-node (since RevComp($P_{\bar{r}}$) ends with A). By Observation 6.3 there is no arc from $P_{\bar{r}}$ to $P_{r'}$ and by line 5, $r_{\text{end}} \ne r'_{\text{start}} - 1$, so these two regions do not form any AA pattern. ◁

**Proof of Claim 10.** See Figure 8. Write $P = P_r = P_{r'}$, $x = \omega[r_{\text{start}} - 1]$ and $y = \omega[\bar{r}'_{\text{end}} + 1]$. If $r_{\text{start}} - 1$ is an unpaired position, then $x = C$ and $P_r$ is an enter-node (by line 6). Otherwise, $r_{\text{start}} - 1$ is the last position of a region $q$ with $P_{\bar{q}} \to P$ (line 5), and $X = P_{\bar{q}}$ is
<figure><img src="image_22.png" /><figcaption>■ **Figure 8** Illustration of Claim 10, providing the core of the non-extension condition: if the prefix of a region matches the suffix of another (including with both regions of an helix, but not necessarily), then the match cannot be extended further beyond the region border, even with `{G, U}` pairs. Here we write $\bar{P}$ for $\text{RevComp}(P)$, and $\bullet$ is used as a wildcard for any base, i.e. $y\bullet$ is a $y$-node of $\mathcal{N}$.</figcaption></figure>

an $\bar{x}$-node in $\mathcal{N}$. Similarly, $\bar{r}'_{\text{end}} + 1$ is either unpaired so $y = C$ and $P = P_{r'}$ is an exit-node (line 7), or it is paired to a region $q$ where $Y = P_q$ is a $y$-node and $P \to Y$ in $\mathcal{N}$.

556

557

Overall, if both are unpaired then $x = y = C$; if only $r_{\text{end}} + 1$ is unpaired then $P$ is an exit-node and by Observation 6.8, $y = C$ and $\bar{x} \neq C$ so $\{x, y\}$ is not valid. If only $\bar{r}'_{\text{start}} - 1$ is unpaired then $P$ is an enter-node and by Observation 6.7, $x = C$ and $y \neq G$ so $\{x, y\}$ is not valid. Otherwise we have $X \to P \to Y$ in $\mathcal{N}$ with $X$ an $\bar{x}$-node and $Y$ a $y$-node. By observation 6.4, $\bar{x} \neq y$ so $\{x, y\}$ is neither {A, U} nor {G, C}. By Observation 6.5 $(x, y) \neq (G, U)$ and by Observation 6.6, $(x, y) \neq (U, G)$. Overall $\{x, y\}$ is not a valid base pair.

558

559

560

561

562

563

564

565

∠