New Algorithm Maps Every Optimal RNA Fold Across Branching Parameters

For decades, the workhorse tool for predicting these shapes has been the nearest-neighbor thermodynamic model, or NNTM, which estimates the folding free energy of an RNA structure by breaking it into small components such as base pair stacks, hairpin loops, and multiloops, then summing the contributions of each. The structure with the lowest total free energy, the minimum free energy or MFE structure, is taken as the prediction. Yet a persistent weakness of this framework lies in how it treats multiloops, the substructures where three or more helices branch from a central loop, and a new study published in Algorithms for Molecular Biology presents an algorithmic advance that turns this weakness into an opportunity for exploration.

In the NNTM, the energy assigned to a multiloop consists of an initiation term, interpreted as an entropic penalty, plus a stacking term capturing favorable interactions inside the loop. The initiation penalty is modeled as a simple linear function governed by three branching parameters: a, an initiation penalty; b, a penalty per unpaired nucleotide in the loop; and c, a penalty per branching helix. This linear form was introduced originally for computational efficiency, and later comparisons against penalties based on logarithmic Jacobson-Stockmayer theory and on polymer theory concluded that the simplest model performed best. Unlike many NNTM parameters derived directly from experiments, the historical values of these branching parameters were partly learned by maximizing prediction accuracy over known structures, and they have shifted considerably across model releases, from the triple (4.6, 0.4, 0.1) in 1989 to (3.4, 0, 0.4) in 1999 and then to (9.3, 0, -0.6) in 2004, in units of kcal/mol.

Previous work had shown that reparametrizing the multiloop entropic penalty can yield markedly improved MFE structures, but directly searching for a single better parameter set risks overfitting. The alternative strategy is to treat the branching parameters as variables and collect the ensemble of distinct structures that are optimal under different parameter choices. The obstacle was computational: the earlier algorithm for exploring this space of alternative branching structures did not scale to long RNA sequences, where each MFE calculation takes time that grows cubically with sequence length. The new study, by Svetlana Poznanović, Owen Cardwell, and Christine Heitsch, addresses exactly this bottleneck with an algorithm that is optimal in the number of MFE computations it requires.

The key insight rests on geometry. For a fixed sequence, every possible secondary structure carries a branching signature, a quadruple recording the number of multiloops, the number of unpaired nucleotides within those multiloops, the number of branching helices, and the residual free energy from all other structural components. Because the number of possible structures is finite, the three-dimensional branching parameter space can be partitioned into a finite set of convex regions, each corresponding to one distinct optimal branching signature. Any parameter combination within a region yields the same optimal structure, and boundaries between regions mark parameter values where multiple optimal structures coexist. The new algorithms compute this partition for the plane where the b parameter is fixed, using repeated MFE calculations, each paired with a parse of the resulting structure, a combined operation the authors call querying an oracle.

The algorithm works iteratively. It begins with signatures discovered at the corners of a chosen search frame, then divides the frame into dominating sets, the areas where each discovered signature beats all others. It queries the oracle at the vertices of these sets; if a known signature comes back, the vertex is confirmed as part of the final partition, and if a new signature appears, the partition is updated and new vertices are added for inspection. The authors prove that the procedure finds every signature whose region has non-zero area, and that the total number of oracle calls equals the number of regions discovered plus the number of vertices in the partition. A second, integer-precision version of the algorithm interfaces with the widely used ViennaRNA package, which scales parameters by a factor of 100 and performs its calculations in dekacal/mol; spot-checking against exhaustive enumeration of lattice points found no non-zero regions missed.

To measure performance, the team applied the integer-precision algorithm to the Archive II benchmarking dataset, a collection of 3,948 sequences from 10 RNA families curated by the Mathews Lab at the University of Rochester. The search frame, drawn from earlier work, was a rectangle in transformed coordinates spanning minimum per-multiloop penalties from 300 to 800 and helix penalties from -350 to 100, in dekacal/mol. Because both the number of multiloops and the maximum number of branches scale linearly with sequence length, the number of regions grows quadratically, and empirically each region introduces roughly two new vertices, so total oracle calls also grow as O(n squared) and total runtime as O(n to the fifth), with small coefficients. For sequences up to 200 nucleotides, 60 percent of the dataset, total compute time was well under one second, with a maximum of 0.68 seconds, and even for the nine sequences of 700 nucleotides or more, the geometric overhead beyond the oracle calls itself totaled just 1.95 seconds.

The biological payoff came from comparing the most accurate structure in the branching ensemble, measured by F1 score against known target structures, with standard predictions. For every one of the 10 families, the best structure from the search rectangle was statistically significantly more accurate than the standard MFE prediction, under a Bonferroni correction. Moreover, when compared against both the MFE structure and the Boltzmann centroid, the current standard for alternative structure generation, the best branching structure improved on both for 47 percent of sequences. The improvement exceeded 0.05 in F1 score for 34 percent of sequences and 0.10 for 24 percent, indicating that reparametrization generates a rich ensemble of candidate folds that standard methods never surface.

The analysis also revealed that achievable accuracy depends in non-obvious ways on auxiliary modeling decisions that most users never question. The default predictions used ViennaRNA’s no-lonely-pair option and its d2 dangling end model, in which dangling energies are counted on both sides of a helix. Toggling the lonely-pair option off and switching to the d1 dangling end model, where unpaired bases participate in at most one dangling end, produced additional improvements for 1,592 of the 1,843 sequences where the default best structure already beat both standard predictions. Expanding the search rectangle into the positive unskewed plane, where the minimum per-multiloop penalty is positive, added further gains for 1,005 sequences. In other words, the potential accuracy hiding in the parameter landscape grows when both the model options and the search space are widened.

Detailed case studies of two mid-length sequences, an RNase P RNA from P. alcalifaciens and an SRP RNA from Gibberella zeae, illustrate both the promise and the remaining difficulty. For the RNase P sequence, whose target includes a pseudoknot that NNTM optimization cannot predict, the MFE structure scored 0.49 in F1 accuracy and the centroid 0.58, while the best structure from the branching ensemble reached 0.647, with 13 of 18 alternative structures beating the MFE prediction. Yet some near-miss helices, such as a predicted stem offset from the target by a single nucleotide with different loop structure, are invisible to strict base-pair matching metrics, prompting the authors to suggest structural edit distances or compatibility-adjusted scoring. For the SRP RNA, the best alternative structure scored 0.72 against MFE’s 0.60, but a close relative scored only 0.56, and without knowing the true structure there is no way to tell which candidate is more reliable.

That identification problem, which the authors call finding the needle in the haystack, remains open, alongside challenges including pseudoknots, sparsely paired regions, sparse long-range contacts that likely represent tertiary interactions, and the possible relevance of coaxial stacking and more sophisticated branching energy models. Still, the study establishes that the conventional choice of multiloop parameters may itself limit prediction accuracy, and that exhaustively mapping which structure is optimal for every parameter combination is now tractable even for longer sequences. The authors also note that their half-plane intersection framework applies beyond RNA, to any linear function minimized over a finite set of rational vectors, suggesting the partitioning algorithms could serve wherever parameters shift and optimal solutions must be tracked. Chemical probing data as constraints, Bayesian optimization, and probabilistic scoring are flagged as promising routes toward finally picking the right fold out of the newly accessible ensemble.

Subject of Research: An efficient algorithm for exploring RNA branching conformations under the nearest-neighbor thermodynamic model

Article Title: An efficient algorithm for exploring RNA branching conformations under the nearest-neighbor thermodynamic model

Article References: Poznanović, S., Cardwell, O., & Heitsch, C. (2026). An efficient algorithm for exploring RNA branching conformations under the nearest-neighbor thermodynamic model. Algorithms for Molecular Biology, 21(1), Article 6. https://doi.org/10.1186/s13015-025-00296-4

Image Credits: AI Generated

DOI: 10.1186/s13015-025-00296-4

Keywords: RNA secondary structure, nearest-neighbor thermodynamic model, multiloop parameters, minimum free energy, parameter space partition, ViennaRNA, Archive II dataset, algorithm, prediction accuracy, computational biology, efficient, exploring