RNA Secondary Structure Prediction Explained

RNA secondary structure prediction is one of the few places in computational biology where a textbook dynamic program is still close to the production method. Understanding why the obvious objective function is wrong is most of what you need.

What secondary structure is

An RNA molecule is single-stranded, so it folds back on itself and pairs with itself. Secondary structure is the list of which positions are paired with which: the set of pairs (i, j) where base i is hydrogen-bonded to base j. It is a much coarser description than the three-dimensional coordinates, and that coarseness is the point — secondary structure is predictable from sequence with a small, well-parameterised model, and tertiary structure largely is not.

The pairs that count are Watson-Crick, A with U and G with C, plus the wobble pair G with U, which is genuinely a pair in RNA and is the detail most often left out of an informal explanation. Two constraints define the search space. Each base is in at most one pair. And in the standard formulation, pairs may not cross: if i pairs with j and k pairs with l, then either the intervals nest or they are disjoint, never interleave. Crossing pairs are called pseudoknots and they are real in nature; excluding them is a modelling choice that makes the problem tractable, and it is the single largest source of missed structure.

The nesting constraint also forbids very tight turns. A hairpin loop needs at least three unpaired bases between the two paired ends, because the backbone cannot bend more sharply than that, so any pair (i, j) requires j minus i to be at least four.

Maximising pairs, and why it fails

The classical starting point maximises the number of pairs. Because pairs nest, the problem decomposes: the best structure on the subsequence from i to j either leaves j unpaired, in which case the answer is the best structure on i to j minus 1, or pairs j with some k in between, in which case the answer splits into the best structure inside the new pair and the best structure to its left. That is a recursion over all sub-intervals, filled in order of increasing interval length:

M[i][j] = max(
    M[i][j-1],                                  // j unpaired
    max over k in [i, j-4] with pair(k, j) of:
        M[i][k-1] + 1 + M[k+1][j-1]             // j pairs with k
)

There are order n squared intervals and each takes order n work to evaluate the inner maximisation, so the algorithm is order n cubed in time and order n squared in memory. That complexity has not changed; every production folding program is still a cubic dynamic program over intervals with the same nesting decomposition.

What has changed is the objective, because maximising pairs gives structures that are wrong in a specific and instructive way. It has no preference for consecutive pairs over scattered ones, so it happily produces a fold full of isolated single pairs separated by loops. Real RNA does the opposite: pairs occur in stacked helices, because the energy that actually stabilises a fold comes largely from the stacking of adjacent base pairs on top of each other, not from the hydrogen bonds of a pair in isolation. An objective that counts pairs cannot express that, because the contribution of a pair in that model does not depend on its neighbours.

The nearest-neighbour energy model

The fix is to score stacked pairs rather than pairs. In the nearest-neighbour model, the free energy of a structure is a sum of contributions from local motifs: each stack of two adjacent base pairs contributes a negative (stabilising) term whose value depends on which two pairs are stacked, and each loop — hairpin, bulge, internal, multibranch — contributes a positive (destabilising) term that grows with the loop’s size. Terminal mismatches, dangling ends and coaxial stacking each get their own terms.

Those terms are measured, not fitted to structures. They come from optical melting experiments on synthesised oligonucleotides, and the parameter set in general use is the Turner 2004 set. The ViennaRNA package documents that its algorithms use the Turner 2004 nearest-neighbour model and ships the table compiled in; the ViennaRNA energy-parameter documentation is where to read the exact values rather than take them from a secondary source. The recursion is then the same shape as above with the pair count replaced by an energy sum and the maximisation replaced by a minimisation.

A worked fold on nine bases

Take GGGAAACCC, positions 1 to 9. The nesting and minimum-loop constraints leave very little room, which is what makes it a good example.

G G G A A A C C C
1 2 3 4 5 6 7 8 9

candidate pairs (G-C only here, j - i >= 4):
  1-7  1-8  1-9  2-7  2-8  2-9  3-7  3-8  3-9

structure A: 1-9, 2-8, 3-7   -> three stacked pairs,
                                hairpin loop = positions 4,5,6 (size 3)
structure B: 3-7 only        -> one pair, hairpin loop size 3
structure C: 1-9, 3-7        -> two pairs, not adjacent:
                                an internal loop between them

Under pair-counting, structure A wins with three pairs and structure C scores two — the right ordering, by luck. Under the energy model the reasoning is different and more robust. Structure A contains two stacks, the 1-9 pair stacked on 2-8 and the 2-8 pair stacked on 3-7, each contributing a negative term, plus one hairpin-loop initiation penalty for the three unpaired A residues. Structure C contains no stacks at all — its two pairs are not adjacent — so it collects two loop penalties and no stabilising stacking term, and comes out substantially worse than A even though it has two pairs to structure B’s one.

That is the mechanism worth keeping: a stack of two pairs is worth more than two pairs, and the energy model is the thing that can say so. Read the actual kilocalorie-per-mole values for the GG/CC stack and for a size-three hairpin off the Turner table rather than reproducing them from an explanation like this one, including this one.

One structure is the wrong answer

The minimum-free-energy structure is a single point estimate, and for many sequences it is barely better than thousands of alternatives. RNA in solution occupies a Boltzmann-weighted ensemble of structures, and the useful output is that distribution rather than its mode.

McCaskill’s partition-function algorithm computes it with the same cubic recursion run with sums and Boltzmann factors instead of a minimisation, giving a probability for every possible pair. Two derived outputs are more informative than the single fold: the base-pair probability matrix, which shows which helices are confident and which are one of several options, and the ensemble free energy, whose gap from the minimum-free-energy value tells you how concentrated the ensemble is. A sequence whose predicted structure is worth trusting is one whose ensemble is concentrated.

Three limits are worth stating plainly. Pseudoknots are excluded by the recursion, and general pseudoknotted folding is computationally hard, so programs that predict them use restricted classes or heuristics. The thermodynamic parameters were measured under specific salt and temperature conditions, and a cell is not those conditions. And RNA folds while it is being transcribed, so the structure a molecule actually adopts can be a kinetic trap rather than the thermodynamic optimum — which is why chemical probing data such as reactivity profiles, folded in as pseudo-energy constraints, improves predictions more than any change to the algorithm does.