Characterizing the bending and flexibility induced by bulges in DNA duplexes

Advances in DNA nanotechnology have stimulated the search for simple motifs that can be used to control the properties of DNA nanostructures. One such motif, which has been used extensively in structures such as polyhedral cages, two-dimensional arrays, and ribbons, is a bulged duplex, that is two helical segments that connect at a bulge loop. We use a coarse-grained model of DNA to characterize such bulged duplexes. We find that this motif can adopt structures belonging to two main classes: one where the stacking of the helices at the center of the system is preserved, the geometry is roughly straight and the bulge is on one side of the duplex, and the other where the stacking at the center is broken, thus allowing this junction to act as a hinge and increasing flexibility. Small loops favor states where stacking at the center of the duplex is preserved, with loop bases either flipped out or incorporated into the duplex. Duplexes with longer loops show more of a tendency to unstack at the bulge and adopt an open structure. The unstacking probability, however, is highest for loops of intermediate lengths, when the rigidity of single-stranded DNA is significant and the loop resists compression. The properties of this basic structural motif clearly correlate with the structural behavior of certain nano-scale objects, where the enhanced flexibility associated with larger bulges has been used to tune the self-assembly product as well as the detailed geometry of the resulting nanostructures.


I. INTRODUCTION
DNA is one of the most important molecules in biology. By virtue of the specificity of Watson-Crick base-pairing, 1 DNA is highly programmable, and beginning with the work of Nadrian Seeman in the 1980s, 2 it has been identified as a major player in nanotechnology, with many structures and devices already realized. 3,4 For example, by starting with appropriate oligonucleotides, DNA can be made to self-assemble into many structures with high yields by cooling solutions from high temperatures. The target structure is usually designed to be the global free energy minimum by virtue of containing the largest number of base pairs. In many examples, the structures are made up of a few basic components, including double helical sections and various types of junctions where the double helices meet. Examples of these nanostructures include DNA polyhedra, such as tetrahedra, 5-7 cubes, 7,8 octahedra, 9,10 icosahedra, 11 dodecahedra, 6 buckyballs, 6 and nanoprisms. 7,12 Alternatively, the "DNA origami" or "DNA brick" techniques allow for the construction of an enormous range of shapes built from closely-packed helices. [13][14][15] Making many of these structures is only possible because of DNA's physiochemical properties, which can be manipulated and therefore controlled. Thus, understanding the biophysical and chemical properties of DNA is of vital importance for realizing future nanodevices and nanostructures.
A common motif used in many nanostructures is a 2-way junction in which two double helices are connected by a noncomplementary bulge loop of varying size. [16][17][18][19][20][21] Adjusting the bulge size in these 2-way junctions provides a way to control their flexibility. For example, a 3-arm DNA "star-tile" a) Electronic mail: john.schreck@chem.ox.ac.uk b) Electronic mail: jonathan.doye@chem.ox.ac.uk design, 22 illustrated in Fig. 1(a), contains a combination of 2-and 4-way structural elements linked together by doublehelical sections. [23][24][25][26] Where the arms meet at the center of the tile, the inner strand (red in Fig. 1(a)) contains bulge loops in order to facilitate the bending of the arms with respect to one another. When a solution of identical star tiles is prepared and cooled, the tiles can self-assemble by linking together at the single-stranded 'sticky' ends at the end of each arm. The assembly product is largely controlled by the bulge size. 12 For example, 3-arm tiles containing five nucleotides in each bulge region assemble into tetrahedra (illustrated in Fig. 1(b)), while the reduced flexibility of 3arm tiles with three nucleotides in the bulge regions leads to the formation of dodecahedra or buckyballs, depending on the tile concentration. 6 Another example is the DNA nanoprism illustrated in Fig. 1(c), where the sizes of the bulges between the different arms can be used to control the detailed geometry of the prism, in particular the relative twist of the top and bottom faces. 12 These examples illustrate that small changes in sequence design can clearly have a strong influence on the final product as well as the product structure. Even though bulged duplex structures are already widely used in nanodesigns, their effects on the structure and flexibility have not been thoroughly studied and the underlying principles have not been explored. Meanwhile, strong bending of DNA duplexes in the absence of bulges continues to be an area of considerable interest. 27 Early experimental studies showed that both DNA and RNA duplexes are subject to bending in the presence of bulges and internal loops, with bulges containing ∼5-7 nucleotides inducing larger bends than smaller bulges. 16,18,[28][29][30][31][32] Subsequent studies showed that the induced bend is not rigid and the junction can adopt a variety of conformations with similar bend angles, which signifies an increase in flexibility caused by the presence of a bulge. 33,34 Following initial computational studies of A-form RNA duplexes in two-way junctions, 35 Bailor et al. 21,[36][37][38][39] have recently applied three inter-helical (Euler) angles that describe the bending and twisting flexibility at the two-way junctions: one angle quantifies the bend at the bulge, and the other two describe the relative twist of each of the duplex arms meeting at the bulge. Their studies of the Protein Data Base (PDB) have shown that the geometric secondary structure of 2-way junctions restricts the overall 3D orientation of helices in RNA, 21 where these constraints arise from the steric and connectivity constraints imposed by the junctions. 38,39 The Euler angles provide a convenient and consistent way to characterize large numbers of configurations of 2-way junctions, which can be used to elucidate the level of influence the constraints imposed by the secondary structure have on the global conformation of DNA and RNA structures. In addition, several experimental groups have carried out similar analyses of two-way junctions in DNA, 29,40,41 where, for example, in recent work by Wozniak et al., 42 the three Euler angles for bending and twisting have been inferred from FRET experiments for bulged DNA duplexes with varying bulge size. However, even with these recent experimental advances, much of the literature discussing experimental and theoretical investigations of two-way junctions has to date focused primarily on RNA.
In this article, we use oxDNA, a coarse-grained model at the nucleotide level, to study in detail the structure and thermodynamics of B-DNA double helices containing a bulged loop, a motif that is henceforth referred to as a bulged duplex. The results of our simulations are compared with recent experiments in which the bend and twist angles of bulged duplex systems have been measured. 29,[40][41][42] The degree of coarse-graining in oxDNA allows us to study not just the equilibrium structure, but also the flexibility of these junctions.
The oxDNA model has been highly successful at reproducing structural, mechanical, and thermodynamic properties for single-and double-stranded DNA. 43 Moreover, applica- tions to study the fundamental biophysics of DNA, including the kinetics of hybridization, 44,45 toehold-mediated strand displacement, 46 the response to mechanical stress such as the over-stretching transition of dsDNA under tension, 47 the formation of cruciform structures under negative twist, 48 and the role of topology in the formation of kissing hairpin complexes, 49 have confirmed the robustness of the model. Furthermore, the model has proved useful in providing physical insight into the action of DNA nanodevices, such as nanotweezers, 50 and walkers, 51,52 and is starting to be applied to characterize large DNA nanostructures. 53 oxDNA captures the relative flexibility of single strands, which can adopt a variety of helical and non-helical structures due to a strand's ability to stack and unstack, 43,54 and the comparatively stiff duplexes. Because oxDNA simultaneously captures these thermodynamic and geometric effects for both single and double strands, it is well-suited to studying how the interplay of such fundamental factors shapes the overall behavior of bulges.

A. oxDNA model
In oxDNA, a single strand of DNA is modeled as a chain of rigid nucleotides, where each nucleotide contains one interaction site for the backbone and two more interaction sites associated with the stacking, coaxial stacking, and crossstacking interactions. The interactions in the model are illustrated in Fig. 2. Base-pairing interactions contribute to the overall potential energy only when bases obey Watson-Crick specificity, e.g. A-T or G-C pairs. In addition to base pairing and stacking, there are interactions associated with backbone connectivity and excluded volume. A detailed description of each interaction can be found in Ref. 43. Simulation codes for oxDNA are publicly available from the oxDNA website. 55 The model does have have some simplifications that are important for the current study. First, the model was fit to experiments using a salt concentration of [Na + ] = 0.5 m, where the electrostatic interactions are strongly screened. Even though some studies have discussed the effect of salt concentration on the overall flexibility of bulged duplexes, the concentration [Na + ] = 0.5 m is in the high-salt regime that is relevant to most DNA nanotechnology experiments, and which is our primary interest. Secondly, the oxDNA doublehelix is symmetrical with major and minor grooves being of the same size. While this may be important for certain motifs, we do not expect this approximation to change the dominant physics underlying bulge behavior, or the generic trends as a function of bulge size. Additionally, we choose to use the "average-base" parametrization of oxDNA, in which the hydrogen bonding associated with base pairing and the stacking interactions have the same strength independent of the bases involved, rather than its sequence-dependent parametrization. 54 This parametrization is advantageous for studying general properties of DNA unmodulated by sequence-dependent effects.

B. DNA bulged duplex systems studied
In our simulations we first consider a system where two DNA strands associate to form a bulged duplex, where each duplex arm flanking either side of the bulge contains 10 hybridized base pairs (red regions are complementary to each other, as are green regions), and the bulge region contains consecutive thymine (T) bases. This structure is very similar to the bulged duplexes found in the star tiles discussed above, 6 which contain 10 and 11 nucleotides in the arms, respectively. To investigate the effects of the size of the bulge loop on the structural properties we will consider ranging from zero to fifteen in our simulations. Additionally, the simulation temperature was set at 23 ∘ C, which is near the temperature (25 ∘ C) where the DNA nanoprism was found at high yield in experiments. 12 In addition, we also consider C and Z-tiles, 56 which are duplexes containing two bulge loops. The following strand can be used to form a Z-tile, where the red regions are complementary, as are the black regions. The central green region is a palindrome. Two of these strands can hybridize through the red and green regions, where the final structure contains one green duplex section, two red duplex sections, and 4 single-stranded black sections (two bulges and two sticky ends). Note that the green duplex section is about one turn in length. In general, the number of helical turns in the green duplex section determines the motif structure, where an even or odd number of half-turns yield Zor C-tiles, respectively. The red sticky-ends are specifically designed to be complementary to the nucleotides in the bulge loop. Under certain conditions, the sticky ends may bind with loops from nearby tiles and can form a T-junction, that is, a region in which 3 helical arms meet at a bulge that resembles a T shape, which facilitates the self-assembly of the Z or C tiles into 1D or 2D nanostructures, respectively. The Z-tile is further discussed in Section III B. The simulation temperature for Z-tiles was set at 22 ∘ C, the temperature at which the assembled tiles produced 1D nanostructures in high yield. 56

C. Simulation details
To calculate free energies of structures and investigate flexibility, bulged duplexes for several bulge sizes are simulated by employing virtual-move Monte Carlo (VMMC), as introduced by Whitelam and co-workers. 57 We use the variant introduced in the appendix of Ref. 57. VMMC is a clustermove algorithm that efficiently samples from the canonical ensemble for systems of strongly-interacting particles, which we found particularly useful for DNA. 53 The relative free energy of the bulged duplexes was sampled as a function of two order parameters: (1) the total number of base pairs in the system, and (2) the end-to-end distance, ee , defined as the distance between the center of mass of the bases at the 5 ′ and 3 ′ ends of the strand not containing the bulge. Two bases are considered to be paired when the hydrogen bonding interaction between them is 0.093 times its well-depth (0.596 kcal mol −1 at 23 ∘ C). Small variations in this cutoff do not significantly affect the results. OxDNA predicts that the base pairs at the ends of double helices may temporarily break, a process termed fraying, even for systems well below their melting temperatures. Since this phenomenon may occur at either terminal end of a bulged duplex, or in the middle near the bulge loop, the effects that fraying may have on bending can be monitored by following changes in the number of base pairs. The end-to-end distance, ee , is useful for monitoring how much the system bends, with strong bending occurring when ee is small compared to the length of a relaxed duplex.

D. Stacking and bulge classes
In our simulations, the bulged duplex can adopt a variety of conformations for different bulge sizes, some of which are illustrated in Fig. 3, where the stacking at the center of the duplex system may or may not be interrupted. For example, in Fig 3(a) the stacking interaction between the two bases that flank either side of the bulge, which are not neighbors along the sequence, is intact. The stacking interaction between the two bases directly opposite to the bulge, which are neighbors, is also intact and the system is mainly straight. In the model, the stacking interaction between the bases flanking the bulge is referred to as a coaxial stack to distinguish it from the stacking interaction between two neighboring bases, such as the stack opposite to the bulge. Fig 3(b) illustrates a configuration where a 1-base bulge has become inserted into the helix and stacks with its neighboring bases that are adjacent to the bulge. The stack opposite the bulge is intact and the system is mainly straight. In Fig 3(c), a configuration is shown where the coaxial stack is broken and a base from the loop is inserted into the helix while the rest of the bulge bases are outside of the helix. The stack opposite from the bulge is intact, and only marginal bending of the system away from the bulge is observed. In Fig 3(d), both the coaxial stack and the stack opposite to the bulge are broken and the system is bent away from the bulge. Other illustrated states of the system include a configuration in which fraying disrupts the base pairs flanking the bulge as well as interrupting stacking along the helix (Fig 3(e)), and a configuration which bends towards the bulge (Fig 3(f)). Similar observations are made for the Z-tile at each of the two bulged regions.
We see four main classes of conformations for the bulged duplex. The classifications are: A: The bases in the loop disrupt the duplex as little as possible and are flipped out, which results in almost no bending.

B:
The bases from the bulge loop are inserted into helix while maintaining stacking opposite the bulge, resulting in some degree of static bending away from the bulge.
C: A combination of A and B in which some bases are inserted into the helix, while others are flipped out.

D:
The stacks opposite and across from the bulge are broken, resulting in increased flexibility and a large static bend.
In order to apply these classifications, we focus on the stacking interactions between the pairs of bases near the bulge at the center of the system as well as base pairs that are adjacent to the bulge, as illustrated in Fig 4. To determine whether a base stacks or coaxially stacks with another base, a lower bound for the stacking interaction between two bases is defined to be the same as the hydrogen bonding cutoff value of 0.596 kcal mol −1 at 23 ∘ C. The choice of this value is somewhat arbitrary, however, small changes in this value do not change the results significantly.
In Fig 4(a), the state of the coaxial stacking interaction between bases across the bulge is denoted by the symbol , while the stacking states of the two bases on either side of the bulge are denoted −1 and +1 , respectively. Likewise, the state of the stack opposite the bulge is denoted , while the state of the stacks between neighboring bases on the same strand on either side of the bulge are denoted −1 and +1 , respectively. Finally, the state of the stacks between the bases with squares and the bases with triangles on either side of the bulge are denoted 1 and +1 , respectively, and are referred to as "stem-loop" stacks, while the stacks between the bases in the loop are denoted 2 , 3 . . . , . If the magnitude of the stacking interaction is greater than or equal to the lower bound, = 1, otherwise the stack is taken to be broken and = 0, and similarly for and . If the stack opposite the bulge, , is present, we take the two duplex arms to be "stacked" and group these configurations into a stacked set. Configurations where is disrupted are taken to be "unstacked" and are grouped into an unstacked set. According to our classification scheme, a stacked configuration may fall into either A, B or C, while an unstacked configuration is classified as D.  stacks from one stem-loop stack to the other, i.e. 1 . . . +1 ; and C if greater than one stack in the loop is broken but stacking is intact at one or both of the stem-loop stacks (more than one of 2 , . . . , broken and at least one of 1 or +1 intact). We note that configurations that have frayed base pairs at the junction can be difficult to classify using our scheme. In the supplemental materials we discuss several extensions to the classification scheme to properly deal with frayed configurations. 58

E. Bend and twist angles of bulged duplexes
We define the bend and twist angles using a scheme similar to that used in Ref. 21. To measure the bend angle, , for a given configuration we place unit vectors labeled 1 and 2 along the helical axes of each of the duplex arms flanking the bulge. These vectors are illustrated by blue arrows in Fig. 4(a) and in more detail in supplemental Fig. S1(a). The duplex arms have the freedom to twist about the vectors 1 and 2 as characterized by the angles and , respectively, as illustrated in Fig. S1(a). 58 The vectors 1 and 2 are defined by finding the longest stretch of base pairs in each arm, and then drawing a line from the center-of-mass of the base pair at one end of the duplex arm to the center-of-mass of the base pair at the opposite end of the same duplex arm. For normalized 1 and 2 , the bend angle is then calculated using a definition that follows the convention for the bend angle in the recent literature. 21 In the supplemental material, 58 we develop a simple convention to determine whether the system is bent away from (0 ∘ < ≤ 180 ∘ ) or bent into the bulge (−180 ∘ ≤ < 0 ∘ ). Lastly, we explicitly define the duplex twisting angles and . Each angle can be calculated by first computing a vector that points from the base flanking the bulge to its complementary partner directly across from the bulge. The two vectors are referred to as 1 and 2 and point from square to circle in Fig. 4(a) and are illustrated in Fig. S1(a). The twist angles and can be calculated using respectively, where = 1 × 2 is a vector normal to the plane of the bulged duplexes. Similar to , we develop a convention for determining when the angles and take on the values between 0 . This is discussed in section S1 in the supplemental materials. The relative twist between the duplex arms flanking the bulge is taken as − . For reference, in a relaxed duplex as represented by oxDNA, − ≈ 32 ∘ , i.e. the twist per base pair rise in a duplex.

III. RESULTS
A. Bulged duplex systems

Stacking classes
We first discuss how the balance between structural classes A, B, C, and D changes as a function of the length of the bulge, , in oxDNA. Each bulged duplex system with a given bulge length was simulated at least 10 times and configurations were collected until the standard error of the mean for the points computed in the average bend angle versus bulge size was less than 1% of the the computed mean value. The results are shown in Fig. 5.
Systems with small bulge loops ( = 1 and = 2) are dominated by class B. In this scenario the bases in the bulge are inserted into the helix, as was illustrated in Fig. 3(b). Insertion is favored over flipping the bases outside of the helix because more stacking interactions are preserved. For = 3 it is no longer feasible to insert 3 bases into the duplex and maintain stacking opposite the bulge. Therefore, the probability of B decreases dramatically, and instead configurations where one or two of the bases in the bulge are inserted into the duplex, while the others predominately flip out, become prevalent. Additionally, = 3 is the first case where unstacked states make a significant contribution. For > 4, class D (e.g. unstacked configurations) becomes the most probable class and of the stacked states, A becomes increasingly dominant over class C as increases. Also very apparent is that the probability of stacked states goes through a maximum at = 6 before plateauing off at 11 where cases A and D are roughly equally probable with some contributions still arising from case C.
The most obvious trend is that although unstacked configurations are rare for small , they become common at larger . In stacked configurations, the endpoints of the bulge loop are strongly constrained; unstacked configurations offer greater freedom. Whether or not unstacking occurs depends on the relative benefit of this freedom compared to the cost of disrupting junction stacking.
Loops with small do not benefit as much from the additional freedom of the unstacked ensemble for two reasons. Firstly, they are too short for all relative orientations of the duplex arms to be explored. As the loop gets longer, the duplex arms have more freedom and the unstacked state becomes more favorable. Secondly, steric penalties associated with the stacked configuration are smaller for small .
The bulged duplex systems that we study are somewhat similar to a system with one or two single-stranded dangling ends at a junction in a duplex, except that in our cases the strands are connected and form a bulge loop. Duplexes with such dangling ends can occur during toehold-mediated strand displacement when an invader strand is displacing an incumbent strand during branch migration, as described in Ref. 46. In these systems, once the invading strand is bound to the toehold, there is a free-energy penalty to initiating displacement, even though the number of base pairs is unchanged, because it is unfavorable to have two single-stranded overhangs at the junction. The free-energy penalty arises from the overcrowding of nucleotides at the junction, and saturates once both overhangs have at least 3 or 4 nucleotides, because further bases are sufficiently far enough away from the junction that their contribution to overcrowding is minimal.
Similar to the displacement system, for the bulged duplex to maintain a stacked state (i.e. one of cases A, B, and C) the bulge loop must arrange itself to minimize steric clashing amongst the bulge bases and also minimize duplex disruption caused by inserted bulge bases which can compete with the coaxial stack. Alternatively, a bulged duplex may unstack and bend away from the bulge gap, freeing the bulge bases to spread out into space and decrease overcrowding. As with displacement intermediates, the benefit of spreading out is greater when there are more nucleotides at the junction, providing the second cause for the the increase in unstacking with increasing .
Several aspects of the data, however, are not explained by this analysis. Firstly, why does case C become less favorable compared to case A as loop length increases? Secondly, why does case D increase with respect to case A, and then subsequently decrease (which will subsequently give a nonmonotonic variation of the bend angle as shown in the next section)? Both of these questions can be understood in terms of how changes in configuration at the junction are related to the typical physical length of the bulge. In order for the duplex arms to be stacked (in class A, B or C), any bases not incorporated into the helix must be compressed so that the single-stranded loop region adopts a conformation with a short end-to-end distance. When the bulge loops are short ( ≤ 4), inserting one or two bases into the duplex helps to reduce the compression substantially: a 2base loop is less constrained than a 3-base loop because it needs to be compressed less. As the loop gets larger, however, this difference becomes less substantial and so it is not as advantageous for a stacked duplex to incorporate bases. Thus case C is favorable for short loops, but less so for longer loops.
Next we address the competition between cases A and D. Medium-sized loops (i.e. 5 ≤ ≤ 10) benefit more from unstacking than longer loops ( ≥ 11) which resist compression less. This is because medium sized loops are only slightly longer than the persistence length of ssDNA, and unstacking of duplexes at the junction allows them to stretch out more easily. By contrast, for longer loops the relative cost of being bent is lower, and so they benefit a lot less from being unstacked at the junction.
This typical behavior of polymers can be seen, for example, in the end-to-end probability distribution of a worm-like chain. 59 The analytical formula for the probability density of the end-to-end distance in Ref. 59 clearly shows that a shorter polymer benefits a lot more from having its end-toend extension increased by a fixed absolute distance when compared with a longer polymer, as illustrated in Fig. S3, where chains with contour lengths that are two and four times the persistence length are compared. Similar results can also be obtained for a freely-jointed chain.
We checked whether the non-monotonic behavior was peculiar to our model of ssDNA, or a generic polymer effect, by switching off nearest-neighbor stacking interactions between loop nucleotides (while maintaining them in the stem and at the interface of stem and loop). The results, shown in Fig. 6(a), which plots the free-energy difference between the combined classifications A, B, C and classification D, show that stacking between bases in the loop does not significantly change our results and therefore that our results are a robust consequence of generic polymer behavior. (a) The relative difference in free energy between unstacked and stacked configurations. us and s are the probabilities that a bulged duplex is found in the unstacked and stacked states, respectively. (b) The average bend angle for stacked (red circles) and unstacked (blue squares) populations is plotted against bulge size . The average bend angle for the complete set of configurations is also plotted (black triangles). In both figures the error bar on each data point are smaller than the symbol size. Fig. 6(b) shows the average bend angle ⟨ ⟩ as a function of bulge size, as well as the bend angle for stacked (classes A, B, and C) and unstacked states (class D). The average bend angle for a bulged duplex for small bulge sizes is small when most configurations are stacked, but it quickly grows with as unstacked configurations become more favorable. The average bend angle peaks at 80 ∘ at = 7 before decreasing again and finally plateauing at about 55 ∘ for ≥ 11, this behavior mainly reflecting the variation of the probability of being unstacked (Fig. 6(a)).

Effects of stacking on bending angle
The average bend angle for unstacked states increases monotonically from = 1 and reaches a maximum at = 8 where ⟨ ⟩ = 117 ∘ , but then starts to decrease as grows, eventually leveling off at ⟨ ⟩ ≈ 105 ∘ for ≥ 11. The drop in the average bend angle before leveling off is also a signal of the compression effect. The configurations belonging to the stacked population start off slightly bent at 20 ∘ due to the high probability that the coaxial stack is broken because a few bases have inserted into the helix, but the bending angle tends towards 10 ∘ for longer loops because there is less incentive for the bases to insert into the helix.
Another noteworthy feature is that the bend angle in stacked and unstacked junction configurations is essentially independent of the nearest-neighbor stacking between the bases in the loop (see Fig. S2). 58 The interactions at the stem-loop interface, however, are important. These stemloop stacks influence the bending angle when configurations are classified as D (configurations in which duplex arms are unstacked at the junction). We demonstrated this by switching off the nearest-neighbor stacking interactions between bulge loop bases, the stem-loop stacks, the coaxial stack, and the stack opposite to the bulge. We compared this system to systems in which the stem-loop interactions were not removed. The results, presented in Fig. S2, show that for unstacked configurations the stem-loop interactions cause the system to bend significantly more than if the stacks were not present. 58 This is because such stem-loop interactions direct the loop bases to carry on in the direction of the stems, thus requiring a larger bend angle to avoid steric repulsion, as is clear from the configuration in Fig. 3(d).
For a more detailed look at how changes in bulge size influence the flexibility of bulged duplexes, in Fig. 7 we plot the bend angle and end-to-end distance distributions for various . For stacked configurations there is a clear asymmetry between bending into, and away from the bulge loop, with configurations that bend away from the bulge more favored because these configurations reduce steric clashing between the bulge loop and the duplex. This asymmetry is most pronounced for small where base insertion is prevalent.
The separation of the stacked and unstacked states is clear from the scatter plots in Fig. 7. As unstacked configurations have a slightly increased effective contour length relative to their stacked counterparts, the two distributions overlap more when projected onto the end-to-end distance, ee , than . The range of bend angles that are available to unstacked states becomes wider as the bulge size increases, and underlies the increase of ⟨ ⟩ with bulge size for = 1 − 8. Comparing the bend angle distributions for = 8 with = 10 for the unstacked states, it is clear that the distribution for = 10 is less sharply peaked than the distribution for = 8, and underlies the slight decrease in the average bend angle ⟨ ⟩ for the unstacked states when increasing the bulge size from = 8 to = 10.

Effects of stacking on duplex twist angle −
To quantify the relative twisting of the duplex arms, − is calculated for the same set of configurations. In Fig. 8, the probability that a configuration occupies a state with a given value of ( , − ) is plotted as a 2D histogram. As with the analysis of the bend angle, , we split the configurations into stacked and unstacked sets. A very pronounced (and expected) feature common in all bulged duplexes we considered is a stacked duplex population with a relative twist that is found to lie in a very narrow range, centered at approximately 32 ∘ , and is due to configurations in class A.
As illustrated in Fig. 8, once a bulge is introduced into the Data is plotted for a selection of bulge sizes in the range = 0 to = 10. In the scatter plots, red crosses indicate unstacked states and blue squares represent stacked states. The probability plots retain the same color scheme as the scatter plots. Equivalent plots for other bulge sizes can be found in the supplemental materials. 58 duplex, the relative twist at the bulge maintains a signal at 32 ∘ for all values of studied, however a second signal at approximately − ≈ 70 ∘ appears for stacked configurations that are mainly bent away from the bulge. These configurations mostly correspond to the scenario where some number of the bases from the loop are inserted into the duplex at the bulge, as illustrated in Fig. 3(b) and (e). As we had shown earlier for all bulged duplexes, one, and sometimes two bases can be inserted into the duplex. Clearly, extra bases at the center of the duplex cause the bulged duplex to twist more at the bulge than stacked configurations belonging to class A. The spread for the relative twist angles in these states is also quite narrow and the effects due to base insertion is greatest for = 1 loops when the base in the bulge is much more likely to be found inserted rather than flipped outside of the duplex. Base insertion significantly persists up to = 4, then the effect decreases for > 4 and plateaus. These plots nicely illustrate the fact that there remains a finite probability at large that a base from a longer loop can still become inserted into a duplex, as can be seen in Fig. 3(c). Fig. 8 also shows some interesting features of the unstacked bent states. Initially, when the unstacked population first appears with a significant probability compared to the stacked state ( = 4), the relative twist angles roughly fall in a circular distribution in the − ( − )-plane, centered around ≈ 70 ∘ and − ≈ 100 ∘ . However, deviations from the center are modest and are increasingly less probable due to the constraint of the short loop. There is also some overlap between stacked and unstacked populations for = 4, which corresponds to those configurations where the coaxial stack is broken, but not the stack opposite the bulge, and those configurations where both stacks are broken, but the bending angle is still constrained by the geometry of the loop, respectively. Upon increasing the bulge size further from = 6, the circular distribution elongates along the − axis, because the longer loop allows the bulged duplex arms considerably more freedom to twist relative to each other. Continuing this trend, once the bulge size increases up to ≥ 11, the duplex arms can be oriented at almost any relative twist angle with relative ease.
FIG. 8. The probability that a configuration assumes some value of ( , -) is plotted for various bulge sizes . Configurations were separated based on the state of the stacking between the double stranded regions, with the stacked set in blue and the unstacked set in red. The probability distributions for the stacked and unstacked sets were normalized to one separately to enhance clarity. Equivalent plots for other bulge sizes can be found in the supplemental materials. 58

B. Z-tile structure
As mentioned in the Section II, the Z-tile, illustrated in Fig. 9(a), is a symmetric DNA building block containing two bulge regions, three short duplex regions, and two sticky regions at the terminal ends of the tile that are complementary to the bulges. The bulge size can be exploited to control the angle in between the duplex regions. Under certain conditions, many identical copies of the Z-tile may link together where the sticky regions of one tile hybridize with the loops from another tile to form T-junctions. 56 An example of a 1D ribbon built using oxDNA and containing 20 Z-tiles is illustrated in Fig. 9(b). Assembled T-junctions can be seen in the structure in Fig. 9(c).
The bending and twisting angles for each bulged region in the Z-tile, as well as end-to-end distance for each duplex region flanking a bulged region, can be defined similarly to the same quantities for the bulged duplex system. A junction is defined to be unstacked if the stack opposite the bulge is broken, and stacked otherwise. We use the stacked/unstacked convention to determine if the Z-tile has 0, 1, or 2 bends. The angles and end-to-end distances are defined and illustrated in Fig. S1 in the supplemental materials. 58 We show the results for the angles by computing the average bend angle at the bulge in each of the two bulge regions, 1 2 ( 1 + 2 ), and the total average twist measured at the bulges as 1 2 ( 1 − 1 + 2 − 2 ). The results for the bend and twist angles for the experi-mental Z-tile that we consider here, which contains 5 bases in each bulge loop, are plotted in Fig. 10. The plots illustrate that the bending and twisting features of the Z-tiles are similar to the same features found in bulged duplexes, namely, a Z-tile can be approximately straight (indicated by "I" in the figure) with a frequency of 16%, one arm can be bent while the other arm is approximately straight and vice versa ("L") with a frequency of 47%, and finally both arms can be bent ("Z") with a frequency 37%. Predictably, the "Z" configurations also display the most amount of flexibility as is evidenced by the widths of the distributions for the bend and twist angles, while the bend and twisting angles for the "I" configurations cluster more tightly around an average value when compared with the "Z" configurations. Thus, when both arms are straight, the structure is quite stiff, but it can gain a significant amount of flexibility when the system breaks the stacks opposite to both of the bulge loop regions. The configurations where both arms are bent have a similar geometry to Z-tiles that have self-assembled into larger structures.

C. Comparison with experiments
FRET experiments carried out by Woźniak et al., 42 yielded bend angles of 32 ∘ , 56 ∘ , and 73 ∘ , for = 1, 3, 5 respectively, which are broadly consistent with previous experimental results. 29,40,41 We see the same trend in the average bend angle predicted by oxDNA, namely, 18 ∘ , 20 ∘ , and 60 ∘ for the same bulge sizes, but the angles are noticeably smaller. Interestingly, the experimental values are actually closer to the average oxDNA bend angle for unstacked states, namely, 21 ∘ , 58 ∘ , and 88 ∘ , for = 1, 3, 5, respectively. This perhaps suggests that oxDNA somewhat overestimates the probability of being in the unstacked state, and could be due to overestimating the coaxial stacking strength or a neglect of local geometric factors that hinder stacking across the bulge, although we note that there are also significant uncertainties in the experimental measurements. Another possibility is that we underestimate the bend caused by base insertion, or the number of bases that can be inserted. A recent study of bulged duplexes in DNA observed an increase in flexibility of the motifs (a wider range of bend angles was observed) upon an increase in bulge loop size, findings which are broadly consistent with the predictions of oxDNA for the average bend angle distributions. 60 Interestingly, Bailor and co-workers 21,39 have reported results for RNA where the bend angle in bulged RNA duplexes monotonically increased up to about = 7 and then plateaued. This is again similar to the behavior we see for the oxDNA bend angles for the unstacked configurations which rises and then plateaus at = 8. Our results are also very relevant to the self-assembly of multi-arm star tiles into polyhedra that was mentioned in the Introduction, and examples of oxDNA representations of these structures are shown in Fig. 1. In particular, the assembly products can depend sensitively on the bulge size with larger bulges favoring structures with larger bend angles at the vertices, 6 and the bulge size can also be used to influence the angles at the vertices of the final polyhedra when not fully constrained by the topology of the polyhedron, as is the case for the prism in Fig. 1(c). Firstly, our results show that the free-energy cost of breaking stacking at the bulge, which is necessary in the assembly of polyhedra, is relatively small for ≥ 3 (i.e. < 2 B ), and in particular is small compared to the free-energy gain from the basepairing associated with inter-tile assembly. Secondly, as the bulge loop gets longer, the range of bend angles available clearly increases until a point is reached where the bulge is sufficiently big that all bend angles are feasible. Simply put, larger bulges give rise to greater flexibility once the stacking at the junction is broken. For example, in the unstacked state, the free-energy cost of having a bond angle of 120 ∘ (as is required for the triangles in the tetrahedron) is significantly smaller for the = 5 than the = 3 bulge (by roughly 4 B ), and so helps to explain why tetrahedra are not formed from 3-arm tiles with bulges of size = 3, but only for = 5. In future work we intend to study how bulge size explicitly controls the rate of closure of linear trimers to forming triangles. Similarly, in the nanoprism the lower freeenergy cost of bending for larger bulges can explain how the bulge size can be used to control the angles at the vertices and hence the relative twist of the two triangular faces. A more extreme example is that 3-, 4-, 5-, and 6-arm motifs have been found to fully dimerize to form nanotubes with well-defined diameters and lengths when the bulge size is large ( = 9), and when the tile contains hairpins in the tile arms, which further increases flexibility. 61 In such structures, the bend angle is likely to be close to 180 ∘ .

IV. CONCLUDING REMARKS
We have studied the structural properties of bulged duplexes and Z-tiles using the coarse-grained oxDNA model, with the aim of understanding and characterizing the flexibility of these basic motifs that are widely used in DNA nanotechnology experiments as a means to tune the selfassembly and the equilibrium structure of the final product. We find that bulged duplexes typically adopt one of four configurations, three of which involve the duplex arms stacking at the junction and one of which does not. When loop sizes are small, the bulged duplex systems mainly prefer to be in a stacked state with some or all of the bases in the loop inserted into the helix. Base insertion has the effect of increasing the amount of twist in the duplex because of the added stacking interaction sites at the center of the system. The resultant average bending angle for small loops is modest and does not deviate too far away from the average bend angle of duplexes in the absence of bulges, which is about 10 ∘ to 20 ∘ .
Once the bulge loops are about four or more nucleotides in length, configurations are increasingly found to be unstacked. In both stacked and unstacked configurations, the loop tends to be found on one side of the bulged duplex and with either the coaxial stack alone broken, or both stacks at the center of the system broken. The latter case allows the system to access greater bending angles. The model also predicts that medium sized bulges that are on the order of the persistence length of ssDNA will resist compression more strongly than longer bulges. This effect causes the system to significantly favor the unstacked configurations over the stacked configurations, with the unstacked configurations exhibiting a large static bend that is partially driven by stacking of the bulge loop with the stems. However, the range of the twisting between the duplex arms is somewhat restricted. Systems with large bulge sizes are less affected by a constrained loop, and were found to have the greatest flexibility in which the bending angle can assume values over a wide range, and also the relative twist between two duplex arms are free to take on nearly any value from 0 ∘ to 360 ∘ . We also studied the Z-tile and found similar bending and twisting features that were seen in the duplex system. oxDNA is a coarse-grained model that was derived to represent generic properties of DNA in a computationally efficient way. As such, it would not be expected to quantitatively reproduce all data from experiments on bulges, to which it was not explicitly parameterized. It does, however, serve to highlight the underlying physics that is of relevance to real systems. In particular, we would expect the four distinct classifications observed to be robust. Further, we have identified generic factors that drive changes in occupancy and properties of these configurations with bulge length. These factors are related to basic polymer properties and geometrical/steric constraints, and therefore will also be relevant in experiment.
The properties that we have observed for bulged duplexes also provide insights into the use of this motif in DNA nanotechnology. In particular, the relatively small free-energy cost for unstacking for ≥ 3 and the greater flexibility in the unstacked state allowed by increasing bulge size can help to rationalize further the design rules for controlling the selfassembly product and final structure of polyhedra assembled from multi-arm star tiles. As is illustrated in Fig. 1, oxDNA is efficient enough to allow us to study these nanostructures, and understanding their self-assembly and structure will be the subjects of future work.
1 and 2 were defined in Section IIE in the main text. If the quantity is greater than zero, we let → 360 − . Similarly, if 2 · ( 2 ×^) < 0, we let → 360 − . Despite these simplistic methods for determining the sign of , and , our simulations show that the results are reasonably reliable.

Appendix SII: Effects of Fraying on Stacking Classes
Here we explain how we take into account the effect of frayed base pairs near the bulge at the center of the duplex on our classification scheme. Even though most configurations analyzed are not frayed at the center, those that have become frayed can influence the interpretation of the results. If only the states of the stacks and are considered, our convention, as it stands, incorrectly associates some frayed configurations with stacked states when they should be counted as being unstacked (for example, the configuration in Fig. 3(e) in the main text would be counted as stacked). Additionally, fraying also interferes with our ability to properly to sort configurations into the classes A, B, C, and D.
To assist in properly sorting most frayed configurations, we instead consider the system to be "stacked across the bulge" when −2 −1 +1 +2 = 1, and the system is said to be "stacked opposite the bulge" whenever −2 −1 +1 +2 = 1 (these symbols are defined in the main text). By comparing with systems where fraying was prevented, we have found that this scheme accurately captures the state of the stacking at the center of the system in the presence of frayed base-pairs. Additionally, the configurations that are classified as B are subjected to a second test in which the status of the following stacks are checked: FIG. S1. In (a), the cylindrical representation of the duplex arms meeting at the bulge is illustrated, with the twist angles and , the duplex unit vectors 1 and 2, and the twist vectors 1 and 2 all illustrated. The two coordinate systems are defined by crossing the normal vector = 1 × 2 with each vector 1 and 2 yielding 1 and 2, respectively. (b) Schematic illustration of a Z-tile where cylinders represent duplex regions. The end-to-end distances ee,1 and ee,2 that were used in calculating the relative free energy are also illustrated. Note that the unit vectors are not drawn to scale. +1 , +2 , +3 , +4 . This test ensures frayed configurations are properly sorted. If more than one of these stacks is broken, the configuration is reclassified as class C, otherwise it is left classified as B.

Appendix SIII: Loop Stacking Effects
In order to determine the role that stacking between bases in the loop had on system flexibility we ran separate simulations for several bulge sizes in the three main regimes seen in the results sections: (i) small loops, where bases are often found inserted into the duplex, and where the stacked states are significantly more prevalent than unstacked states, (ii) medium sized loops, where unstacked states are more prevalent than stacked states, (iii) long loops, where there is no preference for stacked states over unstacked states, and vice versa. To accomplish this, we ran simulations for = 2, 6, 10. In the simulations, we used the sequence dependent version of the model, 54 which allowed us to introduce a dummy nucleotide that was used to replace the nucleotides in the loop. We were then able to switch off any stacking interaction between two dummy bases. All other interactions between a dummy nucleotide and either A, T, G or C were left unchanged. This means that in these simulations, the stem-loop stacks may still occur. The results from the simulations are plotted in Fig. 6(a) in the main text where the free-energy difference between stacked and unstacked configurations is illustrated for several bulge sizes. When comparing the results for the system where stacking is pre- served in the loop with the system where stacking is switched off in the loop, we conclude that the stacking between bases in the loop only marginally affects the probability of being unstacked at the junction, and does not significantly affect the bend angle (see Fig. S2).
The entropic repulsive effect of the stem-loop stacks on bend angle was also considered in the main article. In the simulations, the stacks 1 · · · +1 , , and were all switched off. These results are shown in Fig. S2. There is a clear reduction in the bend angle compared to that for the unstacked states in Fig. 6(b) in the main text (which are replotted in Fig. S2) and no maximum in the bend angle at intermediate values of .
FIG. S3. The interpolated probability density is plotted for wormlike chains with lengths equal to two (squares) and four (circles) times the persistence length, where = / where is the persistence length and is the contour length of the chain. With this convention, the end-to-end separation of the chain ranges from zero to one.
where the vectors 1 , 2 , 1 , and 2 are illustrated in Fig. S1. The relative twist angle from one arm to the other is defined by projecting vectors 2 and 2 into the plane that is orthogonal to the vector 1 − 1 . The projections of 2 and 2 are denoted 2 and 2 , respectively.
FIG. S4. In each figure there are three plots that illustrate the relationship between bend angle and end-to-end distance ee: a scatter plot in ( ee, ) for sets of configurations generated in our simulations, the probability of a state occupying a particular value of ee (top panel), and the probability of a state occupying some value of the bend angle, (right panel). Data is plotted for a selection of bulge sizes in the range = 0 to = 10. In the scatter plots, red crosses indicate unstacked states and blue squares represent stacked states. The probability plots retain the same color scheme as the scatter plots.
FIG. S5. The probability that a configuration assumes some value of ( , -) is plotted for all bulge sizes studied, = 1 − 15. Configurations were separated based on the state of the stacking between the double stranded regions, with the stacked set in blue and the unstacked set in red. The probability distributions for the stacked and unstacked sets were normalized to one separately to enhance clarity.