SamplingDesign: RNA design via continuous optimization with coupled variables and Monte-Carlo sampling

Domain: chemistry
Published: 2026-02-20
Editorial outcome: minor_revision

Abstract
RNA design aims to find a sequence that can fold into a target secondary structure. It can create artificial RNA molecules for specific functions, with wide applications in medicine. It is computationally challenging due to two levels of combinatorial explosion: the exponentially large design space and the exponentially many competing structures per design. Popular methods such as local search cannot keep up with these combinatorial explosions. We instead employ two techniques from machine learning, continuous optimization and Monte-Carlo sampling. We start from a distribution over all valid sequences, and use gradient descent to improve the expectation of an arbitrary objective function. We define novel coupled-variable distributions to model the correlation between nucleotides. We then use sampling to approximate the objective, estimate the gradient, and select the final candidate. Our work consistently outperforms state-of-the-art methods in key metrics including Boltzmann probability and ensemble defect, especially on long and hard-to-design structures.

Introduction
Ribonucleic acid (RNA) is vital for fundamental cellular processes such as transcription and translation, catalyzing reactions, and controlling gene expression1–3. Its importance is also evidenced by COVID-19, an RNA virus, as well as the two recent Nobel Prizes in Physiology and Medicine for messenger RNA vaccines (2023) and microRNAs (2024).

The problem of RNA design, also known as “inverse folding”, aims to find sequences that are capable of folding into a target secondary structure4–12. This process enables the creation of artificial RNA molecules with specific function, such as artificial ribozymes13,14, artificial miRNAs15, artificial RNA aptamers16, and artificial riboswitches17,18.

Computationally, the RNA design problem is extremely challenging due to two levels of combinatorial explosion. First, the total sequence space contains 4n candidates, i.e., {A, C, G, U}n (n is the sequence length); and secondly, for each sequence, there are exponentially many alternative structures competing with the target structure (Fig. 1). A good design should not only fold into the target structure as its minimum free energy (MFE) structure, but also have a sharp energy landscape (Fig. 1b) which results in high equilibrium probability of the target structure (Fig. 1c); a bad design, by contrast, has a flatter energy landscape with competing structures at similar energies, which results in the target structure having a low probability in the ensemble. This is the same principle guiding protein design19.

Due to the above difficulties, RNA design has been proved NP-hard for a simplified energy model20. Therefore, the most popular approach for this problem uses heuristic methods such as local search4–9, which starts with a sequence and tries to optimize an objective function by revising one or a few nucleotides in each step. However, such methods are only capable of exploring a fixed number of candidate sequences, thus not able to keep up with the exponential growth of the design space. As a result, they tend to perform poorly on longer and harder-to-design structures.

Instead, we borrow two techniques from machine learning, continuous optimization21,22 and Monte-Carlo sampling, for the RNA design problem. The basic idea is to start with a distribution over all possible candidate sequences and use gradient descent to gradually sharpen the distribution, with each step changing all positions simultaneously in contrast to local search methods. We first define novel sequence distributions for any target structure using coupled variables for paired and mismatch positions, which not only rules out invalid sequences but also models the positional correlations explicitly. We then aim to optimize the expectation of an arbitrary objective function over these distributions. However, given the diversity of various objective functions in RNA design (such as Boltzmann probability or ensemble defect), it is computationally prohibitive to compute the exact values of the expected objective function or its gradient over the whole distribution of exponentially many sequences (and the exponentially many structures per sequence). Therefore, in contrast to previous work, we use sampling to approximate the expected objective function and its gradient. At the end, we return the best sample in terms of the objective function among all samples collected, which yields high-quality designs.

When tested on the Eterna100 benchmark, our work consistently outperforms the state-of-the-art RNA design methods4,5 in key metrics such as Boltzmann probability and ensemble defect. The advantage of our work is especially salient on longer and harder-to-design structures, demonstrating the benefit of this distributional approach modeling the design space over local search methods operating on a single sequence.

Results
RNA design problem and valid design space
In RNA design, we are given a target (pseudoknot-free) secondary structure y and want to find an RNA sequence x that can fold into it. We represent y as a well-balanced dot-bracket string y=y1y2...yn with each yi ∈{.,(,)}, and x as a nucleotide string x=x1x2...xn with each xi ∈{A,C,G,U}. For each index i, if yi =., its corresponding nucleotide xi must be unpaired in the folded structure; if yi =(, xi must be paired with some downstream xj; and if yi = ), xi must be paired with some upstream xj. The set of unpaired indices in y is denoted unpaired(y) and the set of paired indices pairs(y). For example, if y=(((.))), we have unpaired(y)={4,5,6} and pairs(y)={(1,9), (2,8), (3,7)}.

It is a common misconception to think that the design space has all 4n sequences; in fact, most of that space is invalid for any target structure due to dependencies in base-paired positions (each base-pair is a Watson-Crick-Franklin or wobble pair). For example in Fig. 2(a), if the first nucleotide is C, then the last one must be G. More formally, for each (i,j) ∈ pairs(y), the choices of xi and xj are no longer independent, where xixj must be in the set of valid pairs P={CG, GC, AU, UA, GU, UG}. Therefore, we refer to those 4n sequences as the “total sequence space” (Fig. 1a) and formulate the “valid design space” X(y) as

X(y) = {x ∈ {A,C,G,U}n | ∀(i,j) ∈ pairs(y), xixj ∈ P} (1)

Although this concept has appeared in previous studies23,24, this is the first time it has been rigorously defined. By ruling out invalid sequences, the valid design space is much smaller than the total sequence space (Fig. 2(b)):

|X(y)| = 4|unpaired(y)| × 6|pairs(y)| ≪ 4n (2)

Fig. 2(c) shows that the ratio of valid ratio in the total sequence space decreases exponentially, being <1/1060 for some Eterna100 structure.

Overview of the SamplingDesign algorithm
Conventionally, RNA design is formulated as a discrete optimization problem of finding an RNA sequence that optimizes a given objective f(x,y) (e.g., Boltzmann probability p(y|x) or ensemble defect NED(x,y), see Methods). We instead formulate RNA design as a continuous optimization problem which starts with a distribution py(·) over candidate sequences and then uses gradient descent to gradually sharpen the distribution in order to optimize the expectation Ex∼py(·) [f(x,y)].

Independent vs. Coupled-Variable Distributions. The simplest way of modeling the sequence distribution py(·) is to assume independence among all positions21. In this model, each position has a distribution over four nucleotides {A,C,G,U}, and the sequence distribution is the product of these independent distributions (Fig. 2(d)). However, this is a distribution over the total sequence space of 4n candidates, which, as discussed above, contains mostly invalid sequences for the target structure. So we instead define a novel distribution for the valid design space using coupled variables (Fig. 2(e)), where each base-pair (such as (1,9) in Fig. 2(a)) has a single distribution over six valid base pairs {GC, CG, AU, UA, GU, UG} (rather than 4×4=16 choices).

More importantly, in the independent-variable model, for a base pair (i,j), the gradients for positions i and j are over 4×2=8 variables (two independent sets). By contrast, in the coupled-variable model, the gradients for i and j are only over the six valid base pairs. As a result, we suspect our coupled-variable model is less likely to get stuck in a local minima.

Moreover, in the standard RNA energy model (Turner rules25), there are many other correlated positions, which we classify as mismatches (such as (4,6) in Fig. 2(a)) and trimismatches (see Methods). Each of these groups of positions is also coupled as a single distribution to model the correlations explicitly.

Sampling for objective evaluation, gradient estimation, and candidate exploration. SamplingDesign then aims to minimize the expectation of the given objective with respect to the distribution, i.e., Ex∼py(·) [f(x,y)], via gradient descent (Fig. 3). However, due to the complexity of objective functions in RNA design, computing the exact value of its expectation is generally intractable. Instead, we use Monte-Carlo sampling to approximate the objective function and estimate its gradient. We update the distribution using the estimated gradient, shifting the probability mass toward more promising candidates (Fig. 3d,e). We repeat this process until the objective function converges (or after reaching the maximum number of steps). We then return the best sample (in terms of the objective function f(x,y)) encountered during the optimization as the final design; we found this sampling-based candidate exploration method to result in better solutions than the commonly used alternative of picking the max-probability sequence argmax x py(x) within the distribution (Fig. 3f).

Evaluation setup
We compared SamplingDesign with other state-of-the-art methods (SAMFEO4 and NEMO5) and another continuous optimization approach (Matthies et al.21) on four secondary structure datasets (Eterna100, Eterna100-v2, Rfam, and RNAsolo-764; see Methods for details) using the following metrics:

1. Average Boltzmann probability across all puzzles (arithmetic mean). The Boltzmann probability p(y|x) is the probability of the target structure y given the design x in the ensemble of all possible structures (see Eq. (6)).
2. Geometric mean of Boltzmann probability across all puzzles (putting more weight on harder-to-design puzzles with lower probability) except those 18 that are proven to be undesignable by our previous work26,27. We remove these undesignable puzzles because their extremely low probabilities bias the geometric mean towards 0.
3. Average normalized ensemble defects across all puzzles. Normalized ensemble defect12 NED(x,y) is average fraction of incorrectly paired nucleotides at equilibrium (over the ensemble of all structures) (see Eq. (11)).
4. Average structural distance between the target structure y and the MFE structure of the design x across all puzzles. The structural distance d(y,y0) between two structures is the number of nucleotides that fold differently (see Eq. (9)).
5. Average free energy gap across all puzzles. The free energy gap ΔΔG(x,y) is the difference in free energy between the target structure y and the MFE structure of the design x (see Eq. (10)).
6. Number of puzzles where an MFE solution is found. An MFE solution is an RNA sequence for which the target structure is one of its MFE structures (see Eq. (7)).
7. Number of puzzles where a unique MFE (uMFE) solution is found. A uMFE solution is a design x where the target structure y is x’s only MFE structure (see Eq. (8)). We use the Turner 2004 parameters as implemented in ViennaRNAfold v2.x; this set of parameters is the most commonly used energy model in recent RNA folding and RNA design work.

SamplingDesign outperforms state-of-the-art methods
See Fig. 4(a) for results of SamplingDesign and other baselines (SAMFEO4 and NEMO5) on Eterna100. On Eterna100, SamplingDesign ranks first on all metrics except average structural distance and the number of uMFE-solution puzzles, where it ranks second. Notably, optimizing Boltzmann probability yields the best results across all other metrics, even better than optimizing those metrics directly. For example, when optimizing for Boltzmann probability, SamplingDesign achieves the lowest average NED of 0.035, outperforming both itself and SAMFEO optimizing for NED (0.041 and 0.037, resp.). The fact that optimizing Boltzmann probability also optimizes other metrics is consistent with the observations of Ward et al.28.

See also Tables S1–S3 for results on Eterna100-v2, Rfam, and RNAsolo. Across all datasets, SamplingDesign consistently performs the best in terms of the more important ensemble-based metrics: Boltzmann probability and NED. Against SAMFEO, SamplingDesign achieves significantly higher Boltzmann probability on Eterna100 (p=0.009) and Eterna100-v2 (p=0.0007) by a one-sided Wilcoxon signed rank test; the differences against NEMO are even more significant on all datasets. SamplingDesign also outperforms all baselines by a wide margin on the geometric mean of Boltzmann probability, which indicates that SamplingDesign is much more effective for long and hard-to-design puzzles (see also Fig. 4b, c). Fig. 4(d) shows a Venn diagram of puzzles solved under the uMFE criterion by our method, SAMFEO, and NEMO. Of the 100 puzzles, 79 are designable and 18 are proven to be undesignable26,27. The remaining 3 puzzles (#68, #97, #100) have yet to be solved or proven undesignable. Notably, our method solved puzzle #89 which was not solved by any RNA design paper, although our previous work26 serendipitously solved it when trying to prove undesignability.

Comparison with NEMO. SamplingDesign surpasses NEMO by one for the number of MFE solutions (the highest in literature) and trails by one for the number of uMFE solutions. However, some of NEMO’s advantage relies on heuristic rules hand-designed for the ViennaRNA energy model. Previous work4 showed that NEMO without those rules solves only 76 and 75 puzzles under the MFE and uMFE criteria, respectively. Thus, applying these heuristics increases the number of solved puzzles by 3 (MFE) and 2 (uMFE). More importantly, NEMO’s solution quality is often much worse than ours in terms of the more meaningful ensemble-based criteria. For example, while NEMO solved puzzles #71 and #79, their solutions have far worse Boltzmann probability and NED due to competitions from alternative structures (Fig. S1).

Comparison with Matthies et al. 2023. Lastly, SamplingDesign is compared with a different continuous optimization approach21 (see Methods for detailed comparisons). In their work, the authors reported the Boltzmann probability of the 18 shortest puzzles (up to 50 nucleotides) using an 80GB NVIDIA A100 GPU. We were able to rerun their code up to the shortest 51 puzzles (up to 104 nucleotides) using 80GB NVIDIA H100 GPU under their default settings. Their work is not able to scale to longer puzzles due to the GPU memory limit. By contrast, SamplingDesign is able to scale to the longest structures in the Eterna100 benchmark (400 nt) and produces designs that outperform theirs in all metrics and by very large margins (Table S5 and Fig. S2).

Substantially better designs for long and hard-to-design puzzles
Figures 4b, c and S3a–c compare SamplingDesign to SAMFEO by puzzle length groups. Although both methods perform similarly on short puzzles, ours excels on the four longest groups (116–192, 200–316, 337–387, and 389–400). Our geometric mean consistently outperforms SAMFEO (Fig. 4b) and lowers NED (Fig. 4c). In addition, Figs. S3a, b show substantial improvements for the free energy gap and structural distance on the three longest groups, and Fig. S3c shows noticeable improvements for the arithmetic mean of Boltzmann probability in the four longest groups.

Figs. S3d–f compare SamplingDesign and SAMFEO on individual puzzles. On the original scale of Boltzmann probability (Fig. S3d), both methods achieve similar performance on most puzzles, but ours substantially outperforms SAMFEO on some of them, increasing the probability from 0.325 to 0.485 for puzzle #74 and 0.550 to 0.701 for puzzle #38. On a log scale (Fig. S3e), however, these differences become more salient for puzzles longer than 280 nucleotides, which inherently have lower Boltzmann probability. Our solutions exceed SAMFEO’s by more than a factor of 10 for nine puzzles (annotated in the figure), and corresponding improvements appear for the same puzzles in NED (Fig. S3f).

Figs. 5–6 and S4–S6 provide detailed comparisons between SamplingDesign and SAMFEO for puzzles #73, #78, #76, #91, and #99, respectively. In all cases, our method substantially outperformed SAMFEO across all metrics, and our solutions fold to structures that are identical or very close to the target structures, while SAMFEO solutions fold to vastly different structures. See Table S6 for a summary.

For puzzle #73 (Fig. 5), SamplingDesign achieved a substantially higher Boltzmann probability than SAMFEO (0.005 vs. 3×10−27). More interestingly, our solution folds into a unique MFE structure that is identical to the target structure, thus satisfying both the MFE and uMFE criteria (Fig. 5b). By contrast, the SAMFEO solution folds into a vastly different structure with a structural distance of 138 (Fig. 5c), i.e., out of all the 370 nucleotides, ~37.3% do not fold according to y, which is consistent with its NED of 41.7% (Fig. 5a). We also visualized the base-pairing probabilities of the two solutions (Fig. 5d, e), where ours shows highly probable correct pairs (in blue) while SAMFEO’s shows many strong incorrect pairs (in red) and some correct pairs with probabilities <0.1 (in orange). This suggests that our solution’s Boltzmann ensemble is sharper, while the SAMFEO solution suffers from competing structures, which make the target structure highly unfavorable in its ensemble (see also Fig. 1). Puzzle #76 (Fig. S4) is similar, where our solution folds into a unique MFE structure that is identical to the target structure, with a Boltzmann probability of 0.046, while the SAMFEO solution has numerous alternative base-pairings (Fig. S4e), resulting in a very low probability of 7×10−8.

Figure 6 visualizes the undesignable puzzle #78, with its undesignable motif highlighted (a motif is a set of consecutive loops). Our solution’s MFE structure is almost identical to the target structure, except for the missing two pairs in that motif. Therefore, we believe this structural distance of 4 is the smallest achievable distance. By contrast, SAMFEO’s MFE structure differs greatly from the target structure, with most correct pairs having low probabilities (<0.1, orange arcs) and many incorrect pairs having high probabilities (red arcs). This competition from alternative structures makes the Boltzmann probability of SAMFEO’s solution extremely low. See Fig. S5 for a similar undesignable puzzle #91, where the MFE structure of our design only missed the four pairs inside two undesignable motifs. Puzzle #99 (Fig. S6) is also undesignable, and our MFE structure has a much smaller structural distance than SAMFEO’s (20 vs. 80).

Methods
RNA folding
The ensemble of an RNA sequence x is the set of all possible secondary structures of x, denoted as Y(x). In thermodynamic RNA folding models, Gibbs free energy change ΔG(x,y) is used to characterize the stability of a structure y ∈ Y(x). The structure with the minimum free energy is the most stable structure in the ensemble, i.e., the MFE structure,

MFE(x) = argmin y∈Y(x) ΔG(x,y) (3)

Note that most secondary structure prediction tools use arbitrary tie-breaking when there are multiple lowest free energy structures. This issue was often neglected in the literature, but it deserves clarification here. To be precise, we define

MFEs(x) = {y | ΔG(x,y) = min y0∈Y(x) ΔG(x,y0)} (4)

and say x has a unique MFE (uMFE) structure.

The partition function sums the contribution of all structures in the ensemble:

Q(x) = ∑ y∈Y(x) e−ΔG(x,y)/RT (5)

where R is the molar gas constant and T is the absolute temperature. Accordingly, the Boltzmann probability of a sequence x folding into a structure y is defined as

p(y|x) = e−ΔG(x,y)/RT / Q(x) (6)

RNA design objectives: MFE vs. ensemble
Given a target structure y, RNA design aims to find a suitable RNA sequence x that can naturally and easily fold into y. There are different ways to quantify how “easily” x folds into y, which we categorize into two broad categories: (a) MFE-based and (b) ensemble-based criteria; see Fig. 1.

MFE criterion and uMFE criterion. A sequence x is said to be an MFE solution of y if y is one of the MFE structures of x:

y ∈ MFEs(x) (7)

Or equivalently, ∀y0≠y, ΔG(x,y) ≤ ΔG(x,y0): As a stricter criterion, a sequence x is said to be a uMFE solution of y if y is the unique MFE structure of x, or equivalently:

∀y0≠y, ΔG(x,y) < ΔG(x,y0) (8)

To search for an MFE or uMFE solution, we can gradually update the sequence to minimize one of the following metrics:

• structural distance d(y,MFE(x)), where d(·,·) is a standard distance metric between two secondary structures, returning the number of differently folded nucleotides:

d(y,y0) = |y| − 2 × |pairs(y) ∩ pairs(y0)| − |unpaired(y) ∩ unpaired(y0)| (9)

• free energy gap ΔΔG(x,y), which is the difference between the free energies of y and MFE(x):

ΔΔG(x,y) = ΔG(x,y) − ΔG(x,MFE(x)) ≥ 0 (10)

If d(y,MFE(x))=0 or ΔΔG(x,y)=0, we have an MFE solution.

Ensemble-based criteria: Boltzmann probability and ensemble defect. The above two criteria only consider MFE structures, neglecting exponentially many competing structures. Even if y is the unique MFE structure, there could still be many highly competitive structures that are very close in energy to y (e.g., y* in Fig. 1b). As a result, the Boltzmann probability p(y|x) could still be arbitrarily small due to competition, which means x is highly unlikely to fold into y in equilibrium (Fig. 1c). So it is better to consider the whole ensemble to minimize the competition from alternative structures. We use two such metrics:

• Boltzmann probability p(y|x). Since each x has exponentially many possible structures in the ensemble, this probability can be arbitrarily small. For numerical stability, we minimize the negative log probability −log p(y|x).

• Ensemble defect, which is the expected structural distance between y and all structures in the ensemble36. This metric not only considers competition, but also how similar the competing structures are from y; we want to penalize those that are very different from y. The value of ensemble defect can be normalized to between 0 and 1, known as normalized ensemble defect (NED):

NED(x,y) = 1/|x| Ex0∼p(·|x) d(y,y0) = 1/|x| ∑ y0∈Y(x) p(y0|x) × d(y,y0) (11)

By Eq. (9) and some simplifications36,37, we get

NED(x,y) = 1 − 1/|x| (∑ (i,j)∈pairs(y) pij + ∑ j∈unpaired(y) qj) (12)

where pij is the base-pairing probability of nucleotides i and j, while qj = 1 − ∑ i pij is the probability of j being unpaired. NED(x,y) can also be derived as the mean of per-nucleotide positional defects, which quantify how incorrect a nucleotide folds in the ensemble (see Suppl. Data).

Now we formulate the RNA design problem as optimizing objective function f(x,y) over the design space X(y):

x* = argmin x∈X(y) f(x,y) (13)

where the objective function can be one of these four:

f(x,y) = {d(MFE(x),y) structural distance, ΔΔG(x,y) free energy gap, −log p(y|x) Boltzmann probability, NED(x,y) ensemble defect} (14)

Continuous optimization with coupled variable distributions
The most commonly used method for RNA design is local search (including random walk)4–9 which starts with a sequence and improves the objective function by modifying one or a few nucleotides per step. This is inadequate, however, for exploring the RNA design space X(y) which grows exponentially with the length of the target structure y (Eq. (2)).

As a result, the design quality degrades for longer and harder-to-design structures (see Results). Can we instead consider all candidate sequences simultaneously and promote the better ones? In other words, we prefer a global search method. Here we replace the discrete representation of a single candidate sequence by a probability distribution py(·) over all possible sequences x in X(y). Essentially, we propose the following continuous relaxation of the optimization problem minimizing the following objective function:

J = Ex∼py(·) [f(x,y)] = ∑ x∈X(y) py(x)f(x,y) (15)

This new objective function is to find a distribution py(·) of RNA candidates whose expectation of the objective function f(x,y) is minimized. If the probability mass concentrates on only one sequence, then this new relaxed objective degenerates to the original discrete objective.

Below we discuss how to model the probability distribution py(·) over the design space X(y), with increasingly more refined representations. Our novel distributions (v1–v3) capture important constraints and correlations in the design space, while still being easy to sample from.

Independent distributions (v0). The most obvious modeling of py(x) is to use the product of independent individual distributions p i(·) over N={A,C,G,U} for each position i (see Fig. 2(d)):

p0 y(x) = ∏ i pi(xi) (16)

This is the same distribution in previous work21. However, as explained earlier, the domain of this p0 y(·) distribution is the total sequence space N|y| which contains mostly invalid sequences, rather than the valid design space X(y) (Fig. 2(b–c)).

Coupled variables for pairs (v1). To model the distribution over the valid design space, we use a coupled variable for each base-pair:

p1 y(x) = ∏ i∈unpaired(y) pi(xi) × ∏ (i,j)∈pairs(y) pij(xixj) (17)

where pi(·) is the local distribution (over N) for unpaired position i, and pij(·) is the local distribution (over 6 choices in P) for paired positions (i,j) (see Fig. 2(e)). This distribution over the valid design space X(y) is a major contribution in our work.

Coupled variables for terminal mismatches (v2 and v3). Distributions v2 and v3 are refinements of the above v1. First, we note that in the standard energy models (Turner rules25), there are terminal mismatch lookup tables. For example, for a hairpin loop defined by the pair (i,j), the first and last nucleotide of the loop, xi+1 and xj−1, are the terminal mismatch, and will be looked up together in the energy table (such as x4 and x6 in Fig. 7(a)). Therefore, it is better to make a coupled variable over N2 (4 × 4 = 16 choices) for each terminal mismatch position-pair:

p2(x) = ∏ i∈unpaired0(y) pi(xi) × ∏ (i,j)∈pairs(y) pij(xixj) × ∏ (i,j)∈mismatches(y) pm(xixj) (18)

where unpaired0(y) = unpaired(y) \ mismatches(y). See Fig. 2(e).

Moreover, consider a special case like the 1-by-3 internal loop in Fig. 7(a). For such 1-by-x (x > 1) internal loops, there is exactly one unpaired nucleotide on one of the two branches, and that single unpaired nucleotide (x2 in our example) is included in two mismatches (x2 and x8 on one side and x2 and x10 on the other). Therefore, it is better to model all these three nucleotides together in a coupled variable over N3 (43=64 choices), which we call a “trimismatch”:

p3 y(x) = ∏ i∈unpaired00(y) pi(xi) × ∏ (i,j)∈pairs(y) pij(xixj) × ∏ (i,j)∈mismatches(y) pm(xixj) × ∏ (i,j,k)∈trimismatches(y) ptm(xixjxk) (19)

where unpaired00(y) = unpaired0(y) \ trimismatches(y). Note that, a 1-by-1 internal loop is a special case of mismatch (see Fig. 7(b), nucleotides 10 and 16).

In this work, we only considered mismatches and trimismatches for hairpins and internal loops. There are no mismatches for bulge loops, and we do not couple unpaired nucleotides in a multiloop or exterior loop; see Fig. 7 for details. Table S7 compares coupled-variable distributions (v1–v3) on Eterna100.

Sampling for objective evaluation, gradient estimation, and design space exploration
Given the complexity and diversity of RNA design problem settings, we aim at a method that can optimize any arbitrary objective function f(x,y). Even though we factorize py(·) into many tractable local distributions, without making any assumptions about the structure of f(x,y), the exact calculation of the expectation J=Ex∼py(·) [f(x,y)] is generally intractable. To address this challenge, we use sampling to approximate the objective, estimate its gradients, and explore integral solutions.

Sampling for Objective Evaluation. We approximate the expectation by averaging over samples S drawn from the distribution:

J = Ex∼py(·) [f(x,y)] ≈ 1/|S| ∑ x∈S f(x,y)

As |S| → ∞, this approximation will converge to the true expectation. In practice, even a small sample size might be sufficient (see Fig. 3(f)); our results use |S|=2500.

Sampling for Gradient Estimation. Next, we want to estimate the true gradient of J with respect to py(·):

∇ py J = ∇ py Ex∼py(·) [f(x,y)]

We again use sampling to estimate the above expectation by deriving the approximate gradient (see Supplementary Methods):

∇ py J(t) ≈ 1/|S| ∑ x∈S ∇ py log py(x)f(x,y) (20)

We then use the approximate gradient to update the distribution:

py(t+1) update py(t), ∇ py J(t) (21)

Sampling-based Design Space Exploration. At the end of this continuous optimization, we still need to return a single sequence from the distribution, i.e., an “integral solution”. This can be done by “rounding” if the distribution is close to one-hot, or more generally by taking the sequence with the highest probability in the final distribution. However, this max-probability sequence is not necessarily the best sequence in terms of the objective function, since the distribution is often not perfectly aligned with objective function. Here we use an alternative approach that simply takes the best sample in terms of the objective function out of all samples collected in the optimization process. This method outperforms the max-probability solution because the samples offer much more diversity in the exploration of the distribution than a single sequence (see Fig. 3(f)).

Parameterization and optimization
We describe two distinct ways to parameterize the distributions py(·) as py(·;Θ). The first method, direct parameterization, simply uses Θ as raw probabilities, but every update of Θ needs to result in probabilities leading to a harder constrained optimization problem. To solve this problem, we use Projected Gradient Descent (PGD). After each gradient step, we project Θ back onto the probability simplex that is closest (in ℓ2 norm) to Θ (details in Supplementary Methods).

The second method, softmax parameterization, models the distribution implicitly by applying softmax on Θ. The resulting problem is an unconstrained optimization problem and we can solve it with vanilla gradient descent (details in Supplementary Methods). We run the gradient decent step until the changes in the value of the objective function become sufficiently small, indicating convergence. See Algorithm 1 for the unified pseudocode for both constrained and unconstrained optimization approaches.

Algorithm 1. Sampling-based RNA Design
function DESIGN(y,f,projection=False) ⊳ f: objective
 x* ← random_init(y) ⊳ current best design
 Θ ← init_params(y)
 while not converged do
  sample sequences, S, from distribution py(·;Θ)
  x ← argmin x0∈S f(y0,y) ⊳ best sample in S
  if f(x,y) < f(x*,y) then x* ← x ⊳ better
  estimate objective J(Θ) and gradient ∇ Θ J(Θ) using S
  Θ ← Θ − α∇ Θ J(Θ)
  if projection then Θ ← proj(θ) | θ ∈ Θ ⊳ PGD
 return x*

Implementation details of SamplingDesign
By default, our method samples 2500 sequences at each step. The convergence condition is met when the 50-step moving average of the objective function (i.e., the mean over the last 50 steps) does not improve for 50 consecutive steps. We also have a hard limit of 2000 steps. Our main program is implemented in C++ and utilizes OpenMP for parallelization. We run our program on a high-performance computing server using 64 CPU cores (AMD EPYC 9554 64-Core Processor).

The default initial learning rate is set to 0.01, which works well for all puzzles under the softmax parameterization. However, we observed that the direct parameterization (projection) requires a smaller learning rate as the puzzle length increases; otherwise, the objective value does not improve. Therefore, we apply adaptive learning rate decay for the projection method, which reduces the learning rate by half when the objective function has not improved for 20 steps. We also use momentum-based optimizers to speed up learning: Adam38 for softmax and Nesterov accelerated gradient39 for projection.

We use three types of initialization:
• Uniform: uniform sequence distribution.
• Targeted: 100% A for unpaired, 50:50 C/G for paired, 50:50 G/U for mismatches, from SAMFEO’s initialization4.
• ϵ-Targeted: A combination of targeted and uniform distributions defined by ϵ·targeted+(1−ϵ)·uniform.

For the projection method, we use both uniform and targeted initializations. For the softmax parameterization (our default method), we use uniform and ϵ-targeted initializations with ϵ=0.75. The final solutions are selected from both initializations. There are 14 puzzles where uniform initialization yields higher Boltzmann probability than ϵ-targeted on the Eterna100 dataset (see also Table S8). For example, puzzle #74 improves from 0.202 to 0.457, and puzzle #77 improves from 0.334 to 0.400.

Our main results use the softmax parameterization. Although using the direct parameterization achieves similar results in terms of the arithmetic mean of Boltzmann probability (Table S9), it substantially underperforms in terms of the geometric mean (0.104 vs. 0.502). This aligns with our observation that using projection requires a smaller learning rate for longer puzzles, resulting in much slower optimization.

We use our previous work, LinearPartition37, to compute the partition function and base-pairing probabilities in linear time with beam pruning. We use beam size b=250 for optimizing p(y|x), and b=100 for optimizing NED(x,y). We use a larger beam size for p(y|x) because LinearPartition tends to underapproximate the partition function, resulting in p(y|x)>1 occasionally. For optimizing d(MFE(x),y) and ΔΔG(x,y), we find the MFE(x) using LinearFold40 with b=100.

At each step, we save the best among the 2,500 samples. After completing all steps, we reevaluate those samples using ViennaRNA v2.6.47, and output the best sequence for each metric. The values reported in this paper are thus unaffected by the approximation error from beam search.

Baseline RNA design methods
Baseline 1: SAMFEO. It is a random walk method that selects nucleotides to mutate based on positional defects (see Supplementary Methods) and structural information4. It achieved the strongest results in the literature across multiple metrics such as p(y|x), outperforming many competitors including RNAinverse7, NUPACK 4.012, and MODENA41. Like our approach, it can optimize any arbitrary objective function. We ran SAMFEO five times using their default settings, and report the best solution among them.

Baseline 2: NEMO. It is another local search method that applies Nested Monte Carlo Search and uses heuristics tailored to the energy model5. It is the strongest baseline in the literature in terms of MFE-based criteria. It maximizes:

score(x,y) = K(1+ΔΔG(x,y)) − sign(K) |pairs(y) ∩ pairs(y0)| (22)

We ran NEMO five times and report the best solutions.

Baseline 3: Matthies et al.21. Although it also uses continuous optimization for RNA design, our approach is vastly different from and substantially outperforms theirs in both scalability and design quality (by all metrics).
• First, their sequence distribution is a simple product of independent distributions for each position (same as our distribution v0). Instead, we use coupled variables for paired and mismatch positions (our distributions v1–v3), which rules out invalid sequences and explicitly models the dependencies between correlated positions.
• Second, our sampling framework can work with arbitrary objective functions while their work is specifically designed the Boltzmann probability.
• Third, our sampling yields an unbiased approximation to the expectation of an arbitrary objective function over the distribution of sequences. For example, for the case of Boltzmann probability, our sampling approximates the expected Boltzmann probability, which converges to the true expectation as sample size increases:

1/|S| ∑ x∈S p(y|x) = Ex∼py(·) [p(y|x)] (23)

By contrast, they optimize a different objective (in red below) that deviates from the true expectation of Boltzmann probability with a bias (E[X/Y] ≠ E[X]/E[Y]):

Ex [e−ΔG(x,y)/RT] / Ex [Q(x)] ≠ Ex [e−ΔG(x,y)/RT / Q(x)] = Δ Ex [p(y|x)] (24)

Deep learning baselines. There is also work using deep reinforcement learning for RNA design42,43. However, their results are inferior to the above baselines (see Table S10).

RNA secondary structure datasets
To evaluate SamplingDesign against other RNA design methods, we use four pseudoknot-free secondary structure datasets:
1. Eterna10044 is the standard benchmark for evaluating RNA design programs. It contains 100 secondary structures (i.e., “puzzles”) of up to 400 nts, varying in design difficulty from simple hairpins to intricate multiloop structures.
2. Eterna100-v245 modifies 19 puzzles from Eterna100 such that all puzzles are (MFE) solvable under the Turner 2004 energy model25 (using ViennaRNA v27). This dataset also contains solutions from Eterna players for all puzzles except for puzzle #100 (whose solution is held out).
3. Rfam-Taneda-2741 includes consensus secondary structures from seed alignments in Rfam 9.0, with each entry corresponding to an RNA family.
4. RNAsolo-76446,47 is a collection of secondary structures curated from a database of experimentally determined RNA 3D structures. See Fig. S8 for the statistics of these datasets.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript presents an innovative Monte Carlo sampling-based method for RNA inverse folding, integrating joint sequence distribution optimization with the objective function to enhance efficiency. While the core concept of probabilistic optimization via gradient updates is novel and addresses computational bottlenecks, the study currently relies on a single dataset (Eterna100) for validation. Broader experimental validation across diverse RNA structure types and datasets would significantly strengthen the claims regarding generalizability and practical impact.\nMajor:\n1. On page 6, in \"Comparison with NEMO,\" it's mentioned that NEMO's effectiveness relies on heuristic rules manually summarized. Is there any quantification or demonstration of the improvement from these rules? Also, can SamplingDesign use these rules to enhance itself?\n2.On page 10, Fig. 3(g) is mentioned twice, but it is absent in the article.\n3.As far as I know, the reason why Eterna 100 has some undesignable RNA puzzles is that RNAFold, which is based on Turner 1999 parameters, was used to assess the designs. Now, RNAFold has adopted Turner 2004 parameters. Have you tried using RNAFold with Turner 1999 parameters to design Eterna 100 puzzles? Meanwhile, there is also a dataset specifically improved for the RNAFold with Turner 2004 parameters, i.e., Eterna 100 v2 (https://github.com/eternagame/eterna100-benchmarking). This dataset has modified 19 puzzles that were previously undesignable in Eterna 100 v1, ensuring feasible sequences exist under the evaluation of RNAFold with Turner 2004 parameters. Would you consider conducting comparative experiments and analysis on Eterna 100 v2?\n4.Have you thought about optimizing RNA distribution probability by combining multiple objective functions?\nMinor:\n1.On page 3, in the last sentence of the \"Independent vs. Coupled - Variable Distributions.\" section, it's better to introduce Figure 7 earlier to explain mismatch and tri-mismatch, easing readers' initial reading burden.\n2.On page 9, in the second paragraph of \"Continuous optimization with coupled variable distributions,\" there's an improper line break after \"Here.\"", "reviewer_2": "The paper presents a sampling-based continuous optimization method (SamplingDesign) for RNA design problems and demonstrates its superiority through the Eterna100 benchmark. The method outperforms existing state-of-the-art approaches on several key metrics, particularly for long and complex RNA structures. The idea and practicality of the paper are commendable. However, there is still room for improvement in the experimental setup, theoretical exposition, and evaluation methods.\n1.The experimental setup relies solely on the Eterna100 dataset, which lacks diversity. Additional datasets with more complex RNA structures and real-world biological RNA molecules should be included to better validate the method's applicability. The comparison with other methods should also encompass emerging deep learning-based RNA design approaches for a more comprehensive benchmark.\n2.The theoretical exposition of the coupled-variable model is insufficient. A more detailed explanation, supported by mathematical formulas and diagrams, is needed to clarify how the model is constructed and updated. The optimization process, particularly the gradient estimation via Monte Carlo sampling, requires a more in-depth mathematical derivation to enhance the paper's rigor.\n3.The evaluation metrics could benefit from the inclusion of additional criteria related to RNA functionality, such as thermodynamic stability and kinetic properties. Statistical analysis to demonstrate the significance of differences between the proposed method and other methods would strengthen the findings. The results could be more effectively communicated through additional visualizations, such as 3D structure diagrams of the designed RNA molecules and detailed energy landscape plots.\n4.The presentation of results lacks sufficient visual aids. Including 3D structure diagrams of the designed RNA molecules and detailed energy landscape plots would help illustrate the method's performance more intuitively. A sensitivity analysis of key parameters, presented visually, would provide valuable insights into the robustness and tuning of the method.\n5.The discussion section should more thoroughly address the limitations of the current approach, such as potential computational costs associated with the coupled-variable model and the accuracy of gradient estimation. A more detailed outlook on future work, particularly potential improvements to the model and integration with other biological constraints, would offer a clearer direction for subsequent research.\n6.The paper does not adequately address the scalability and computational efficiency of the proposed method. Given the complexity of RNA structures, especially longer sequences, it is essential to discuss how the method performs in terms of runtime and memory usage. Including a section on computational efficiency and potential optimizations would be beneficial.\n7.While the method shows promise in computational benchmarks, its biological relevance is not fully explored. The paper should include a discussion on how the designed RNA sequences could be validated experimentally and their potential applications in real-world biological systems. This would bridge the gap between computational design and practical utility.\n8.The results are presented in terms of several metrics such as Boltzmann probability and NED. Are these metrics equally important for all applications of RNA design, or do different applications place different emphasis on certain metrics? How does this impact the overall evaluation of the method's effectiveness?\n9.The evaluation focuses on the ability of the method to design RNA sequences that fold into target structures. Are there any evaluations of the biological feasibility or functionality of the designed sequences? For example, have any experiments been conducted to test the designed RNA molecules in a biological context?", "reviewer_3": "Tang et al. present the method SamplingDesign, which combines several computational procedures into a compelling RNA design method. Namely, it uses sampling of sequences and their associated structures to estimate gradients as well as distributions that explicitly model pairs, mismatches, tri-mismatches, and unpaired nucleotides. The combination of these innovations greatly helps to reduce the possible design search space. This manuscript’s Supplementary Methods are done well and helps with understanding SamplingDesign’s algorithm. The authors do a good job benchmarking their new method, which performs well compared to the previous methods SAMFEO, NEMO, and Matthies et al. Benefits of the SamplingDesign method include greater flexibility of objective function potentially for different applications, ability to design longer RNA structures, and better performance in many metrics on the Eterna100 dataset.\nMajor points:\n1. Authors should make clear that SamplingDesign is meant for RNA secondary structure design in its current form.\n2. A little more clarification under “Evaluation setup” #6 and #7. Does #6 mean the number of global minimum free energy structures, and #7 mean only counting if there is exactly 1 global minimum free energy structure? It is more precisely stated around Eqn. 2 and 8 in the Methods but would be good to help the reader more in the main text.\n3. Fig. 6:\n- caption starts with panel b and not a\n- Panel b: pink coloring is used to denote undesignable region, but the color is not in caption.\n4. Item #2 in Discussion discusses how they test Matthies et al. on GPUs and ran into GPU memory issues with longer RNA lengths. The authors do not mention what resources they used for running SamplingDesign. Was SamplingDesign also run on an 80GB NVIDIA H100 GPU, and for this type of resource can SamplingDesign reasonably jump to designing mRNAs that can be several kb long?\n5. Item #4 in Discussion: The distribution categories used within SamplingDesign seem designed for RNA (pair, unpaired, mismatch, tri-mismatch) so the jump to protein design should be further qualified since these four categories are not directly applicable to protein.\n6. How are asymmetric bulges with only paired nts on one side of the helix encoded in the distribution versions? Ex. 1 nt bulge. Similarly in Fig.7 right side, the unpaired nts in the 3-way junction are labeled as unpaired even though they could be interpreted as a “mismatch”? Please provide more description.\n7. The description around the possible design space and quantifying this design space got confusing through the manuscript. The authors’ major innovation is reducing the search path through their choice in distribution (i.e. base pairs, mismatch, trimismatch, unpaired), but they simultaneously say that the possible design space is 4N in Fig. 1 and that 4N is incorrect in other parts of the manuscript (ex. “RNA design space: full vs. valid”). Please go back through manuscript and make these distinctions clearer. Placing equation 6 earlier in the manuscript may help.\n8. Introduction:\n“However, such methods are only capable of exploring a fixed number of candidate sequences, thus not able to keep up with the exponential growth of the design space.”\nHere, the authors could clarify the rationale for their method vs. local search. The wording might imply that SamplingDesign can search an exponential number of sequences, whereas local search only explores a fixed number. It should be made clear that SamplingDesign’s discrete sampling allows the method to escape local minima and explore more of the search space as the number of samples increase.\nSince the objective function is arbitrary, this is an advantage over local search methods that reach a local minimum when the objective function is non-convex with regards to the search space.\nMinor points:\n1. Fig. 4d: ‘*’ is not described in caption or main text.\n2. In “Softmax parameterization and unconstrained optimization” between equations 31 and 32, there is mixed usage of ‘S’ and ‘Z’.\n3. Is random initialization of SamplingDesign possible? The authors use 3 other types of initialization methods and curious about random initialization because it is used frequently in ML.\n4. Would it be possible within the framework of SamplingDesign to vary sample size between steps? Earlier steps may benefit more from larger sample size to estimate the gradient and promote exploration.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to Reviews\nDear Reviewers,\nThank you very much for the detailed feedback. We have revised our manuscript accordingly. In this\ndocument, we detail all the changes, where reviewer feedback is provided verbatim in black and our\nresponses in blue. Major revisions in the manuscript are shown in quoted italic text.\nSincerely,\nLiang Huang (corresponding author)\nQuick Links to:\nReviewer #1\nReviewer #2\nReviewer #3\nOther changes\nReviewer #1 (Remarks to the Author)\nThe manuscript presents an innovative Monte Carlo sampling-based method for RNA inverse folding,\nintegrating joint sequence distribution optimization with the objective function to enhance efficiency.\nWhile the core concept of probabilistic optimization via gradient updates is novel and addresses\ncomputational bottlenecks, the study currently relies on a single dataset (Eterna100) for validation.\nBroader experimental validation across diverse RNA structure types and datasets would significantly\nstrengthen the claims regarding generalizability and practical impact.\nMajor:\n1. On page 6, in \"Comparison with NEMO,\" it's mentioned that NEMO's effectiveness relies on heuristic\nrules manually summarized. Is there any quantification or demonstration of the improvement from these\nrules? Also, can SamplingDesign use these rules to enhance itself?\nThank you for the suggestion. NEMO heuristics are handwritten rules based on the experience of Eterna\nplayers, including:\n1. Sample base pairs with bias (GC > AU > GU);\n2. Use “blocking” nucleotides to prevent unwanted base-pairing in triloops and bulges;\n3. Favor specific mismatches by loop type (e.g., G/G in 1×1 internal loops).\n(a) Re: “quantification of improvement from NEMO’s rules.” In the main results, we stated the\nnumber(s) of Eterna100 puzzles solved by NEMO without the heuristic rules: “Previous work (4)\nshowed that NEMO without those rules solves only 76 and 75 puzzles under the MFE and uMFE\ncriteria, respectively.” Basically, applying these heuristics increases NEMO’s solved counts by 3\n(MFE) and 2 (uMFE), respectively. We added the previous sentence to the manuscript.\n(b) Re: “Can SamplingDesign use these rules to enhance itself?” Applying NEMO’s rules to\nSamplingDesign is non-trivial because (1) they are tailored to the MFE criteria, and (2) NEMO\nhard-codes them in their Monte Carlo playout policies, whereas SamplingDesign is a generic\nalgorithm that can optimize any arbitrary objective function. That being said, in our ε-targeted\ninitialization, we also favor GC base pairs and unpaired A nucleotides. We could further\nincorporate their loop-specific mismatch priors in place of our current uniform choice for\nmismatches in future work.\n2. On page 10, Fig. 3(g) is mentioned twice, but it is absent in the article.\nThank you for catching this. Both references should point to Fig. 3(f), not 3(g). We have corrected the two\ninstances on page 10.\n3. As far as I know, the reason why Eterna 100 has some undesignable RNA puzzles is that RNAFold,\nwhich is based on Turner 1999 parameters, was used to assess the designs. Now, RNAFold has adopted\nTurner 2004 parameters. Have you tried using RNAFold with Turner 1999 parameters to design Eterna\n100 puzzles? Meanwhile, there is also a dataset specifically improved for the RNAFold with Turner 2004\nparameters, i.e., Eterna 100 v2 (https://github.com/eternagame/eterna100-benchmarking). This dataset has\nmodified 19 puzzles that were previously undesignable in Eterna 100 v1, ensuring feasible sequences\nexist under the evaluation of RNAFold with Turner 2004 parameters. Would you consider conducting\ncomparative experiments and analysis on Eterna 100 v2?\nThank you for the thoughtful comment. We use the Turner 2004 parameters as implemented in Vienna\nRNAfold v2.x and LinearFold/LinearPartition (-V mode); this set of parameters is the most commonly\nused energy model in recent RNA secondary structure prediction and RNA design work [Runge et al.,\nBioinformatics 2024; Matthies et al., NAR 2023; Eastman et al., PLOS 2018]. We have added the previous sentence to\nthe manuscript.\nThe Eterna100 v2 dataset modifies 19 puzzles from v1 so that these puzzles are MFE-solvable under the\nTurner 2004 parameters. The v2 dataset provides human player solutions for 18 of them, while the\nsolution for puzzle #100 (“Teslagon”) is undisclosed. We verified that all these 18 solutions fold to their\ntarget structures using Vienna RNAfold v2.6.4. However, only 13 of the 18 solutions satisfy the unique\nMFE (uMFE) criterion, meaning that some of these puzzles may still be undesignable under the uMFE\ncriterion we used.\nWe ran SamplingDesign and baselines on the Eterna100 v2 and report the results in Table S1 (reproduced\nbelow). With the modified puzzles, SamplingDesign still achieves the strongest performance in the\nensemble-based metrics, which are more important than MFE-based criteria (see Fig. 1(a) and Reviewer\n#2, comment 8). It obtains the highest Boltzmann probability (both arithmetic and geometric means) and\nthe lowest ensemble defect.\n4. Have you thought about optimizing RNA distribution probability by combining multiple objective\nfunctions?\nThank you for the suggestion. Yes, SamplingDesign is capable of optimizing a composite objective, such\nas the weighted sum of Boltzmann probability and ensemble defect. However, our main results on\nEterna100 show that optimizing Boltzmann probability alone outperforms optimizing other metrics\ndirectly (e.g., normalized ensemble defect and free energy gap). Furthermore, optimizing a single\nobjective avoids the cost of evaluating additional objectives. For these reasons, we focus on optimizing\nBoltzmann probability in this work, but optimizing combined objectives remains a promising direction for\nfuture work.\nMinor:\n1. On page 3, in the last sentence of the \"Independent vs. Coupled - Variable Distributions.\" section, it's\nbetter to introduce Figure 7 earlier to explain mismatch and tri-mismatch, easing readers' initial reading\nburden.\nThank you for the helpful suggestion. However, we think that moving the entire Fig. 7 forward would add\ntoo much method-specific detail to the Introduction. The hairpin stem-loop in Fig. 2(a) already introduces\nthe three most commonly used position types (unpaired, pairs, and mismatches), while trimismatches only\noccur in the specific cases of 1 x n internal loops with n > 1. Instead, we added a forward reference to Fig.\n7 in the Results subsection “Independent vs. Coupled-Variable Distributions” (reproduced below).\n“Moreover, in the standard RNA energy model (Turner rules (23)), there are many other correlated\npositions, which we classify as mismatches such as (4,6) in Fig. 2(a) and trimismatches (Fig. 7 in\nMethods).”\n2. On page 9, in the second paragraph of \"Continuous optimization with coupled variable distributions,\"\nthere's an improper line break after \"Here.\"\nThank you for catching this. We have removed the line break in the manuscript.\n(Remarks on code availability)\nThe code is clearly documented and reproducible. The README provides sufficient instructions.\nReviewer #2 (Remarks to the Author)\nThe paper presents a sampling-based continuous optimization method (SamplingDesign) for RNA design\nproblems and demonstrates its superiority through the Eterna100 benchmark. The method outperforms\nexisting state-of-the-art approaches on several key metrics, particularly for long and complex RNA\nstructures. The idea and practicality of the paper are commendable. However, there is still room for\nimprovement in the experimental setup, theoretical exposition, and evaluation methods.\n1. The experimental setup relies solely on the Eterna100 dataset, which lacks diversity. Additional\ndatasets with more complex RNA structures and real-world biological RNA molecules should be included\nto better validate the method's applicability. The comparison with other methods should also encompass\nemerging deep learning-based RNA design approaches for a more comprehensive benchmark.\nThank you for the suggestion.\n(a) Re: “Additional datasets with more complex RNA structures and real-world biological RNA\nmolecules should be included.” We expanded our evaluation beyond Eterna100 and now compare\nSamplingDesign to baseline methods on three additional RNA secondary structure datasets:\nEterna100 v2 (see also Reviewer #1, comment 3), Rfam-Taneda-27, and RNAsolo-764.\nEterna100 (v1/v2) datasets are specifically constructed to include particularly challenging,\ncomplex structures, whereas both Rfam-Taneda-27 and RNAsolo-764 datasets comprise\nsecondary structures derived from naturally occurring RNAs. Rfam-Taneda secondary structures\nare determined by comparative sequence analysis, while RNAsolo secondary structures are\nderived from experimentally determined 3D RNA structures [Badura et al., Scientific Reports, 2025].\nIn this revision, we added a new Methods subsection, “RNA secondary structure datasets,” and\nFig. S8, which shows the data curation process and histograms of secondary structure lengths for\neach dataset. Fig. 4(a) and Table S1–S3 report the results of SamplingDesign and the baselines on\neach dataset. The results show that SamplingDesign consistently outperforms other baselines on\nthe more important ensemble-based metrics (i.e., Boltzmann probability and NED).\n(b) Re: “Deep learning-based RNA design approaches.” For deep learning-based approaches, there\nare EternaBrain [Koodli et al., PLOS, 2019], LEARNA and its extension Meta-LEARNA [Runge et al.,\nICLR, 2019], and Eastman et al. [PLOS, 2018]. EternaBrain uses a convolutional neural network\n(CNN) trained on Eterna players’ moves, but its training and evaluation rely on the outdated\nTurner 1999 parameters (Vienna 1.8.5). LEARNA/Meta-LEARNA and Eastman et al. use\nreinforcement learning to design sequences for a target secondary structure; however, their\nreported Eterna100 results (68/100 and 60/100 solved under the MFE criteria, respectively;\nreproduced below) remain below the local search-based baselines that we used, such as NEMO\n(79/100). Since we were unable to reproduce their results using their code, in Table S10, we only\nreported the number of MFE-solved puzzles from their publication and re-evaluated their\nsolutions (if provided) for the number of uMFE solutions.\n2. The theoretical exposition of the coupled-variable model is insufficient. A more detailed explanation,\nsupported by mathematical formulas and diagrams, is needed to clarify how the model is constructed and\nupdated. The optimization process, particularly the gradient estimation via Monte Carlo sampling,\nrequires a more in-depth mathematical derivation to enhance the paper's rigor.\nThank you for the suggestion.\n(a) Re: “The theoretical exposition of the coupled-variable model is insufficient.” For theoretical\nexposition, we already have Eqs. (16)-(19) in the Methods, which gives the precise factorization\nof the distribution: starting from the fully independent model (v0) and progressively introducing\nmore position types (base-pairs, mismatches, and trimismatches), resulting in Eq. (19), which\ndefines the coupled-variable distribution (v3). For diagrams, Fig. 2 motivates why a\ncoupled-variable distribution is needed over a fully independent distribution, and Fig. 7 in\nMethods shows examples of how nucleotides are classified into one of the position types in a\ncoupled-variable distribution. Finally, Fig. 3 and Alg. 1 in Methods explain how the model is\nupdated (i.e., the optimization process) with mathematical detail provided in the Supplementary\nMethods.\n(b) Re: “The optimization process requires a more in-depth mathematical derivation.” For readers\ninterested in the mathematical detail, we refer them to the Supplementary Methods. For example,\nin the Methods subsection “Sampling for Gradient Estimation”, we stated: “We again use\nsampling to estimate the above expectation by deriving the approximate gradient (see\nSupplementary Methods).” The Supplementary Methods provide the full derivation of Monte\nCarlo gradient estimation, including the final expressions for the direct parameterization (Eq. 31)\nand the softmax parameterization (Eq. 36).\nWe hope these pointers address the request for a more rigorous exposition without duplicating material\nalready present in the paper.\n3. The evaluation metrics could benefit from the inclusion of additional criteria related to RNA functionality,\nsuch as thermodynamic stability and kinetic properties. Statistical analysis to demonstrate\nthe significance of differences between the proposed method and other methods would strengthen the\nfindings. The results could be more effectively communicated through additional visualizations, such as\n3D structure diagrams of the designed RNA molecules and detailed energy landscape plots.\nThank you for the suggestion.\n(a) Re: “Inclusion of additional criteria related to RNA functionality.” We do not evaluate\nthermodynamic stability or kinetic properties in this work. First, the thermodynamic stability (i.e.,\n∆𝐺(𝑥,𝑦)) alone is not a reliable RNA design metric because it ignores competition from the\nexponentially many alternative structures (see Fig. 1 and Methods subsection “Ensemble-based\ncriteria: Boltzmann probability and ensemble defect”). In addition, optimizing ∆𝐺(𝑥,𝑦) leads to\ndesigns with many alternative structures with similar free energy; see [Yamagami et al., NAR, 2019, Fig.\n2C] (reproduced below), where low-energy designs have high NEDs (top left corner of the plot).\nSecond, kinetic properties are more important for designing RNA switches, whereas our work\nfocuses on nano-structure design. Therefore, we choose the thermodynamics model, and\nevaluating kinetics is outside the scope of this work.\n(b) Re: “Statistical analysis.” We performed one-sided Wilcoxon signed rank tests between\nSamplingDesign and other baselines for Boltzmann probability across all datasets used. We\nreported the p-values that are significant in the captions (Fig. 4(a) and Table S1 – S3) and\ndescribed the test results in the Results subsection “SamplingDesign outperforms state-of-the-art\nmethods” (reproduced below):\n“Across all datasets (Fig. 4(a) and Tables S2–S4), SamplingDesign consistently performs the best\nin\nterms of the more important ensemble-based metrics: Boltzmann probability and NED. Against\nSAMFEO, SamplingDesign achieves significantly higher Boltzmann probability on Eterna100 (p\n= 0.009) and Eterna100 v2 (p = 0.0007) by a one-sided Wilcoxon signed rank test; the\ndifferences against NEMO are even more significant on all datasets.”\n(c) Re: “Additional visualizations.” This is addressed in the next comment (Reviewer #2, comment\n4).\n4. The presentation of results lacks sufficient visual aids. Including 3D structure diagrams of the designed\nRNA molecules and detailed energy landscape plots would help illustrate the method's performance more\nintuitively. A sensitivity analysis of key parameters, presented visually, would provide valuable insights\ninto the robustness and tuning of the method.\nThank you for the helpful suggestions.\n(a) Re: “3D structure diagrams.” To clarify, our work is about secondary structure-based RNA\ndesign, which is why we did not include 3D structure diagrams. We have revised the abstract and\nintroduction to make that clear (reproduced below, see also Reviewer #3, comment 1).\n“RNA design aims to find an RNA sequence that can fold into a given target secondary structure, ...”\nThe functional relevance of RNA secondary structure is well documented (e.g., hammerhead and\ndouble-pseudoknotted ribozymes) [Dotu et al., NAR, 2014; Yamagami et al., NAR, 2019], and our metrics\nare defined at the secondary structure level, not the tertiary level. For these reasons, we focused\nour figures on secondary structure diagrams and base-pairing probabilities, rather than tertiary\n(3D) structure diagrams. Adding predicted 3D models could overstate claims unrelated to our\nscope. In addition, improving the prediction accuracy of 3D models from secondary structure is\nstill an active research problem.\n(b) Re: “Energy landscape plots.” We did not include energy landscape plots, but our results in the\nensemble-based metrics (e.g., Boltzmann probability) already reflect the energy landscapes. For\nexample, Fig. 1 shows that obtaining a high Boltzmann probability requires designing sequences\nthat sharpen the energy landscape.\n(c) Re: “Sensitivity analysis of key parameters.” SamplingDesign has several tunable options,\nincluding the choice of parameterization, initializations, and the version of coupled-variable\ndistributions (v1–v3). We evaluated these choices on Eterna100 (Tables S7–S9) and found that\nusing softmax parameterization, coupled-variable distribution v3, with combined initializations\n(both uniform and ε-targeted), yields the optimal results.\nIn addition, we now visualize the effect of three hyperparameters: sample size, learning rate, and\nε (in ε-targeted initialization) by plotting learning curves for Puzzle #49 under different settings\n(Fig. S13, reproduced below). We also added a new paragraph to the Supplementary Methods\nsubsection, “Learning curves,” to discuss the effect of tuning these hyperparameters.\nCaption: “Fig. S13. SamplingDesign learning curves on Eterna puzzle #49 (101 nt) under different\nhyperparameters. Learning curves comparing (a) learning rate α, (b) sample size per step |𝑆|, and (c)\nϵ-targeted initializations. Higher ϵ favors GC base pairs and unpaired A nucleotides during initialization;\nϵ = 0 corresponds to uniform initialization. Settings highlighted in blue text indicate the default settings.\nThe target secondary structure (puzzle #49) is shown at the bottom.”\n5. The discussion section should more thoroughly address the limitations of the current approach, such as\npotential computational costs associated with the coupled-variable model and the accuracy of gradient\nestimation. A more detailed outlook on future work, particularly potential improvements to the model and\nintegration with other biological constraints, would offer a clearer direction for subsequent research.\n(a) Re: “Limitations of the current approach.” First, the computational cost associated with the\ncoupled-variable distribution is small because our optimization uses a sampling framework and\nwe couple only a small number of positions in a single variable (a trimismatch has 4 = 64\nnucleotide combinations). We found that using coupled variables for correlated positions with\nlong-range dependencies (i.e., base pair, mismatch, and trimismatch) introduced negligible\ncomputational cost and improved the results substantially (Table S7, reproduced below).\nSecond, the accuracy of gradient estimation depends on the number of samples. In Fig. S13(b)\n(reproduced below), we show how different sample sizes affect the learning curve for a target\nstructure with 101 nucleotides. Even with only 50 samples per step, the Boltzmann probability of\nthe best sample (green circle) continues to improve across steps. However, using a larger sample\nsize results in faster convergence and a final design with a higher Boltzmann probability.\nThus, the main limitation of SamplingDesign lies in its sampling efficiency. In Discussion bullet\n2, we suggested two ways to improve the sampling efficiency by caching subsequences and using\nimportance sampling.\n(b) Re: “Potential improvements to the model.” In the Discussion, we now provide four bullet points\nfor ways that our coupled-variable models could be improved or extended. In bullet 1, we propose\nusing richer distributions, such as discriminative or neural models, though sampling from these\nwould be non-trivial. In bullet 2, we suggest improving sampling efficiency through subsequence\ncaching and importance sampling (as noted in point (a)). Bullets 3 and 4 describe how the\nframework could be applied to related design problems, such as mRNA design (see Reviewer #3,\ncomment 4) and partial RNA design (discussed in point (c)).\n(c) Re: “Integration with other biological constraints.” One commonly used biological constraint in\nRNA design is nucleotide constraints at specific positions. These constraints can be specified\nusing the IUPAC nucleotide code (e.g., ‘Y’ means either C or U), and have been applied in\nprevious RNA design studies with experimental validation [Dotu et al., NAR 2014; Yamagami et al., NAR\n2018]. These constraints are often derived from comparative analyses of natural sequences that\nidentify highly conserved nucleotides in sequence alignments or from well-characterized sites\nimportant for RNA function. Incorporating such nucleotide constraints effectively reduces the\nsize of the design space. Implementing SamplingDesign’s coupled-variable distribution with\nnucleotide constraints is definitely feasible and will be addressed in future work. More generally,\nwe can extend SamplingDesign to the “partial RNA design problem” proposed by Runge et al.\n[NAR 2024], where both the sequence and target structure can be partially specified. We have\nadded this extension as bullet 4 in the Discussion.\nSamplingDesign can also be extended to the related mRNA design problem (see Discussion,\nbullet 2), which takes a protein sequence 𝑝 as input and tries to find an mRNA sequence 𝑥\nencoding 𝑝 that optimizes some objective function 𝑓(𝑥). This extension requires changing the\nobjective function and modifying the coupled-variable distribution. This is particularly suitable\nfor SamplingDesign because it is capable of optimizing any arbitrary objective function.\n6. The paper does not adequately address the scalability and computational efficiency of the proposed\nmethod. Given the complexity of RNA structures, especially longer sequences, it is essential to discuss\nhow the method performs in terms of runtime and memory usage. Including a section on computational\nefficiency and potential optimizations would be beneficial.\nThank you for the comment. We previously addressed scalability and computational efficiency in the\nSupplementary Methods subsection “Time analysis.” We have revised both the subsection and Figure S10\n(reproduced below) to address the runtime analysis more thoroughly, including the resources on which\nSamplingDesign was run and the total runtime under different initializations.\nCaption: “Fig. S10. Computational cost on Eterna100 puzzles. (a) Runtime per puzzle on Eterna100 using\nuniform and ϵ-targeted initializations, and the combined runtime using both initializations (sum of\nruntimes). (b) Total number of steps per puzzle and the step at which the best Boltzmann probability p(y |\nx) solution was found. (c) Peak memory (RSS) usage on Eterna100. All computations were performed on\na high-performance computing server using 64 CPU cores (AMD EPYC 9554 64-Core Processor).”\n7. While the method shows promise in computational benchmarks, its biological relevance is not fully\nexplored. The paper should include a discussion on how the designed RNA sequences could be validated\nexperimentally and their potential applications in real-world biological systems. This would bridge the\ngap between computational design and practical utility.\nThe secondary structure of an RNA molecule can be determined by chemical probing [Deigan et al. PNAS,\n2009]. The idea is that unpaired nucleotides are generally more accessible and reactive than paired\nnucleotides. First, a chemical reagent (e.g., dimethyl sulfate) is used to modify specific bases in the RNA.\nThese modifications are then detected through reverse transcription and sequencing, producing a\nreactivity profile for each nucleotide. Finally, the resulting reactivity profiles are integrated with\nthermodynamic models to infer base pairs and unpaired nucleotides.\nBesides that, RNAs adopt specific structures to carry out their biological functions. As such, the designed\nRNAs can be validated experimentally by examining their activity. For example, previous studies have\nsynthesized self-cleaving ribozymes by designing RNA sequences that fold into known ribozyme\nstructures [Dotu et al., NAR 2014; Yamagami et al., NAR 2018]. In these studies, the designed RNA molecules are\nobtained by in vitro transcription from DNA templates preceded by a T7 RNA polymerase promoter.\nThen, the activities of these RNA molecules are examined through self-cleaving assays.\nWe added a bullet in Discussion (reproduced below), bullet 5, for experimental validation and referenced\nthese studies.\n“While SamplingDesign excels in in silico benchmarks, we plan to validate our designs experimentally in\nfuture work, such as designing self-cleaving ribozymes (13, 14) or using chemical probing to determine\nsecondary structures (33).”\n8. The results are presented in terms of several metrics such as Boltzmann probability and NED. Are these\nmetrics equally important for all applications of RNA design, or do different applications place different\nemphasis on certain metrics? How does this impact the overall evaluation of the method's effectiveness?\nAs discussed in the Introduction and Fig. 1, we believe these metrics are not equally important. While\nmost of the RNA design literature focuses on the MFE metric, it is not sufficient because it ignores the\ncompetition from alternative structures in the ensemble. On the other hand, ensemble-based metrics such\nas the Boltzmann probability and NED are more comprehensive measures of the design quality.\nInterestingly, we found that optimizing for Boltzmann probability also resulted in the best performances\nin other metrics. Incidentally, another paper [Ward et al., NAR. 2023], which we have added in this revision,\nshares our observation that optimizing Boltzmann probability also optimizes other metrics as byproducts\n(reproduced below).\n“Notably, optimizing Boltzmann probability yields the best results across all other metrics, even better\nthan optimizing those metrics directly. For example, when optimizing for Boltzmann probability,\nSamplingDesign achieves the lowest average NED of 0.035, outperforming both itself and SAMFEO\noptimizing for NED (0.041 and 0.037, resp.). The fact that optimizing Boltzmann probability also\noptimizes other metrics is consistent with the observations of Ward et al. (26).”\n9. The evaluation focuses on the ability of the method to design RNA sequences that fold into target\nstructures. Are there any evaluations of the biological feasibility or functionality of the designed\nsequences? For example, have any experiments been conducted to test the designed RNA molecules in a\nbiological context?\nThank you for the comment.\n1. We note that the vast majority of RNA design papers do not have experimental validation [Busch &\nBackofen, Bioinfo., 2006; Garcia-Martin et al., J. Bioinfo. Comp. Bio., 2013; Eastman et al., PLOS, 2018; Zhou et al.,\nBioinfo., 2023; Matthies et al., NAR, 2023]. Our work belongs to this category.\n2. That being said, there are prior studies that have experimentally validated the functionality of\ncomputationally designed RNA sequences [Dotu et al., NAR 2014; Yamagami et al., NAR 2018;\nWayment-Steele, Nat. Methods, 2022]. Dotu et al. validated type III hammerhead ribozymes, while\nYamagami et al. validated double-pseudoknotted HDV ribozymes. In both studies, an RNA\ndesign program was used to generate sequences for specific target structures, followed by\nself-cleavage assays on the synthetic RNA molecules (see Reviewer #2, comment 7 for\nexperimental detail). Dotu et al. reported that all 10 selected designs across different selection\ncriteria showed self-cleaving activities with different cleaving rates. Yamagami et al. found that\n11 out of 12 ribozymes with NED < 0.03 showed self-cleaving activity, compared with only 3 out\nof 10 when NED ≥ 0.03. These results demonstrate that RNA design algorithms can successfully\ndesign active ribozymes.\n3. Because all these prior RNA design efforts (with and without experiments), including ours, use\nthe same nearest neighbor energy model, we believe our design will also work in experiments\nsimilar to those cited above in Point #2.\nReviewer #3 (Remarks to the Author):\nTang et al. present the method SamplingDesign, which combines several computational procedures into a\ncompelling RNA design method. Namely, it uses sampling of sequences and their associated structures to\nestimate gradients as well as distributions that explicitly model pairs, mismatches, tri-mismatches, and\nunpaired nucleotides. The combination of these innovations greatly helps to reduce the possible design\nsearch space. This manuscript’s Supplementary Methods are done well and helps with understanding\nSamplingDesign’s algorithm. The authors do a good job benchmarking their new method, which performs\nwell compared to the previous methods SAMFEO, NEMO, and Matthies et al. Benefits of the\nSamplingDesign method include greater flexibility of objective function potentially for different\napplications, ability to design longer RNA structures, and better performance in many metrics on the\nEterna100 dataset.\nMajor points:\n1. Authors should make clear that SamplingDesign is meant for RNA secondary structure design in its\nits\ncurrent form.\nThank you for the suggestion. We now state this explicitly in the Abstract and the opening paragraph of\nthe Introduction:\n“RNA design aims to find an RNA sequence that can fold into a given target secondary structure…”\n2. A little more clarification under “Evaluation setup” #6 and #7. Does #6 mean the number of global\nminimum free energy structures, and #7 mean only counting if there is exactly 1 global minimum free\nenergy structure? It is more precisely stated around Eqn. 2 and 8 in the Methods but would be good to\nhelp the reader more in the main text.\nThank you for the comment. We have revised the subsection “Evaluation setup” in Results to clarify\nmetrics #6 and #7 and added forward references to their corresponding equations (reproduced below).\n3. Fig. 6:\n- caption starts with panel b and not a\n- Panel b: pink coloring is used to denote undesignable region, but the color is not in the caption.\nThank you for pointing this out. We have revised the captions accordingly. For example,\nCaption: “Fig. 6. (a) Comparison of the best Boltzmann probability solutions from SamplingDesign and\nSAMFEO for Puzzle #78. (b) Target structure, with the undesignable motif highlighted in pink. (c): MFE\nstructure…”\n4. Item #2 in Discussion discusses how they test Matthies et al. on GPUs and ran into GPU memory\nissues with longer RNA lengths. The authors do not mention what resources they used for running\nSamplingDesign. Was SamplingDesign also run on an 80GB NVIDIA H100 GPU, and for this type of\nresource can SamplingDesign reasonably jump to designing mRNAs that can be several kb long?\n(a) Re: “What resources they used for running SamplingDesign.” SamplingDesign was run on CPUs\nand not GPUs. All experiments were run on an HPC server with 64 CPU cores (AMD EPYC\n9554 64-Core Processor). Runtime and peak memory for the Eterna100 puzzles are reported in\nFig. S10(a) and (c) (see also Reviewer #2, comment 6). We have added this detail to the Methods\nsection “Implementation details of SamplingDesign”: “We run our program on a\nhigh-performance computing server using 64 CPU cores (AMD EPYC 9554 64-Core Processor).”\n(b) Re: “Can SamplingDesign reasonably jump to designing mRNAs that can be several kb long?”\nRegarding the design of mRNAs, yes, our sampling framework is well-suited for that problem\nbecause it can optimize arbitrary objective functions (e.g., ensemble free energy or average\nunpaired probability); see Discussion, bullet 2, reproduced below:\n“2. We can use this framework to design messenger RNAs (mRNAs) for a given protein,\nwhich has wide applications in mRNA vaccine and mRNA medicine. While mRNA\ndesign optimizing for minimum free energy can be solved efficiently by dynamic\nprogramming (30), optimizing for other metrics such as ensemble free energy (22) and\naverage unpaired probability (31) is much more challenging. This is a perfect fit for\nour sampling framework which can handle arbitrary objective functions.”\nIn fact, before the submission of this manuscript, we have already done preliminary studies on\nmRNA design using sampling, with promising results. For longer mRNAs of several thousand\nnucleotides, we will need to follow our practice of using LinearFold/LinearPartition beam search\nmode (which is much faster than exact search, and scales linearly with sequence length) during\nthe sampling, and only use the exact search mode for the top candidates among the samples in\neach step (see Methods section “Implementation details of SamplingDesign”). We can further\nspeed up using importance sampling.\n5. Item #4 in Discussion: The distribution categories used within SamplingDesign seem designed for\nRNA (pair, unpaired, mismatch, tri-mismatch) so the jump to protein design should be further qualified\nsince these four categories are not directly applicable to protein.\nThank you for pointing this out. We agree that a direct jump from RNA design to protein design is\nnon-trivial, and so we removed this bullet from the Discussion.\n6. How are asymmetric bulges with only paired nts on one side of the helix encoded in the distribution\nversions? Ex. 1 nt bulge. Similarly in Fig.7 right side, the unpaired nts in the 3-way junction are labeled as\nunpaired even though they could be interpreted as a “mismatch”? Please provide more description.\nThank you for the clarification question. In our work, mismatches and trimismatches are only defined for\nhairpins and internal loops.\n(a) In bulge loops, there are no mismatches (because there are unpaired nucleotides on only one of\nthe two sides), and we treat all the unpaired nucleotides in a bulge loop independently. See Fig.\n7(b) in Methods for an example (reproduced below, with the 1-nt bulge circled in red).\nFig. 7(b). Examples of unpaired nucleotides, base pairs, mismatches, and trimismatches in\ncoupled-variable distribution v3.\n(b) For a multiloop (n-way junction), we do not use coupled mismatches either, for simplicity\nreasons, because the mismatch positions might coincide with paired positions in a multiloop. For\nexample, consider the following 3-branch multiloop in (a):\n“Fig. S7. For a multiloop (n-way junction), we do not use coupled mismatches because the\nmismatch positions might coincide with paired positions. (a) The outermost (1,12) pair’s\nmismatch positions are (2,11), which coincide with the (2,6) pair and the (7,11) pair, respectively.\nSimilarly, the (2,6) pair’s mismatch positions are (1,7), which coincide with the (1,12) pair and\nthe (7,11) pair, respectively; and the (7,11) pair’s mismatch positions are (6,12), which coincide\nwith the (2,6) pair and the (1,12) pair, respectively. (b) Examples from Eterna100, where\nmismatch positions that coincide with a paired position are indicated by orange dashed lines.\nThere are 39 puzzles from Eterna100 that show at least one such overlap.”\nTo be clearer, we added the above figure in the Supplementary Figures as Fig. S7 and a clarification\nstatement (see below) in the Methods subsection “Coupled variables for terminal mismatches (v2 and\nv3)”:\n“In this work, we only considered mismatches and trimismatches for hairpins and internal loops.\nThere are no mismatches for bulge loops, and we do not couple unpaired nucleotides in a\nmultiloop or external loop; see Fig. S1 for details.”\n7. The description around the possible design space and quantifying this design space got confusing\nthrough the manuscript. The authors’ major innovation is reducing the search path through their choice in\ndistribution (i.e. base pairs, mismatch, trimismatch, unpaired), but they simultaneously say that the\npossible design space is 4N in Fig. 1 and that 4N is incorrect in other parts of the manuscript (ex. “RNA\ndesign space: full vs. valid”). Please go back through manuscript and make these distinctions clearer.\nPlacing equation 6 earlier in the manuscript may help.\n𝑛\nThank you for pointing this out. To clarify the two different concepts, we now refer to the set of all 4\nsequences as the “total sequence space” and the set of valid sequences 𝑋(𝑦) as the “valid design space,”\n|𝑢𝑛𝑝𝑎𝑖𝑟𝑒𝑑(𝑦)| |𝑝𝑎𝑖𝑟𝑠(𝑦)|\nwhich has size |𝑋(𝑦)| = 4 × 6 . By ruling out invalid sequences, this latter size is much\nsmaller than that of the total sequence space. To make this point clearer in the manuscript, we added a\nnew subsection in Results, “RNA design problem and valid design space,” and moved Eqs. (5) and (6)\nhere. Additionally, we added Fig. 2(b) to illustrate the idea that the valid design space is a small subset of\nthe total sequence space (reproduced below).\n8. Introduction:\n“However, such methods are only capable of exploring a fixed number of candidate sequences, thus not able to keep up with the exponential growth of the design space.”\nHere, the authors could clarify the rationale for their method vs. local search. The wording might imply that\nSamplingDesign can search an exponential number of sequences, whereas local search only explores a\nfixed number. It should be made clear that SamplingDesign’s discrete sampling allows the method to\nescape local minima and explore more of the search space as the number of samples increase.\nSince the objective function is arbitrary, this is an advantage over local search methods that reach a local minimum\nwhen the objective function is non-convex with regards to the search space.\nThank you for the very deep question. We want to first clarify that SamplingDesign works on\ndistributions of all (exponentially many) valid sequences. Even if the sample size is small, in each step we\n|𝑢𝑛𝑝𝑎𝑖𝑟𝑒𝑑(𝑦)| |𝑝𝑎𝑖𝑟𝑠(𝑦)|\nstill choose candidates from an exponentially large space of size |𝑋(𝑦)| = 4 × 6 . By\ncontrast, at each step, local search works on a specific sequence, and can only choose candidates from the\nset of neighboring sequences by changing one unpaired or a pair of base-paired nucleotides. This\nneighbor set has size 3 × |𝑢𝑛𝑝𝑎𝑖𝑟𝑒𝑑(𝑦)| + 5 × |𝑝𝑎𝑖𝑟𝑠(𝑦)| which is linear in |𝑦| instead of exponential\n(although sophisticated methods like SAMFEO can change a bit more nucleotides at a time due to\nstructured mutation, but the size of the neighbor set is still 𝑂(|𝑦|)). This is the fundamental advantage of\nSamplingDesign over local search. In addition, in SamplingDesign, each update changes the whole\ndistribution of exponentially many sequences. We have added the clarification to the Discussion.\nMinor points:\n1. Fig. 4d: ‘*’ is not described in caption or main text.\nEach starred puzzle ‘★’ in Fig. 4(d) links to its corresponding structure diagram. We have updated the\ncaption of Fig. 4(d) (reproduced below) to explain the stars:\n“...Each starred (⋆) puzzle has a secondary structure diagram in Figs. 5–6, S1, S4–S6, and S12.”\n2. In “Softmax parameterization and unconstrained optimization” between equations 31 and 32, there is\nmixed usage of ‘S’ and ‘Z’.\nThank you for pointing this out. We have corrected the mixed usage in the Supplementary Methods.\n3. Is random initialization of SamplingDesign possible? The authors use 3 other types of initialization\nmethods and curious about random initialization because it is used frequently in ML.\nYes, running SamplingDesign with random initialization is possible and should in principle further\nimprove the results. However, doing this (with multiple random initializations) is time-consuming.\nBecause of that, we only used uniform initialization and ε-targeted initialization (favoring GC base pairs\nand unpaired A nucleotides) in this work. We leave the use of random initialization to future work.\n4. Would it be possible within the framework of SamplingDesign to vary sample size between steps?\nEarlier steps may benefit more from larger sample size to estimate the gradient and promote exploration.\nYes, it would be possible to vary the sample size between steps in SamplingDesign, and we agree that\nearlier steps may benefit more from a larger sample size. Our sensitivity analysis shows that the learning\ncurve of the average Boltzmann probability converges faster when we use a larger sample size (in the new\nfigure, Fig. S13(b), reproduced below).\nThat being said, sampling in later stages is also important for exploring promising candidates, and the\ntop-scoring solution often appears in later steps. We can also implement adaptive sample sizes based on\nthe uncertainty (i.e., entropy) of the sequence distribution (a flatter distribution requires more exploration)\nand the amount of shift (i.e., KL-divergence) between two consecutive distributions (using more samples\nwhen the distribution shifts substantially). However, we leave the varying sample sizes to future work.\nOther Changes\n● We updated the results of SamplingDesign on Eterna100 (Fig. 4a) after re-running the evaluation\nwith a fixed random seed for reproducibility. The new results (table shown below; not included in\nthe manuscript) show improvement across most metrics (e.g., the geometric mean of Boltzmann\nprobability and the number of MFE solutions), though the average structural distance is slightly\nhigher.\n● We added a subplot in Fig. S11 (reproduced below) to show the reduction in entropy during\noptimization. The grey solid curve shows the average entropy of sequence distributions across the\nEterna100 dataset.\n⋆\n● We introduced math notations in the main text and replaced the symbol 𝑦 for the target structure\nwith 𝑦. We made these changes to improve the manuscript clarity and address some of the\nconcerns raised (e.g., Reviewer #3, comment 7; Reviewer #3, comment 8).\n● In the main results (Fig. 4a), we added clarification in the caption that each column is evaluated\non different sequences because SamplingDesign and SAMFEO output the best sequence for each\nmetric found during optimization. We added Table S4 (reproduced below), where each column is\nevaluated using the same sequence."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Thanks for the revisions and detailed responses. The provided clarifications and supplementary experiments have adequately addressed my previous concerns and strengthened the manuscript.", "reviewer_2": "The author carefully responded to most of my questions, and I have no further questions.", "reviewer_3": "The authors did a great job in addressing the points raised, and we appreciate their effort. However, a few points remain with the current manuscript:\nMajor:\n1. It is not accurate to consider the “valid design space” |χ(y)| alone as a “novel concept”. However, the authors here leverage the valid design space definition very well in their compelling algorithm and describe the valid design space in an easy-to-understand manner.\nExamples of prior literature that describe the equivalent valid design space concept:\n- Yang, Xiufeng, et al. (2017) BMC Bioinformatics – “Sequence space” in Methods\n- Churkin, Alexander, et al. (2017) Briefings in Bioinformatics – second paragraph in Introduction includes:\n“This upper bound can be refined by noting that paired positions have to form valid base pairs under the standard A-U, C-G, G-U base pairing scheme. This implies that 6^(p/2)*4^u sequences are compatible with a secondary structure having, respectively, u unpaired and p paired nucleotides.”\nExample in manuscript in Results section:\n“…and define a novel concept of “valid design space” X(y)”\nMinor:\n1. Please add y-axis label to Fig. 2c.\n2. This phrase from the introduction is now outdated and replacing “the last two” with “two recent” will be more accurate:\n“the last two Nobel Prizes in Physiology and Medicine for messenger RNA vaccines (2023) and microRNAs (2024).”\n3. The “almost” is an important reflection of the study’s results in this sentence from the introduction. Please remove the parentheses or reword:\n“When tested on the Eterna100 benchmark, our work consistently outperforms the state-of-the-art RNA design methods (4, 5) in (almost) all metrics such as Boltzmann probability and ensemble defect.”\n4. Typo of “puzzles” in second sentence of SI – Learning curves\n5. Could be easier for the reader to remove mixed usage of “continuous optimization” and “continuous relaxation” to describe the same method.\nIn introduction paragraph 4: \"Instead, we borrow two techniques from machine learning, continuous optimization (21, 22) and Monte-Carlo sampling, for the RNA design problem.”\nBut in abstract: “We instead borrow two techniques from machine learning, continuous relaxation and Monte-Carlo sampling, to the RNA design problem”\n6. Should be “ranks second.” under “SamplingDesign outperforms state-of-the-art methods”:\n“On Eterna100, SamplingDesign ranks first on all metrics except average structural distance and the number of uMFE-solution puzzles, where it ranks seconds.”\n7. “learning curves” can be capitalized under “Runtime Analysis” to match referred section title:\n“Other optimization options such as using smaller sample sizes or a larger learning rate can further reduce runtime but also make finding better solutions more challenging (see Supplementary Methods subsection “learning curves”).”", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Thank you for the thoughtful comment. While these prior literatures have used the concept of the valid design space, these descriptions were generally informal. To the best of our knowledge, our work is the first to provide a rigorous definition of this set of all valid designs (Eq. (1)) and apply it directly to continuous optimization. Accordingly, we have revised the paragraph as follows:\n“Therefore, we refer to those 4𝑛 sequences as the “total sequence space” (Fig. 1(a)) and\nformulate the “valid design space” X(y) as\nX(y) = {x ∈ {A, C, G, U}n | ∀(i, j) ∈ pairs(y), xx ∈ P} (1)\ni j\nAlthough this concept has appeared in previous studies (23, 24), this is the first time it has been\nformally defined.”\nMinor:\n1. Please add y-axis label to Fig. 2c.\nThank you for the suggestion. We have added the y-axis label (“valid ratio”) to Fig. 2(c) (reproduced\nbelow).\n2. This phrase from the introduction is now outdated and replacing “the last two” with “two recent” will\nbe more accurate:\n“the last two Nobel Prizes in Physiology and Medicine for messenger RNA vaccines (2023) and\nmicroRNAs (2024).”\nThank you for catching that. We have replaced the phrase in the sentence.\n“... the two recent Nobel Prizes in Physiology and Medicine for messenger RNA vaccines (2023) and\nmicroRNAs (2024).”\n3. The “almost” is an important reflection of the study’s results in this sentence from the introduction.\nPlease remove the parentheses or reword:\n“When tested on the Eterna100 benchmark, our work consistently outperforms the state-of-the-art RNA\ndesign methods (4, 5) in (almost) all metrics such as Boltzmann probability and ensemble defect.”\nThank you for the suggestion. We have revised the sentence.\n“When tested on the Eterna100 benchmark, our work consistently outperforms the state-of-the-art RNA\ndesign methods (4, 5) in key metrics such as…”\n4. Typo of “puzzles” in second sentence of SI – Learning curves\nThank you for pointing this out. The typo has been corrected.\n5. Could be easier for the reader to remove mixed usage of “continuous optimization” and “continuous\nrelaxation” to describe the same method.\nIn introduction paragraph 4: \"Instead, we borrow two techniques from machine learning, continuous\noptimization (21, 22) and Monte-Carlo sampling, for the RNA design problem.”\nBut in abstract: “We instead borrow two techniques from machine learning, continuous relaxation and\nMonte-Carlo sampling, to the RNA design problem”\nThank you for the suggestion. We have updated the abstract to use “continuous optimization”.\n6. Should be “ranks second.” under “SamplingDesign outperforms state-of-the-art methods”:\n“On Eterna100, SamplingDesign ranks first on all metrics except average structural distance and the\nnumber of uMFE-solution puzzles, where it ranks seconds.”\nThank you for catching this. We have corrected “rank seconds” to “rank second.”\n7. “learning curves” can be capitalized under “Runtime Analysis” to match referred section title:\n“Other optimization options such as using smaller sample sizes or a larger learning rate can further reduce\nruntime but also make finding better solutions more challenging (see Supplementary Methods subsection\n“learning curves”).”\nThank you for noting this. We have capitalized “Learning Curves” to match the title."}]

Decision letter