DNA hairpins primarily promote duplex melting rather than inhibiting hybridization

The effect of secondary structure on DNA duplex formation is poorly understood. We use a coarse-grained model of DNA to show that specific 3- and 4-base pair hairpins reduce hybridization rates by factors of 2 and 10 respectively, in good agreement with experiment. By contrast, melting rates are accelerated by factors of ~100 and ~2000. This surprisingly large speed-up occurs because hairpins form during the melting process, stabilizing partially melted states, and facilitating dissociation. These results may help guide the design of DNA devices that use hairpins to modulate hybridization and dissociation pathways and rates.

OxDNA and its interaction potentials have been described in detail elsewhere [1][2][3]. The model represents DNA as a string of nucleotides, where each nucleotide (sugar, phosphate and base group) is a rigid body with three interaction sites. The potential energy of the system can be decomposed as where the first sum is taken over all nucleotides that are nearest neighbors on the same strand and the second sum comprises all remaining pairs. The interactions between nucleotides are schematically shown in Fig. 1 in the main article. The backbone potential V b.b. is an isotropic spring that imposes a finite maximum distance between backbone sites of neighbors, mimicking the covalent bonds along the strand. The hydrogen bonding (V HB ), cross stacking (V cr.st. ), coaxial stacking (V cx.st. ) and stacking interactions (V stack ) are anisotropic and explicitly depend on the relative orientations of the nucleotides as well as the distance between the relevant interaction sites. This orientational dependence captures the planarity of bases, and helps drive the formation of helical duplexes. The coaxial stacking term is designed to capture stacking interactions between bases that are not immediate neighbors along the backbone of a strand. Base and backbone sites also have excluded volume interactions V exc and V exc . Hydrogen-bonding interactions are only possible between complementary (A-T and C-G) base pairs. In the sequence-dependent parameterization that we use for all simulations, the strengths of interactions V stack and V HB further depend on the identity of the bases involved [3].
Interactions were fitted to reproduce melting temperatures and transition widths of oligonucleotides, as predicted by SantaLucia's nearest-neighbor model [4]. Note that our approach is significantly more complex than the nearest-neighbor model. We simply treat the latter as a high-quality fit to experimental data. For the purpose of parametrization, structural and mechanical properties of both double-and single-stranded DNA are also taken into account in the fitting procedure. In DNA the double helical structure emerges because there is a length-scale mismatch between the preferred inter-base distance along the backbone, and the optimal separation of bases when stacking. It is exactly this feature that drives the helicity of oxDNA, rather than an imposed natural twist on the backbone. Overall, the emphasis in our derivation of oxDNA was on physics relevant to the duplex formation transition. As discussed in the main text, oxDNA has been extensively tested for other DNA properties and systems to which it was not fitted. Our success in describing all these phenomena gives us confidence to use it to study the dynamics of hybridization in the presence of hairpins.

SII. SIMULATION METHODS
A. Thermodynamics

Virtual Move Monte Carlo
A standard approach for calculating thermodynamic properties of computational models is the Metropolis algorithm [5]. A drawback with this approach is that only moving single particles at a time results in slow equilibration for systems with strong attractions. This is true for DNA strands, where collective diffusion is strongly suppressed if nucleotides are moved individually. Simulations can be made more efficient by using the Virtual Move Monte Carlo (VMMC) algorithm proposed by Whitelam and Geissler which allows for collective diffusion using cluster moves of particles [6]. Specifically, we use the variant presented in the appendix of Ref. 6. Initially a particle is selected, and a move is chosen at random as in the Metropolis algorithm. The particle's neighbors are then added to a co-moving'cluster' with probabilities determined by the energy changes that would result from the move. Consequently, multiple particles tend to move at once. To use VMMC, we must select 'seed' moves of a single particle. For all VMMC simulations reported here, the seed moves were: • Rotation of a nucleotide about its backbone site, with the axis chosen uniformly on the unit sphere and the angle drawn from a normal distribution with a mean of zero and a standard deviation of 0.22 radians.
• Translation of a nucleotide, where the displacement along each Cartesian axis is drawn from a normal distribution with a mean zero and a standard deviation of 0.15 simulation units of length (0.1277 nm).
To improve efficiency, if the algorithm generates a cluster move involving more than 7 particles the move is automatically rejected.

Umbrella Sampling
An important concept is that of a reaction coordinate (or order parameter) Q, which groups together microstates of a system that share some macroscopic property (for example, all configurations of strands with a certain number of base pairs). The free-energy profile as a function of Q can provide useful information about the reaction, provided an appropriate choice has been made. Free-energy barriers can make certain regions of configuration space hard to reach, which prevents efficient sampling of all of the states of interest. The freeenergy landscape can be artificially flattened by weighting states with different values of Q appropriately, a technique known as umbrella sampling [7]. Thermodynamic properties of the system can then be extracted from simulations by unweighting the resulting distributions.
In particular, for an unweighted simulation a particular microstate with coordinates q N and energy E(q N ) is sampled with probability The equilibrium average of some variable A(q N ) is then given by the sum over all states, weighted by their Boltzmann factors: By applying a weighting w = w(Q(q N )) to each value of the order parameter, we change the sampling frequency to where the subscript w indicates a property of the weighted system. So we can artificially ensure that our simulation samples all states equally by making P w constant for all microstates. Equilibrium thermodynamic properties are then obtained by unbiasing afterwards, as can be seen by rewriting Eq. (S3) as follows: Throughout this article, it makes sense to use the number of base pairs in our definition of Q. This is the usual choice for studying hybridization processes.

Single Histogram Reweighting
To determine melting temperatures for structures, we implemented the temperature extrapolation method known as single histogram reweighting, based on the method introduced by Ferrenberg and Swensden [8]. The method of single histogram reweighting allows extrapolation of results from simulations at a particular temperature T 0 to other temperatures [8]. States of the system are grouped by their value of some quantity A and their energy E, so a histogram p(A, T 0 , E) can be produced. The temperature-independent density of states Ω(A, E) can then be inferred via where β 0 = 1/RT 0 . The proportionality constant is unknown because we can only ever know the relative ratios of states in our simulations. Ω(A, E) can then be used to calculate the average value of A at any temperature T by integrating over all possible states: We can rewrite this using Eq. (S7) as We point out that our potential energy function (Eq. (S1)) depends explicitly on the temperature through V stack . It is straightforward to extend Eqs. (S8) and (S9) to our case. It must be pointed out that the extrapolation can only go so far from the temperature T 0 because sampled states will not be representative of the dominant states at other temperatures. For this reason we explicitly calculate free energies at T = 20 • C for the main data of the manuscript and in subsequent sections here. We only use histogram reweighting for melting curves in the vicinity of the T m for a particular system, which are of less importance. In all VMMC simulations employing umbrella sampling (for single strand and double strand systems), we estimated the error of the computed relative free energies corresponding to different states by computing the standard error of the mean value of multiple independent simulations. The details of each simulation, including the order parameters used as well as the number of independent simulations, are discussed in Section SIII A.
B. Kinetics

Molecular Dynamics
Kinetic simulations were performed using an Anderson-like thermostat, similar to the one described in appendix A of Ref. 9. The Newtonian equations of motion for the system are integrated by Verlet integration [10] with a discrete time-step δt, so that the positions, velocities, orientations, and angular velocities of the nucleotides are recalculated at each time-step. This alone would give the DNA strands constant energy and cause ballistic motion. In reality, DNA in a solvent is being bombarded by water particles and thus undergoes Brownian motion. To model Brownian motion, the velocity of each nucleotide is resampled with a probability p v = 0.02 from a Maxwell-Boltzmann distribution at the temperature of the solvent every 103 time steps. The algorithm also resamples angular velocities with a different probability p ω = 0.0067. The solvent thus acts as a large heat bath at a fixed temperature, ensuring that the simulated system samples from the canonical ensemble. On time scales longer than N N ewt δt/p v , where δt is the integration time step, the dynamics is diffusive. We choose δt = 1.52 × 10 −14 s for all dynamics simulations in this study. In oxDNA this time step gives a diffusion constant D sim for a 14-mer duplex that is about 100 times higher than experimental measurements [11] of D exp = 1.19 × 10 −10 m 2 s −1 .
This artifical increase in D sim is a common procedure for coarse-grained models where higher diffusion constants can be used to accelerate diffusion. Accelerated diffusion can also speed up certain processes by smoothing out, on a microscopic scale, energy profiles [12]. This can be advantageous because it means simulations utilizing coarse-grained models can be used to study more complicated systems. In a previous study using oxDNA, the hybridization kinetics of a non-repetitive sequence was considered [13]. In that study it was shown that using higher friction constants (smaller diffusion coefficients) in simulations utilizing Langevin dynamics at 300K slowed down hybridization, but did not otherwise qualitatively affect the results. In particular, the tendency for initial base pairs to melt away rather than lead to a full duplex was found to be preserved. Our systems are similar to those studied in Ref. 13, possessing similar numbers of total base pairing between the strands, and using a similar simulation temperature. Additionally, many approximations of real DNA have already been made in the construction of the oxDNA model, and we expect that running simulations with a diffusion coefficient that is larger than the experimentally measured value should preserve the effects that hairpins in single strands have on the relative hybridization and dissociation rates.

Forward Flux Sampling
'Brute force' dynamics simulations using an Andersonlike thermostat are not efficient enough to generate a representative ensemble of trajectories that start in the single-stranded state and end in the duplex state. Thus, we resorted to using Forward Flux Sampling (FFS) to more efficiently calculate fluxes between local free-energy minima as well as sample the transition pathways. The term 'flux' from (meta)stable state A to state B has the following definition: Given an infinitely long simulation in which many transitions are observed, the flux of trajectories from A to B is where N AB is the number of times the simulation leaves A and then reaches B, τ is the total time simulated, and f A is the fraction of the total time simulated for which state A has been more recently visited than state B.
FFS requires use of an order parameter, Q, which provides a descriptive measure of the extent of the reaction between states A and B. Additionally, the order parameter must be chosen such that non-intersecting interfaces λ Q Q−1 can be drawn between consecutive values of Q. At the beginning of an FFS implementation, a brute force simulation is run starting from states described by Q = -2, and the flux of trajectories crossing the surface λ 0 −1 is measured. The total flux of trajectory from Q = -2 to another free-energy minimum Q = Q max can be calculated as the flux of trajectories crossing λ 0 −1 multiplied by the probability that trajectories subsequently reach Q = Q max , all before returning to Q = -2. The probability can be factorized as (S10) The first term in the product on the right-hand size of Eq. S10 is calculated by loading random configurations that have just crossed λ 0 −1 , which are used to estimate P (λ 1 0 |λ 0 −1 ). The process is then iterated for successive interfaces, and the flux as well as the trajectories that successfully reach Q max from the distribution of pathways can be measured.
We estimated the random error in the FFS simulations in the following way. In Table 1 in the main text we report the mean value for the hybridization rates from 5 identical and independent implementations of FFS for each system. The error reported for each system is the standard error of the mean value. In Tables SIII, SIV, and SV, we report the mean and the standard error of the mean for each individual interface. Each calculation of the flux in the three systems was repeated 240 times in total (48 simulations were used to obtain the flux in each independent calculation of the rate), while the probability of crossing interface λ Q Q−1 was computed from the 5 independent calculations of the rate.

SIII. SIMULATION PROTOCOLS
In this section we discuss the implementation of the algorithms of Section SII for both single-stranded and duplex systems. As mentioned in Section SII B 1, we simulated the three duplex systems using molecular dynamics and VMMC simulations. Unless otherwise stated, the temperature in a simulation was taken to be T = 20 • C, which is the same temperature used by Gao et al. in the experiments. Additionally, for simulations of the duplex systems, we used a simulation box with a volume of 3.96 × 10 −23 m 3 which corresponds to a concentration of 42 µm, and is 21 times larger than the experimental concentrations of 2 µm used for each system. We also used two types of order parameter in the simulations that we combined into multi-dimensional order parameters, which are discussed individually for each simulation in the sections to follow. Specifically, a 'distance order parameter' measures the minimum distance between hydrogen-bonding sites over correct pairs of bases in the two strands. A 'bonds' order parameter measures the total number of base pairs, which can be specified to be intra-or inter-strand base pairs. The definition of a bonded base pair in our simulations is two bases with a hydrogen bonding energy below 0.596 kcal mol −1 . This value for the selected cutoff corresponds to about 15% of typical hydrogen-bond energy.
A. Thermodynamics

Single-strand Thermodynamics
Melting properties for "monomers" (secondary structure of isolated strands) can be calculated from Φ, the ratio of bound to unbound states in a simulation of a single strand. For self-interacting strands, the fraction of folded states in a hypothetical bulk system is concentration-independent and can be inferred from The melting temperature is taken to be the point where f mon bulk = 1/2. We determined the approximate location of the melting temperatures of the hairpins to ensure that they were stable at T = 20 • C. OxDNA has been shown to reproduce the dependence of hairpin melting temperature on stem-length and loop-length [1]. We ran 10 independent VMMC simulations for 6.3 × 10 10 steps at temperatures of T = 20 • C for both the P 0 and T 0 strands, 1.1 × 10 11 and 8.8 × 10 10 steps at T = 45 • C for the P 3 and T 3 strands respectively, and 1.1 × 10 11 steps at T = 60 • C for the P 4 and T 4 strands. Defining melting temperatures of hairpins is complicated because strands may exist in multiple stable structures. We were interested in the point where the strands did not have significantly stable intra-strand base pairs, so we counted all states with at least one intra-strand base pair as 'bound' and states with no intra-strand base pairs as 'unbound', then calculated the yields from Eq. S11. The yield curves for P and T strands are obtained by single histogram reweighting and are shown in Fig. S3, and in Fig. 2 in the main text for just the P strands.
In addition to the melting curves, we also determined free-energy landscapes at T = 20 • C from VMMC simulations with umbrella sampling for the P 0 , P 3 , and P 4 . For P 0 strands we used an order parameter that kept track of any intra-strand base pair. For the P 3 and P 4 strands, we used a two-dimensional order parameter where the two coordinates describe (1) the intended base pairs according to Gao et al. [14] ('correct' base pairs) of the structures, and (2) all possible other intra-strand base pairs (total base pairs minus correct base pairs). For a typical strand there was ∼80 possible different intrastrand base pairs. The free energies for the strands were calculated from cumulative distributions from 10 parallel runs with 5 × 10 10 , 5.5 × 10 10 , and 2.5 × 10 10 steps for P 0 , P 3 and P 4 , respectively. The results for the free energies and yields of the single strands are shown in Fig. 3 in the main article and in Fig. S3, respectively, and are discussed in Section SIV A 1.

Duplex Thermodynamics
We first computed the melting temperatures of the three duplex systems. For structures consisting of two molecules care must be taken in extrapolating from a simulation of two strands to a bulk solution with many more strands, because fluctuations in local concentrations play an important role. If Φ is the ratio of bound to unbound states in a simulation of two molecules, the yield of a non-self-complementary duplex in a bulk solution (with the same average concentration of reactants) is given in Ref. 15 as The melting temperature occurs when f dim bulk = 1/2, which corresponds to a simulation yield of φ = 2. To compare simulations of single duplexes with experimental data, as in Table SII for the melting temperatures of the 3 systems, Φ was measured in simulations and then scaled to the experimental concentration (Φ is proportional to the concentration, so scales from a concentration c 1 to another concentration c 2 by the factor c 2 /c 1 ). The bulk yield was then calculated by using Eq. S12. Note that this approximation only works if the systems are essentially ideal -the accuracy of this approximation has been previously established for oxDNA under similar conditions [15,16]. We ran 10 VMMC simulations with umbrella sampling for each system using an order parameter that measured the number of inter-strand bonds between the two strands. The simulations for the three systems were carried out at T = 77.5 • C, which is near the melting temperature of each system. For P 0 T 0 , P 3 T 3 , and P 4 T 4 , each of 10 simulations ran for 1.6 × 10 10 steps, 2.9 × 10 10 steps, and 1.8 × 10 10 steps, respectively. The results for the duplex yields are shown in Fig. S4 and discussed in Section SIV A 2.
Next, the computations of the relative free energies of the P 0 T 0 , P 3 T 3 , and P 4 T 4 systems at T = 20 • C were carried out using VMMC moves along with umbrella sampling with a multi-dimensional order parameter that measures (1) the number of intra-strand base pairs in the P strand, (2) the number of intra-strand base pairs in the T strand, and (3) the number of inter-strand base pairs between P and T strands. In the order parameter, any complementary bond is taken into account. This means we include all secondary structural base pairs in the single strands. We ran 10 simulations for 2.3 × 10 11 steps for the P 0 T 0 system and for 1.1 × 10 11 steps for the P 3 T 3 system. The results are shown in Fig. 5 in the main article for P 0 T 0 , P 3 T 3 , and P 4 T 4 , and discussed in Section SIV A 2.
The same order parameter used for P 0 T 0 and P 3 T 3 was initially used for P 4 T 4 . Ten simulations utilizing umbrella sampling each ran for 1.1 × 10 11 VMMC steps. These results are shown in Fig. 2(c) in the main text. During the course of the simulations, we noticed a difficulty in sampling the 'pseudoknotted' intermediate states (in which the strands were bound by both tails and loops of the hairpins, illustrated in Fig. S6) that were observed in kinetic simulations; the order parameter was not efficient in driving their formation. We sampled these intermediate states in separate simulations with a distinct multi-dimensional order-parameter depending on (1) only the intended 4-stem hairpins in the P strand and (2) in the T strand, (3) the number of base pairs between the loops of the hairpins, (4) the number of base pairs between the two strands not including the loop-loop base pairs, and (5) the correctly aligned duplex base pairs between the strands. We sampled only those states in which at least one intended hairpin base pair was present in each strand, and also only states where coordinate (4) had at least one base pair formed. These simulations allow an estimate of the free energy of the pseudoknotted state compared to the tail-bound state. We ran a set of 10 simulations using this order parameter for 2 × 10 11 VMMC steps each. These results are shown in Fig. S5 and are discussed in Section SIV A 2.
B. Kinetics

FFS Simulation Details
In this section we discuss the implementation of the FFS algorithm, discussed in Section SII. As mentioned in Section SII B 1, we simulated the three duplex systems using molecular dynamics at the experimental temperature of T = 20 • C. Additionally, in all kinetics simulations we used a simulation box with a volume of 3.96 × 10 −23 m 3 which corresponds to a concentration of 42 µm that is 21 times larger than the experimental concentration, as was noted in Section SIII. We also use the same definition of a bonded base pair that was discussed in Section SIII.

Order Parameter Used in FFS Simulations
The order parameter used in simulations is detailed in Table SI. Specifically, we use a combination of distance and bond criteria as outlined in Section SIII. Distance criteria are used to define states Q = −2 → 2, and bonding criteria for states Q = 2 → 7. For the Q = 2, 3 states, we allowed the bond criteria to track any inter-strand bond between the two strands, which allowed us to monitor the number of non-intended inter-strand base pairs (i.e. mis-aligned base pairs) and also the number of correctly aligned inter-strand base pairs that have formed during the initial association events. The bond criteria for states Q = 4 − 7 track only correctly aligned inter-strand base pairs.

Initial Equilibration of Single-strand States
Before implementing FFS, we performed lengthy equilibration simulations to ensure that the single strands were initialized in thermodynamically representative states. Here we describe the procedure used to select these states. In Section SIV A 1 the relative free energies for each single strand were computed using VMMC with umbrella sampling. We performed similar simulations except that both P and T strands were simulated in the same box corresponding to a concentration of 42 µm. A 3dimensional order parameter was used that measures (1) the minimum distance between any pairs of nucleotides that are intended to be base pairs in the final duplex, (2) the number of intra-strand base pairs in the P strand, and (3) the number of intra-strand base pairs in the T strand. The strands were prevented from coming within 5.1 nm of each other, as measured by coordinate (1). We ran 10 simulations for 5 × 10 9 VMMC steps and saved configurations every 5 × 10 6 VMMC steps, which ensured that any two saved configurations were energetically decorrelated from each other. In total we collected a set containing 10000 configurations. For each state described by the 3-dimensional order parameter there is an umbrella bias w(Q). We randomly selected a configuration from the set and saved it to be used in FFS simulations if w(min)/w(Q) ≤ R, where R is a random number selected within the range 0 < R ≤ 1 and w(min) was the smallest biasing weight applied in the simulation. This step was repeated until 200 configurations were obtained for each of 5 independent FFS simulations. At the start of each flux generation simulation an initial configuration from the saved set of 200 was selected at random and set to be the starting configuration.

Second-order Kinetics Approximation
An important question is whether or not the approximation of instantaneous reactions (i.e. second-order kinetics) is valid for oxDNA. Such an approximation is reasonable if the time taken from first interaction to full duplex formation or separation of strands is small compared to the diffusional time scale governing the first contact between strands. We did observe simulations spending significant computational times in states where the partially hybridized strands had formed kissinghairpins. Theoretically, FFS should account for intermediates states with long lifetimes during flux generation. However, as the formation of such states is rare during flux generation, the sampling is poor. As an alternative to the brute-force approach, we assumed second-order kinetics, reducing the sampling challenge during flux generation, and then checked the accuracy of the assumption from the resultant data.
During flux generation, we therefore restarted (from Q = −2) trajectories that reached Q = 3, the first state in which a bond is present between strands. Consequently, any time spent in configurations with bonds between strands was not measured in our simulations. Technically our overall FFS protocol measured the flux from Q = −2 to Q = 2, and the subsequent probability of reaching Q = Q max before returning to Q = −2. This approximates the flux from Q = −2 to Q = Q max provided that the time spent in the intermediate states is small. In this limit, the measured flux is also proportional to the second-order rate constant.
To justify this assumption, we measured the time taken for shooting trajectories launched from intermediate values of Q = Q to reach Q = Q + 1 or Q = −2. We could therefore determine the typical time taken for a configuration starting in the state Q to either rearrange and proceed to a full duplex (Q = Q max ), or to dissociate (taken to be when the system reaches the state Q = −2), for comparison with the diffusional time scale of first contact. These results are shown and discussed in Section SIV B 1. We find that the second-order approximation is reasonable at the concentration used in the simulations (and would be even better at the experimental concentration, 21 times lower) and therefore the relative fluxes estimated by our approach are decent predictions for the relative rate constants in oxDNA.

SIV. RESULTS
A. Thermodynamics

Single-strand Thermodynamics
Gao et al. used the mFold software [17] to design the strands listed in Table 1 in the main article, and assumed the predicted lowest energy structures to be the only important structures in their investigation. mFold uses the nearest-neighbor thermodynamic model developed by Santa Lucia [4] to analyze secondary structure. However, this model cannot yet incorporate more complex structures like pseudoknots or multiple internal loops. It cannot take into account forces that may result from the FIG. S2. Schematic representations of the significant hairpin states of the P3 and T3 strands. The (0, 3) state is the hairpin predicted by Nupack to be the most stable, the (2, 3) state is the predicted hairpin with the shorter tail folded in, and the (3, 0) state is a hairpin which appears to be slightly more stable than the predicted hairpin in oxDNA simulations, with the stem at the other end of the strand. three-dimensional structure, which can be important in some cases [18]. mFold predicts that at T = 20 • C the P 0 T 0 strands, while designed to minimize their secondary structure, had multiple possible structures with free energies close to zero relative to the hairpin-free case, showing the difficulty of eliminating hairpins completely from long strands. From the simulations, we also found that the P 0 and T 0 strands were not dominated by any particular hairpin, but did frequently have some limited secondary structure. These transient base pairs should have a limited effect on hybridizationas they can melt easily.
To determine the prevalence of stable secondary structure at room temperature, we calculated the yields of the three P strands, which were plotted in Fig. 2 in the main article. Comparable results for the T strands are shown in Fig. S3, which show almost no difference between P and T strands.
OxDNA and mFold disagree slightly about the relative stabilities of similar hairpins within the P 3 T 3 system. However, these subtleties are likely to be relatively unimportant, as both predict that hairpins are much more stable in P 3 T 3 than P 0 T 0 , and also much more stable in P 4 T 4 than in P 3 T 3 . Further, the main hairpin stem is always predicted to be of 3 and 4 base pairs in P 3 T 3 and P 4 T 4 respectively.   [4], and simulated with oxDNA.

Duplex Thermodynamics
Yield curves for bulk solutions in the region of the melting temperature are plotted in Fig. S4 for the three duplex systems, and the melting temperatures are listed in Table SII alongside the experimental values as measured by UV absorbance spectroscopy by Gao et al. [14], and the predicted values that were calculated by the Santa Lucia model [4]. All three methods are in agreement as to the order of the melting temperatures, although oxDNA appears to underestimate the true value by around 3 • C, corresponding to an error of <1%. This is unlikely to be of significance. What is important for this investigation is that all duplexes melt at temperatures well above room temperature, the order of stability is reproduced, and the differences between the curves is small.
We also computed separately the free energies of the intermediate states in the P 4 T 4 system relative to the state (5,8), which is shown in Fig. S5, using an order parameter designed to allow the sampling of the pseudoknotted intermediate with inter-strand base pairs between stem and loop. An example of a pseudoknotted intermediate containing 15 inter-strand base pairs and 8 intra-strand base pairs is shown in Fig. S6. The plot in the main text clearly shows a metastable intermediate at (5,8), corresponding to two hairpins bound by their tails. This is also visible in Fig. S5, as is a second local minimum at (13,8) corresponding to the pseudoknotted state. This plot indicates that the pseudoknotted state is less stable than the tail-only state. Kinetic results in Section SIV B show that interchange between the two minima is reasonably fast.

B. Kinetics
The hybridization rate constants, k + , the melting rate constants, k − , and the equilibrium constants K eq , for each system are listed in Table 2 in the main article. The melting rates, k − , were not computed using FFS, but rather by using Eq. 1 in the main article combined with our calculations of k + and exp(∆G 0 /k B T ) from the freeenergy calculations in Section SIV A 2. The cumulative statistics of the FFS simulations for P 0 T 0 , P 3 T 3 , and P 4 T 4 systems are presented in Tables SIII, SIV and SV, respectively.

Considerations of the Kinetic Intermediate States
The formation of the P 4 T 4 duplex is suggested by Gao et al. to have two different kinetic regimes which cannot be fitted to a simple two-state model. They propose a 'fast' regime where the tails of the hairpins bond (or perhaps their loops kiss); and a 'slow' regime where the hairpin stems are displaced by inter-strand base pairs as the strands zip up. The fast regime has a rate constant smaller by a factor of 6, and the slow regime has a rate constant smaller by a factor of 25, than the P 0 T 0 duplex. In order to see non-second-order behavior, the intermediate state needs to be long-lived so that both dissociation and completion of the reaction are slow relative to the association rate.
We find the rearranging time for metastable states to proceed to a full duplex for the majority of simu-lations to be typically less than ∼1 × 10 −7 seconds for P 0 T 0 (Fig. S7(a)), and ∼1 × 10 −5 seconds for P 3 T 3 and P 4 T 4 (Figs. S7(b) and (c)), while the longest rearrangement times that lead to dissociation events took less than ∼1 × 10 −5 for P 3 T 3 (Fig. S8(b)) and ∼4 × 10 −4 seconds for P 4 T 4 (Fig. S8(b)). Comparing the longest rearrangement and disassociation times with the diffusion time, we find for all duplex systems studied that the most longlived metastable states (in P 3 T 3 and P 4 T 4 systems) have a lifetime slightly smaller than the rate at which they are produced. Thus reactions are second order to a reasonable approximation in oxDNA at concentrations of 42 µm (justifying our simulation procedure). Further, our model, consistent with Nupack, suggests that regardless of completion rate, the 6-base pair toeholds (i.e. the tails of the intended 4-stem hairpins) are not stable enough to give rise to two kinetic regimes at the experimental concentration used. The strands just fall off too quickly, which is unsurprising given the known physics of DNA. We find that the kissing hairpin loops are even less stable, as is the pseudoknotted configuration formed when both the loops and tails bind. Our simulations do not support the claim of two different kinetic regimes for the formation of the P 4 T 4 duplex. Therefore, non-second-order behavior in Gao's experiment, if real, must be due to some unknown aspect of DNA thermodynamics that is not incorporated into the oxDNA model.

Additional Results for Hybridization Pathways
In the top panel of Figs. S9 and S10 we plot the frequency that a given base in the P strand is involved in base pairing for configurations that have crossed interfaces λ Q Q−1 , for Q = 3, . . . , 6, respectively. In the bottom panel of Figs. S9 and S10 we plot the probability that the base pairs between the strands led to a fully-bound duplex state, given that a configuration crossed the λ Q Q−1 interface, for Q = 3, . . . , 6, respectively. For example, 2 → 3 refers to the configuration having crossed interface λ 3 2 that started from λ 2 1 . The quantity plotted on the y-axis is actually a probability density. Note that the times for configurations that crossed λ 7 6 coming from λ 6 5 are not plotted as they were found to be neglible.  Fig. S7, the labels in each legend indicate a configuration having crossed a particular interface λ Q Q−1 . For example, 6 → −2 refers to a configuration having intially crossing λ 6 5 , but crossed λ −1 −2 before crossing λ 7 6 . As in Fig. S7, the quantity plotted on the y-axis is actually a probability density. Note that the times for configurations that crossed λ −1 −2 coming from λ 6 5 are not plotted as they were found to be neglible.