Secondary Structure Design for Cotranscriptional 3D RNA Origami Wireframes

Concurrently to the advances in DNA nanotechnology, there has been increasing interest in using RNA as the fabrication material for self-assembling bionanostructures. In comparison to DNA, the appeal of RNA is that the strands can be produced by the natural process of polymerase transcription, and the structures can thus be created in room temperature in vitro, and possibly eventually in vivo, from genetically engineered DNA templates. The challenge, on the other hand, is that the folding process of RNA is kinetically more complex and hence less predictable than DNA helix formation, at least at the present stage of RNA engineering.

The first design technique applied in this area of RNA nanotechnology was modular “RNA tectonics”, in which naturally occurring RNA structures are connected together to create bigger target complexes using specific connector motifs such as sticky-end pairings and a variety of pseudoknots [11, 12]. A complementary top-down method of “RNA origami”, in which a task-specific strand is rationally designed to fold into the desired target structure, was then introduced in a pioneering work by Geary et al. in 2014 [8]. Geary et al. demonstrated the feasibility of their method by synthesising 2D “RNA tiles” of different sizes, and this approach has since then been further developed with new design motifs, techniques, and tools [14, 6]. (For an overview, see [17].)

One emphasis in the work of Geary et al. [8, 6] has been the cotranscriptional nature of the polymerase transcription process, that is, the way the transcribed RNA strand starts to fold into secondary structures already while being spooled off the polymerase enzyme (Figure 1). This characteristic of natural RNA generation introduces new challenges and also opportunities for the rational design process, some of which we shall explore in the present work.

In the following, Section 2 presents a spanning-tree based framework for self-assembly of wireframe structures by co-transcriptional folding, and introduces the topological folding obstacle of polymerase trapping. Section 3 then demonstrates how this obstacle can always be avoided by using a depth-first-search (DFS) based design scheme. Next, Section 4 introduces the notion of the KLX number of a DFS tree, which corresponds to the maximum number of kissing loops that are simultaneously “open” in the folding process, and hence need different designs in order to avoid nonspecific pairings. Minimising this number provides the possibility of efficiently reusing KL designs, although, as proved in Section 5, the KLX minimisation problem is in the general case NP-hard. Since an efficient minimisation algorithm is thus unlikely, Section 6 provides a branch-and-bound type enumeration approach to the problem. Section 7 provides some examples of using the DNAforge tool to compute the DFS tree based designs and KLX minimisation. Section 8 summarises the results and suggests some directions for further work.

An extension of the RNA origami method to the design of 3D wireframe structures was presented by Elonen et al. in [3]. We conduct our discussion in this framework, but the basic ideas apply, mutatis mutandis, also to the task of designing 2D RNA origami tiles (cf. [15]). The general spanning-tree based 3D wireframe design scheme is outlined in Figure 2.

In this scheme, one starts from the targeted wireframe, which in the case of Figure 2(a) is a simple tetrahedron. (Or more precisely the wireframe skeleton of a tetrahedral mesh.) In the first design step (Figure 2(b)) one chooses some spanning tree $T$ of the wireframe graph $G$,¹¹1A spanning tree of a graph $G$ is an acyclic subgraph that connects all the vertices of $G$. and designs the primary structure of the RNA strand so that it folds to create a twice-around-the-tree walk on $T$, covering each edge of $T$ twice in antiparallel directions. In the second design step (Figure 2(c)) one then extends the walk halfway along each of the co-tree (= non-spanning tree) edges of $G$ into a hairpin loop, and designs the base sequences at the termini of the hairpins so that pairwise matching half-edges are connected by the 180° kissing-loop design motif introduced in [8], thus constituting the co-tree edges. Figure 2(d) presents a nucleotide-level model of the eventual nanostructure.

One potentially significant topological obstacle to cotranscriptional folding in this framework is the phenomenon of polymerase trapping, first identified by Geary and Andersen in [9] and also addressed in the recent 2D RNA origami design tool ROAD by Geary et al. [6, pp. 551–552, SI p. 102 ff]. Our mathematical model of this phenomenon is basically the same as introduced and analysed by Mohammed et al. in [15], but adapted to the present setting of 3D wireframe origami designs.

To explain this concern, let us review the previous tetrahedron design, presented in more detail in Figure 3. Figure 3(a) shows the tetrahedral wireframe as a Schlegel diagram, that is, as a planar projection from a point above one of the tetrahedron’s faces. The edges of the chosen spanning tree, which in this case is a 3-pointed star, are indicated by solid black lines, and the co-tree edges by dashed green lines.

Figure 3(b) depicts again the corresponding twice-around-the-tree strand routing (blue) and the complementary kissing-loop pairings (green and red). The helix junctions in the design, which constitute the vertices of the eventual 3D nanostructure, are now indexed according to their first-visit order in the strand routing.

The schematic in Figure 3(c) presents the design as a domain-level arc diagram, where the strand is laid out along a line in the 5´ to 3´ direction, the vertex visits are marked by the corresponding indices, the domain-to-domain helical pairings are indicated by solid blue arcs, and the kissing-loop pairings by green and red dashed arcs. (For simplicity and clarity, the pairings of the half-edge stem domains flanking each kissing-loop hairpin are not shown.)

Consider now how a cotranscriptional folding process of this structure might proceed. Instead of thinking of the RNA strand being spooled out of the polymerase starting at the 5´ end and folding as the appropriate base pairings become available, it may be easier to visualise the large polymerase–DNA-template complex as traversing the 5´-3´ strand route outlined in Figure 3(c) and transcribing the nucleotide domains as it goes.

First the domains 1-2 and 2-3 are transcribed, and the RNA strand stays linear until the transcription of domain 3-2 begins. (For simplicity, we are ignoring any transient nonspecific nucleotide pairings that arise during the folding process.) Between the completion of domain 2-3 and the initiation of domain 3-2, the two opening hairpins for the kissing loops 3-4 and 3-1 are transcribed. (The best relative ordering of these two transcriptions is a geometric and sequence-design issue, and we leave this choice open in this high-level view.) After (or while) the complementary domains 2-3 and 3-2 hybridise, domain 2-4 is transcribed, and after that the closing hairpin of the 3-4 kissing loop and the opening hairpin of the 4-1 kissing loop, in some order.

Consider now what happens when the polymerase reaches domain 4-2 (marked with a red cross in diagram 3(b)), where it should create a double-stranded helix with domain 2-4, by winding the strand around it in antiparallel direction. If the 3-4 kissing loop (drawn in red and marked with a red arrow in 3(b)) has already closed, the strand with the big polymerase–DNA complex coupled to it cannot achieve this, since the kissing-loop pairing is blocking the pathway. (This is of course also a time-scale issue, and depends among other things on the strand distance between the closing hairpin of the kissing loop and the closing domain of the helical pairing; but let us again ignore these lower-level details at this presentation.²²2Furthermore, moving from our strand-centric to a polymerase-centric view, the polymerase of course does not stop transcribing even if the folding is temporary blocked. Eventually the polymerase detaches from the fully transcribed strand and leaves it to fold the unfolded 3´ tail the best it can. From this perspective, the polymerase trapping phenomenon is more of a kinetic than a topological obstacle.)

Schematically, one can see that the risk of this kind of “polymerase trapping” obstacle emerges when a kissing-loop pair, initiated (opened) before a given helical pairing, closes between the opening and closing of the helical pairing; or in terms of our arc diagram when a “dashed arc” that has been initiated before a “solid arc” terminates inside that solid arc.

As another example, let us consider the tetrahedron design presented in Figure 4. As shown in Figure 4(a), in this case the spanning tree is a simple 4-vertex path. Figure 4(b) again outlines the corresponding strand route and kissing-loop pair arrangement, with the helix junctions numbered according to their first-visit order. As witnessed by Figure 4(c), this time there is no risk for the polymerase trapping obstacle. That is, every kissing loop closes only after the completion of all the helical pairings that have been initiated after the kissing loop was opened.

Such complete absence of polymerase traps seems like a very particular property, and one may wonder for which kinds of wireframe models this situation can be achieved. As we shall see in the next section, however, such an arrangement of the helical and kissing loop pairings can in fact be found for any connected wireframe graph, by an application of the fundamental algorithmic method of depth-first search [1, Sec. 20.3].

To streamline the presentation, we assume henceforth that any graph under consideration is connected and undirected, unless stated otherwise. The depth-first search (DFS) method for systematically traversing and labelling a (connected, undirected) graph is presented as Algorithm 1.

Consider a graph $G=(V,E)$ with $|V|=n$ vertices and $|E|=m$ edges. Then a DFS traversal of $G$, starting from any chosen root vertex $r\in V$, partitions the set of edges $E$ in time $O(m)$ in two disjoint classes: $n-1$ tree edges and $m-n+1$ back edges. The tree edges constitute a spanning tree $T$ of $G$, which can be considered to be rooted at $r$ and oriented accordingly, whereas the back edges (which constitute the corresponding co-tree $E\setminus T$) have the important property that they can always be oriented to point “upward” towards the root of the tree [13]; in other words, there are no “cross edges” connecting two different branches of the directed tree, as exemplified in Figures 5 and 6 (c), and also in Figure 7.³³3The simple reason for the absence of cross edges is that if for an edge $e=\{u,v\}$, vertex $v$ is still unmarked when the traversal first considers $e$ at vertex $u$ (or vice versa), then $e$ becomes a tree edge.

To make this precise and to introduce another important notion, consider such a DFS (spanning) tree for a graph $G=(V,E)$ to be a four-tuple $T=(V,S,r,\delta )$, where $S\subseteq E$ is the set of tree edges, $r\in V$ is the chosen root which determines the orientation of the tree, and $\delta :V\to [1..n]$ is a pre-ordering that labels the vertices in order of their first visits in the traversal. As discussed above, any edge of $G$ is either part of the stem $S$ or connects a vertex and its ancestor along the unique path in $T$ from the vertex to the root. We often consider $T$ as embedded in the plane so that the children of each vertex are ordered from left to right in increasing order of their $\delta$ values. Thus, for example, vertex $d$ in Figure 6(c) is presumed to have been visited earlier than vertex $h$ in the traversal that created the tree.

A plane embedded DFS tree $T=(V,S,r,\delta )$ provides an initial blueprint for designing a cotranscriptionally folding RNA wireframe nanostructure. In the first stage of the design process, the contour of the tree $T$ is traced by an RNA strand $w$ so that each edge of $T$ becomes assembled as an RNA helix, hybridised from two complementary antiparallel domains of $w$. In the second stage the co-tree edges of $G$, which are now back edges of the DFS tree $T$, are constituted as kissing loops made of two half-edge hairpins that extend to meet from the end-vertices of the respective edges.

Let us call the sequence of vertex labels encountered during a walk around the contour of any labelled, plane embedded tree $T$ a contour trace or linear arrangement of $T$.⁴⁴4Also known as a “twice-around-the-tree walk” and full walk in [1, p. 1112]. In the special case that $T$ is a DFS tree we refer to this sequence as a DFA arrangement based on $T$. By construction, any contour trace of a tree $T$ accommodates every edge $e$ of $T$ exactly twice: firstly, when the traversal crosses $e$ in a forward direction, and secondly when the traversal crosses $e$ in the opposite direction. Since any DFS tree (or more generally any spanning tree) of a graph $G$ with $n$ vertices has $n-1$ edges, any DFS arrangement in $G$ contains $2n-1$ vertex labels. For example, every DFS tree of a tetrahedron is a simple 4-path in its Schlegel diagram, and the corresponding DFS arrangement has the pattern $1,2,3,4,3,2,1$ (see Figures 4(a) and (c)).

To complete our blueprint for designing a cotranscriptionally folding RNA nanostructure, the DFS arrangement needs to be augmented into a (domain-level) arc diagram by adding information that indicates when, for each pair of hybridising domains $(\alpha ,{\alpha }^{\ast })$, the forward domain $\alpha$ is transcribed and when its complementary reverse domain ${\alpha }^{\ast }$. For each such pair, these time points are connected by an edge (“arc”), with the first time point considered as the “opening” and the second the “closing” time for this pair. Since in our design scheme, DFS tree edges correspond to helical domains, we correspondingly say that a tree edge is opened when it is crossed by the contour traversal in the forward direction, and closed when it is crossed in the reverse direction. Note that the chosen DFS completely determines the scheduling of the opening and closing pattern of the helical domains.

Let us then consider scheduling the opening and closing times of the back edges, i.e. kissing loops in the eventual RNA nanostructure. Multiple occurrences of internal vertices (neither a root nor a leaf) provide us with freedom to choose, for each back edge $(a,b)$, at which visits of $a$ and $b$ each of the two constituent hairpins of it is to be formed. For each back edge $(a,b)$, an occurrence of $a$ and that of $b$ can be chosen and connected by an arc; then we say that this edge is opened at the time point which corresponds to the left end of the arc, and closed at the right end.

The freedom in timing opening and closing of back edges brings multiple arc diagrams based on the same DFS arrangement. Which should we choose then? Given an arc diagram of $G$, we say that a cycle in $G$ is closed at an edge $e$ if $e$ is in the cycle, and in the arc diagram, all the other edges involved in the cycle have already been closed before $e$ is closed. It is known to be kinetically unfavourable to close a cycle at a tree edge. This experimental obstacle motivates our “phloem principle” for ordering back edge connections.⁵⁵5Phloem is a pathway for transporting products of photosynthesis from leaves to the rest of a plant. Recall that the vertices $a$ and $b$ are labelled with distinct integers $\delta (a)$ and $\delta (b)$, respectively, and whichever labelled smaller is an ancestor of the other. This principle says that, if , then the assembly of this edge (by a kissing loop) should begin at $b$ and end at $a$.

A set of kissing loop types should be as orthogonal as possible to minimise the risk of crosstalk, that is the risk of mismatching hairpins hybridising. However, sets of kissing loops that have proven orthogonal enough in the laboratory are limited in size (see, e.g., [10]). The size of largest orthogonal KL sets available in reality serves as a standard for deciding whether a specific (rooted, ordered) DFS tree should be chosen or not. If the design leaves more than this number of kissing loops open simultaneously at any point of folding, it increases the risk of crosstalk.

Recall that any DFS arrangement corresponds one-to-one with a rooted and ordered DFS tree. Given an arc diagram $D=((p_1,\mathrm{\dots },p_m),\alpha )$ of a DFS tree $T=(V,S,r,\delta )$ of a graph $G=(V,E)$, the KLX number of a segment $(p_i,p_{i+1})$ is the number of arcs that cross the vertical line drawn between $p_i$ and $p_{i+1}$. It is defined formally as

The maximum of these values across all segments is the KLX number of this diagram $D$, that is, . Finally, the KLX number of the graph $G$, denoted by $\kappa (G)$, is the minimum among the KLX numbers of all the possible arc diagrams of $G$.

The 1-to-1 correspondence between DFS arrangements and pairs of a graph and its rooted and preordered DFS tree justifies the introduction of the notation $\kappa (G,T)$ as an alias of $\kappa (G)$.

The minimum KLX (kissing loop crossing) number problem (MinKLX) asks, given an undirected, connected graph $G$ and a positive integer $k$, if there exists an arc diagram of $G$ whose KLX number is at most $k$. This problem is computationally hard as stated in the next theorem; its proof can be found in the technical appendix A.

Another relevant problem is that of determining the depth of the shallowest DFS tree for a given graph $G$, because finding a shallower DFS tree decreases the number of helical domains kept open in parallel, and may result in a sequence with fewer helical domain types. This quantity is closely related to the tree-depth of $G$, which is defined to be the depth of the shallowest DFS tree for a supergraph of $G$. It is NP-hard to compute the tree-depth even for the class of triangulated graphs [2].

Considering the NP-hardness of the MinKLX problem, the only approach to finding a KLX-minimal, cotranscription-friendly designs for a given graph $G$ may be via full enumeration of all the DFS trees for $G$. In the next section, we present a branch-and-bound approach to this enumeration process.⁶⁶6As an aside, if one gives up the goal of cotranscription-friendly design, a strand routing that minimises the total number of kissing loops needed to complement it can be found efficiently for any graph, and this number is for many typical wireframe models just zero or one [4].

Before proceeding, let us note that Theorem 5 does not exclude the possibility that KLX-optimal DFS trees could still be found in polynomial time for some graph classes of practical importance such as 3-regular or polyhedral graphs. As any graph can be approximated by a 3-regular one by replacing each vertex by a network of vertices of degree 3, it would be quite interesting to know whether the KLX-minimisation task remains hard for this class. As regards the class of polyhedral graphs, it is known that the MinKLX-related cutwidth minimisation problem (described in the proof of Theorem 5 in Appendix A) is NP-hard for the class of planar graphs with maximum degree 3 [16].

Let us now propose a prototype of algorithmic enumeration of DFS trees for computing the KLX number of a given graph $G=(V,E)$. Algorithm 2 is its pseudocode. Starting from the empty forest, it enumerates the DFS trees by inserting edges of $G$ one by one in a predetermined order as a tree edge, that is, as an edge of a DFS tree to be built up. Biconnectedness is utilized, the property of a graph being free from an articulation point, whose removal disconnects the graph. DFS trees of a graph cannot circumvent any biconnected components, or blocks, of $G$. We shall demonstrate how each local branching of a DFS tree orients edges of $G$ globally once it is accommodated inside a block of $G$. Orientations thus insisted by more than one such local branching are highly likely to contradict each other, enabling the incremental enumeration of DFS trees to prune the edge-insertion-based search tree. As a shallower tree involves more branches, the tighter an optional upper bound on the depth of DFS trees to be obtained is set, the more effective this block-based pruning should get.

A graph $G=(V,E)$ and its spanning tree $T=(V,S)$ can be uniquely decomposed into a set of biconnected components, or blocks, $E_i\subseteq E$, along with a spanning tree $T_i=(V_i,S_i)$ that is the intersection of $T$ and $G_i=(V_i,E_i)$, the induced subgraph of $G$ by $E_i$; unless confusion arises, we may call even $G_i$ a block. The blocks can then be uniquely organised into a so-called block-cut tree by connecting two blocks by an edge if these blocks induce the subgraphs of $G$ that share (exactly one) vertex in common, which is an articulation point of $G$ (see Figure 7 for an example); any articulation point of $G$ thus serves as an interface among multiple blocks. We say that $G$ is the underlying graph of this block-cut tree. Each block $G_i$ is one of the following types:

• $\mathrm{■}$

  Branching. if $T$ involves a vertex that is incident to at least three edges in $E_i$ (branch);

• $\mathrm{■}$

  Spinal.   otherwise.

Every block contains at least one edge of $T$ and no blocks share an edge; hence, the block-cut tree is composed of at most $|S|=|V|-1$ blocks. If a node in a block $E_i$ is incident to at least two tree edges of the block, say $e_1,e_2\in S\cap E_i$, then we say the node is internal. Hence, a spinal block can contain at most two non-internal nodes. Note that an internal node in a spinal block can be a branch of $T$; unlike a branch in a branching block, this branch is an articulation point of $G$. Therefore, a spinal block can have three or more interfaces. A spinal block consisting of a single edge is particularly called a bridge in [18], and we will borrow this term when it matters whether a spinal block involves an internal vertex or not.

The cotranscription-friendly DFS-tree based design method presented in Section 3 is implemented and available for use in the online design tool DNAforge (https://dnaforge.org), together with an option for minimising the KLX cost of the design with a preliminary version of the enumeration method presented in Section 6.⁷⁷7Design method ST-RNA, additional parameters “co-transcriptional route” and “minimise the number of kissing loop sequences”.

Figures 11 and 12 illustrate some outcomes from the tool. Figure 11 shows designs of wireframe dodecahedra based on a randomly chosen spanning tree (upper row) and a DFS spanning tree (lower row). The DFS-tree based design has also been KLX-optimised, resulting in a reduction from a KLX number of 9 in the initial DFS tree to 6 in the optimal one. (The spanning tree diagrams in Figures 11(b) and (e)) have been manually reconstructed from the tool-generated diagrams in Figures 11(c) and (f).) Figure 12 displays random-tree and DFS-tree designs for a 66-vertex, 192-edge wireframe model of a bunny. Also here the DFS-tree based design has been KLX-optimised, resulting in a KLX number reduction from 60 in the initial DFS tree to 33 in the optimal one. Table 1 summarises the KLX number reductions for some basic mesh models.

We have presented models and algorithms for addressing two tasks in secondary structure design for cotranscriptionally folding DNA origami wireframe nanostructures: avoiding the topological folding obstacle of polymerase trapping and minimising the number of distinct kissing loop designs (the KLX number). The key tools in this work have been the algorithmic method of depth-first search in graphs and the ensuing DFS spanning trees. Our branch-and-bound approach to the KLX minimisation problem can also be used for any other effectively computable objective function on DFS trees, such as the DFS tree depth of a given graph (the TD number).

Relevant directions for further work include for instance the following:

1. 1.

   Nucleotide-level sequence design for DNA origami wireframes in the cotranscriptional setting.

2. 2.

   Efficient combinatorial algorithms for minimising the KLX and TD numbers in some interesting classes of graphs, such 3-regular or polyhedral graphs, or proving the problems NP-hard in these classes.

3. 3.

   Efficient fixed-parameter or approximation algorithms for minimising the KLX and TD numbers in some relevant classes of graphs.