An Efficient Algorithm to Compute the Minimum Free Energy of Interacting Nucleic Acid Strands

Decades ago, the deep relationship between secondary structures and dynamic programming algorithms was established [47, 46, 25, 24, 39, 23]. If a secondary structure can be drawn as a polymer graph without edge crossings it is called unpseudoknotted (Figure 1(c)). The earliest polynomial time algorithms were for single-stranded unpseudoknotted secondary structures, where the absence of crossings allows for planar decompositions of secondary structures suited to dynamic programming techniques.²²2We do not study pseudoknotted structures here, indeed most literature ignores pseudoknots for both modelling and algorithmic considerations: Energy models for pseudoknots are difficult to formulate for geometric reasons [10]. Pseudoknotted MFE prediction is NP-hard even for a single strand [1, 20, 21]. Early NP-hardness results [1, 20] used a simple energy model (no loops, only consecutive base pairs forming “stacks” contribute to free energy). Already seeing hardness results with simple energy models make it unlikely that MFE prediction problem for more realistic energy models [9] is tractable. Nevertheless, dynamic programming algorithms exist for restricted classes of pseudoknots, for both MFE [28, 37, 7, 18, 27] and partition function [11, 12]. For a single RNA/DNA strand, both MFE and partition function are computable in $𝒪(N^3)$ time5 (Table 1), using the standard nearest neighbour energy model³³3The standard model we use is variably called the nearest neighbour model, the Turner model, or loop energy model. Versions of the model have been implemented in software suites such as NUPACK [10, 12, 15], ViennaRNA [19] and mfold [45], for both RNA and DNA [29, 30]. that will be formally defined in Section 2.

Work in DNA computing [42, 26, 35, 38, 14, 6, 32, 44], and nucleic acid nanotechnology more generally [16, 41], involves building molecular systems and structures with, to date, hundreds, and soon, thousands, of interacting strands, so there is a need for better algorithms for these multi-stranded “inverse design problems” [8, 13]. And, of course, biologists need to understand molecular structure in order to understand and predict molecular interactions. However, when there are multiple interacting strands, the situation becomes significantly more complicated than the single-stranded case for two reasons: First, for a secondary structure to be unpseudoknotted, it implies there should be at least one permutation of the strands without crossings on the polymer graph [10] (Figure 1). Second, if strand types are repeated then so-called rotational symmetries (Figure 2) arise that need to be accounted for in the model to match the underlying statistical mechanics⁴⁴4This fact from statistical mechanics is discussed in some papers [10, 17], although we’ve not found its full derivation in the modern nucleic-acid algorithmic literature. We leave a first-principles derivation for future work., otherwise structures may be over- or undercounted, leading to incorrect probabilities in the Boltzmann distribution, in other words: incorrect predicted free energy of a secondary structure.

For multiple strands, albeit a constant number $c=𝒪(1)$, Dirks, Bois, Schaeffer, Winfree and Pierce [10] gave a polynomial time partition function algorithm that runs in time $𝒪(N^3(c-1)!)$.⁵⁵5We note that in the original literature [23, 10] that established these results, the polynomial is to the power 4 due to interior loop contributions (i.e. $𝒪(N^4)$ for single-stranded and $𝒪(N^4(c-1)!)$) for multistranded). The subsequent improvement in the exponent from 4 to 3 is achieved using techniques in Dirks and Pierce (2003) [11] to handle the energy contributions of interior loops more efficiently in dynamic programming. Interestingly, Boehmer, Berkemer, Will, and Ponty [3] recently reduced the $c$-strand parameterized running time for computing partition function and symmetry-naive MFE (i.e. ignoring rotational symmetry) from factorial, $𝒪(N^3(c-1)!)$, to exponential, $𝒪(N^3{3}^c)$. In future work, it could be possible to similarly improve the $(c-1)!$ factor in our MFE algorithm via this result.

In the case where the number of strands is non-constant, in particular $c=𝒪(N)$, Codon, Hajiaghayi and Thachuk [9] showed MFE is NP-hard, and even APX-hard.⁶⁶6This hardness result holds whether or not rotational symmetries are accounted for in the energy model. For the second, rotational symmetry, problem, in order to compute partition function, Dirks et al. [10] found an algebraic link between the overcounting and the rotational symmetry correction problems, which allowed both to be solved simultaneously, aided by the exponential nature of the partition function. Surprisingly, that trick does not work for MFE: Since MFE prediction is minimization-based, there is no secondary structure overcounting problem in MFE prediction. In other words, unlike the case of summation-based partition function, algorithmic examination of repeated secondary structures will not change the outcome of energy minimization. Hence, the absence of the overcounting problem makes MFE prediction harder to solve, and was left open by Dirks et al. [10]. For the special case of $c=2$ strands, Hofacker, Reidys, and Stadler [17] gave an $𝒪(N^6)$ algorithm.

We give an efficient solution to the $𝒪(1)$-strand MFE problem, the first that runs in polynomial time. Our main result is stated as the following theorem, whose proof is in Section 5:

In Section 5 we give a time-space trade-off for our result, by showing a variation of the algorithm runs in $𝒪((N^4\log N)(c-1)!)$ time but $𝒪(N^3)$ space.

We use the standard [10] definition of free energy (Equation 2) of multistranded unpseudoknotted secondary structures, which includes rotational symmetry, see Section 2 for formal definitions. We first give an extensive overview of the proof and paper structure, followed by future work.

Our algorithm runs in polynomial time $𝒪(N^4(c-1)!)$ for the case of $c=𝒪(1)$ strands, the $(c-1)!$ term coming from the fact that our algorithm, as well as the snMFE algorithm from Dirks et al. [10], is assumed to be called from an outer loop that explicitly tries all $(c-1)!$ cyclic strand permutations. Can we increase the number of strands and still have a polynomial time algorithm? We know “not by much”, since the problem is NP-complete when $c=𝒪(N)$ in a simpler energy model that merely maximises the count of base pairs [9].

Our MFE algorithm exploits a polynomial upper bound, $𝒰$, on the number of so-called symmetric secondary structures, or distinct pizza cuts. That bound is linear in “most” cases (Lemma 26), but quadratic in one special subcase (Lemma 24) of 2-fold rotational symmetry with a central internal loop. Reducing that special case to linear would subtract one from our algorithm’s running time and space exponent.

The computational complexity of partition function for DNA/RNA strands is less well understood than MFE. Table 1 shows an open problem for partition function on multiple strands. Intuitively, it seems that partition function should be at least as hard as MFE, however that intuition is tempered by the fact that Dirks et al.’s approach found an efficient algorithm for partition function but not MFE. Nevertheless: are there settings where partition function, or problems counting numbers of structures, are #P-hard?

The graph representation of a secondary structure $S$ is the graph $G=(V,E)$, where $V$ is the set of bases of each strand $s\in \{1,\mathrm{\dots },c\}$, and $E=E_v\cup E_b$, where $E_v$ is the set of covalent backbone bonds connecting base $i_n$ with base ${(i+1)}_n$ for all bases $i=1,2,\mathrm{\dots },|n|-1$ on all strands $n\in \{1,\mathrm{\dots },c\}$, and $E_b=S$ is the set of base pairs in $S$. $E_v$ and $E_b$ are disjoint.

The set of circular permutations, $\mathrm{\Pi }$, of $c$ strands has size $(c-1)!$ [5] (e.g., for the three strands $\{A,B,C\}$, $\mathrm{\Pi }=\{ABC,ACB\}$), e.g., the orderings $ABC$, $BCA$, and $CAB$ are the same on a circle. Next, we define a polymer graph for each $\pi$, see also Figure 1.

In the nearest neighbour energy model [22, 30, 10], any connected unpseudoknotted secondary structure $S$ can be decomposed into different loop types [36, 30, 22]: namely hairpin loops, interior loops, exterior loops, stacks, bulges, and multiloops as shown in Figure 1. As usual, let $k_B$ be Boltzmann’s constant and $T$ is the temperature in Kelvin (also a constant).⁸⁸8All results hold if we assume these are typical values from physics, or just 1 in appropriate units. The free energy of $S$ is defined as the sum of three terms:⁹⁹9Throughout this paper $\log n={\log }_{\mathrm{e}}n$.

• $\mathrm{■}$

  the first is itself the sum of the (well-defined, empirically-obtained) free energies $\mathrm{\Delta }G(l)$ of $S$’s constituent loops [10], where the energy of each loop $l\in S$, is defined with respect to the free energy of the unpaired reference state.

• $\mathrm{■}$

  $\mathrm{\Delta }{G}^{\text{assoc}}$ is the entropic association [10] penalty applied for each of the $c-1$ strands added to the first strand to form a complex of $c$ strands.

• $\mathrm{■}$

  $R$ is the rotational symmetry of the secondary structure $S$, illustrated in Figure 2, and to be formally defined in Section 2.3. In particular, since favourable free energies are usually negative, the term $k_{\mathrm{B}}T\log R\ge 0$ corresponds to the reduction in the entropic contribution of $S$, as any secondary structure with an $R$-fold rotational symmetry has a corresponding $R$-fold reduction in its distinguishable conformational space [10] as shown in Figure 2.¹⁰¹⁰10This is perhaps counter-intuitive. In a follow-up expanded version of this paper we will give a full statistical mechanics explanation, which uses the symmetry penalty to offset the fact that non-symmetrical structures are overcounted from algorithmic perspective. Dynamic programming algorithms to date for multi-stranded MFE ignored this term for reasons we outline in Section 2.3.

For $c$ strands, we let $\mathrm{\Omega }$ be the set (usually called the ensemble) of all connected unpseudoknotted secondary structures. For any circular permutation $\pi \in \mathrm{\Pi }$ of the $c$ strands, let $\mathrm{\Omega }(\pi )\subseteq \mathrm{\Omega }$ be the subset of $\mathrm{\Omega }$ such that each connected unpseudoknotted secondary structure $S\in \mathrm{\Omega }(\pi )$ is representable as a crossing-free polymer graph with circular permutation $\pi$.

Predicting the minimum free energy means finding a minimum over the ensemble $\mathrm{\Omega }$. The known strategy is to deal with each partition $\mathrm{\Omega }(\pi )$ separately, then finding their minimum:

Here, we formalise rotational symmetry. In the previous section we assigned each one of the $c$ strands a unique identifier, dealing with them as distinct strands even if two or more have the same sequence. But in most experimental settings, strands with the same sequences are indistinguishable in the sense that they behave identically with respect to relevant measurable quantities [10]. Mathematically, we say that two strands are indistinguishable if they have the same sequence. Also, two secondary structures are indistinguishable if there exists a permutation of the implied unique strand ordering (Remark 6), that maps indistinguishable strands onto each other while preserving all base pairs, otherwise, the two structures are distinct [10].

For any $c$ strands, not necessarily distinct, they consist of $k\le c$ strand types, usually denoted by uppercase English letters $X,Y,\mathrm{\dots }$¹¹¹¹11In contrast with some of the literature. we exclude using $\{\mathrm{A},\mathrm{T},\mathrm{G},\mathrm{C}\}$ for strand types; to avoid any confusion between strand and base types. A multi-stranded DNA system $M=\{(t_1,n_1),(t_2,n_2),..,(t_k,n_k)\}$, is a multiset of $k$ strand types $t_1,\mathrm{\dots },t_k$ with repetition numbers $n_1,\mathrm{\dots },n_k\in \mathbb{N}$ such that $n_1+\mathrm{\dots }+n_k=c$.¹²¹²12It is known [31] how to efficiently reduce the circular permutation space by getting rid of circular permutations that are redundant due to indistinguishable strands, which is important to consider when computing the partition function but needed for MFE. For such a multiset $M$ we can think of each circular permutation $\pi$ as a string over strand types such that each strand type $t_i$ appears exactly $n_i$ times (e.g., $M=\{(X,6),(Z,3)\}$ one valid $\pi$ is $XZXXZXXZX$).

For example, $\{1,2,4\}$ are the symmetry degrees of $\pi =XZXZXZXZ$ which follows from $\pi ={(XZXZXZXZ)}^1={(XZXZ)}^2={(XZ)}^4$.

For any circular permutation $\pi$, its maximum symmetry degree is denoted $v(\pi )$, and the corresponding repeating prefix $x$, such that $x^{v(\pi )}=\pi$, is the fundamental component of $\pi$. It can be seen that $x$ is the smallest prefix that repeats over $\pi$. Indeed, $v(\pi )$ is the number of cyclic permutations that map each strand to a strand of the same type. Any repeating prefix of $\pi$ must be a multiple of its fundamental component, see Lemma 33, Appendix A in the full version [33].

We can visualize any ordering $\pi$ by representing it as a regular $v(\pi )$-gon with each of its $v(\pi )$ vertices representing a fundamental component. Let $\rho =(\mathrm{1\; 2\; 3}\mathrm{\dots }v(\pi ))$,¹³¹³13Here we use algebraic cycle notation. The order of $\rho$, denoted by $o(\rho )$, is the length of $\rho$ which is $v(\pi )$. and consider the cyclic group¹⁴¹⁴14$G^{\pi }$ is isomorphic to $C_{v(\pi )}$, cyclic group of order $v(\pi )$. $G^{\pi }$ generated by $\rho$. Intuitively, $G^{\pi }$ is the group of all $v(\pi )$ rotational motions in plane of the regular $v(\pi )$-gon that give the same $v(\pi )$-gon. We can represent $G^{\pi }$ as follows: $G^{\pi }=\{{\rho }^0,{\rho }^1\mathrm{\dots },{\rho }^{v(\pi )-1}\}$, where ${\rho }^i$ represents rotation of the regular $v(\pi )$-gon by the angle of $i\times \frac{{360}^{\circ }}{v(\pi )}$, where $|G^{\pi }|=v(\pi )$.

Now, we are ready to define the rotational symmetry of a secondary structure and a strand ordering $\pi$, intuitively, the number of rotations of its polymer graph that give the same polymer graph, as shown in Figure 2.

The following lemma restricts us to deal with only specific and constant number of different folding rotational symmetries, and hence constant different symmetry corrections in total. The proof of this lemma is in Appendix A in the full version [33].

We want to slice any $R$-fold rotationally symmetric secondary structure $S$, like pizza, to the centre of its $\mathrm{Poly}(S,\pi )$, without intersecting any of its base pairs. First, we formalise (Definition 15) a special type of backbone cut, called an admissible $R$-symmetric backbone cut. Then, we will prove (Lemma 16) its existence for $S$.