Probing the Mechanical Properties of DNA Nanostructures with Metadynamics.

Molecular dynamics simulations are often used to provide feedback in the design workflow of DNA nanostructures. However, even with coarse-grained models, the convergence of distributions from unbiased simulation is slow, limiting applications to equilibrium structural properties. Given the increasing interest in dynamic, reconfigurable, and deformable devices, methods that enable efficient quantification of large ranges of motion, conformational transitions, and mechanical deformation are critically needed. Metadynamics is an automated biasing technique that enables the rapid acquisition of molecular conformational distributions by flattening free energy landscapes. Here we leveraged this approach to sample the free energy landscapes of DNA nanostructures whose unbiased dynamics are nonergodic, including bistable Holliday junctions and part of a bistable DNA origami structure. Taking a DNA origami-compliant joint as a case study, we further demonstrate that metadynamics can predict the mechanical response of a full DNA origami device to an applied force, showing good agreement with experiments. Our results exemplify the efficient computation of free energy landscapes and force response in DNA nanodevices, which could be applied for rapid feedback in iterative design workflows and generally facilitate the integration of simulation and experiments. Metadynamics will be particularly useful to guide the design of dynamic devices for nanorobotics, biosensing, or nanomanufacturing applications.

Introduction

In structural DNA nanotechnology, a collection of DNA sequences is chosen to form a desired structure via molecular self-assembly.[1,2] Such DNA constructs often have a single well-defined free energy minimum, corresponding to geometries like ribbons,[3] tiles,[4] square, or honeycomb arrangements of helices,[5,6] or brick-like voxel arrays of short DNA oligonucleotides.[7] These unimodal structures (i.e., having one primary configuration in space) have been translated to applications where structural rigidity is important—in fiducials for super-resolution microscopy,[8] as scaffolds to visualize biomolecular processes,[9] or as nanopores for single-molecule detection.[10]

With the DNA origami technique, a bacteriophage genome and synthetic oligonucleotides coassemble to create near-arbitrary shapes.[4−6,11] DNA origami has emerged as a dominant approach in nanoscale structural design, and it unlocked the manufacture of nanostructures programmed to perform complex motion,[12−14] for example, hinges,[15] pistons,[16] interlocked axles and sliders,[16,17] and rotors.[18−20] These deformable elements have formed the basis of stimuli-responsive materials,[21] sensors,[22] single-molecule probes,[19,23,24] drug delivery vectors,[25] and nanoreactors.[26] The motion of origami nanomachines can be constrained to occur along given axes, and configurational distributions can feature multiple stable states separated by energy barriers.[12,27]

Molecular modeling has become a key element in the design workflow of DNA nanostructures, with the two most common approaches being finite-element modeling and Molecular Dynamics (MD). Finite-element frameworks, such as Cando[28] and SNUPI,[29] describe DNA helices as elastic rods and apply continuum mechanics to predict the equilibrium structure and its deformation modes. The latter are, however, only accurate in describing small deformations, and they become poor approximations when the structure deforms significantly or has multiple stable states. Additionally, these continuum approaches lack the resolution to describe molecular processes such as formation and dissociation of base pairing and stacking bonds, which may be critical for the behavior of dynamic devices.

Conversely, MD infers mechanical properties by an explicit simulation of the system’s Newtonian dynamics. Atomistic simulations of DNA nanostructures may take weeks to complete,[30] motivating the development of coarse-grained models such as the multiresolution DNA (MrDNA) framework[31] and oxDNA.[32]

Thanks to its ability to accurately represent nucleotide stacking and base pairing, the oxDNA force field[32,33] has succeeded in replicating various phenomena, including kinking in duplexes[34,35] and force-induced unravelling of origami,[36] and has been applied to predict conformational distributions of origami mechanical elements.[37,38] As a result, oxDNA is now frequently used as part of iterative nanostructure design workflows.[13,39]

However, even coarse-grained simulations can be impractically slow,[38] and without ad hoc biasing techniques they can only sample configurations with free energy within a few kBT away from the minima. Additionally, trajectories can become trapped in local minima, hindering complete sampling. As a result, coarse-grained simulations are often performed merely to check for mechanical strain or undesired deformations, instead of quantitatively assessing the range of motion or the forces required for actuation.

Various biasing techniques can be used to flatten free energy landscapes and accelerate sampling.[40] These approaches use fictitious forces along collective variables, which are low-dimensional representations of conformational states. One such method, previously integrated with oxDNA,[41] combines steered MD with the use of the Jarzynski equality[42] to reconstruct free energy landscapes along a one-dimensional (1D) reaction coordinate. This method is, however, unsuitable for acquiring multidimensional landscapes, and its estimates are dominated by unlikely low-work trajectories, resulting in difficult-to-assess uncertainties.[43]

Shi et al.[44] and, more recently, Wong et al.[45] have demonstrated that the integration of umbrella sampling with oxDNA can enable the exploration of 1D and two-dimensional (2D) free energy landscapes associated with the deformation of origami, while this technique had been previously applied to exploring deformations in smaller nanostructures, including duplex bending[35] and junction flexibility.[46] Umbrella sampling relies on defining multiple (partially) overlapping windows across the space of the relevant collective variables, in order to limit the scale of the free energy features that the system needs to thermally explore. A full free energy surface is then reconstructed by stitching together samples from the individual windows. While successful, this approach requires a system-specific definition of the thermodynamic windows and laborious postprocessing, making it challenging for nonexperts.

Alternatively, a single biasing potential can be designed to globally counteract the free energy profile. However, that ideal bias is unknown at the outset; it must be initially set using intuition and then iteratively refined in subsequent simulations. The fast-iteration limit of refinement is an on-the-fly update, where an optimal bias is progressively learned in a single simulation rather than optimized through separate runs—this is the idea behind metadynamics (MetaD).[47,48]

In MetaD, a bias is constructed from the history of observed configurations, which discourages the revisiting of previously sampled states. This process encourages iteratively wider exploration of the state space, eventually enabling transitions over the free energy barriers separating local minima. Even for systems with a single free energy minimum, MetaD enables sampling of high free energy states, an ability that would be particularly useful to probe force-response in DNA nanomachines and mechanical sensors.[12,49−51] MetaD simulations can benefit from graphics processing unit (GPU) acceleration and a natural parallelization route through multiwalker metadynamics.[52] The well-tempered variant of MetaD[53] limits the maximum correction to the free energy landscape, preventing irreversible disassembly. Finally, there is no need to run simulations with multiple thermodynamic windows, as in umbrella sampling, simplifying execution and postprocessing and eliminating some concerns about hysteresis.[54]

Here we introduce an implementation of well-tempered MetaD in the oxDNA simulation framework, which offers a viable route for the rapid assessment of conformation free energy landscapes in DNA nanotechnology. To demonstrate the validity of the technique we applied it to four case studies where conventional MD would be unable to probe the relevant landscapes: (i) the compression-induced buckling in duplex DNA, (ii) conformer transitions in bistable Holliday junctions[55] and (iii) switchable tiles,[27] and (iv) force response in an origami-compliant joint, where conformation is prescribed by balancing competing forces.[49] For systems (ii) and (iv) we compared simulation outcomes with experimental observations, finding quantitative agreement. Overall, we demonstrated that MetaD, as applied to oxDNA, can effectively sample transitions between multistable systems and facilitate the computational characterization of highly deformable designs, all in an automated fashion that requires limited system-specific user input. This tool could therefore be highly valuable in computer-assisted design and assessment pipelines for reconfigurable DNA nanostructures.

Results and Discussion

Principles of Metadynamics

Here, we give a brief overview of the principles and implementation of MetaD. A complete theoretical description can be found in Bussi et al.[48] The objective of MetaD is to map a free energy landscape from a molecular simulation. As landscapes typically have high dimensionality, for human interpretation the free energy is projected onto a set of lower-dimensional coordinates or collective variables (s⃗), defined as functions of the coordinates of the simulated system (q⃗).

In principle, long trajectories sampled from Monte Carlo (MC) or MD can be used to infer free energy landscapes from state-occupancy histograms. The projection of the free energy onto a discretized coordinate s⃗0 can then be estimated as ΔG(s⃗0) ≃ −kBT logN(s⃗0) + c, where c is an immaterial constant, kB is the Boltzmann constant, T is the temperature, and N(s⃗0) is the number of samples in the histogram bin centered at s⃗0.[56] However, convergence of this unbiased approach is practically unfeasible for many macromolecular and DNA nanosystems owing to the presence of thermally inaccessible configurations that frequently separate multiple metastable minima.

MetaD generates a history-dependent bias that progressively flattens the free energy landscape, thus rendering high free energy regions accessible and enabling efficient sampling.

A MetaD simulation proceeds as follows. The system is initialized and simulated (with either MD or MC algorithms) using a potential defined as Ut(q⃗) = U(q⃗) + Bt(s⃗(q⃗)), where U(q⃗) is the unbiased potential and Bt(s⃗(q⃗)) is the time-dependent bias. The index t indicates the number of MetaD iterations performed, each iteration consisting of τ (MD orMC) time steps. The bias is initialized as B(s⃗(q⃗)) = 0 and is updated after each iteration to counteract the projection of the free energy onto s⃗. To calculate the updated bias, the instantaneous value of s⃗ is evaluated, termed s⃗t. The bias is then updated through the addition of a Gaussian potential centered at s⃗t, which discourages the system from revisiting its current state.

In eq , σ is the width of the deposited Gaussian, while the parameter w controls the rate at which the free energy wells are filled. In the earliest version of MetaD, also known as direct MetaD, w was set to a constant value,[47] resulting in a Bt that oscillates rather than converging.[48] Alternatively, convergence of Bt can be guaranteed by reducing w in areas that are already strongly biased, an approach known as well-tempered metadynamics,[53] which we adopt throughout this work. In well-tempered MetaD, the time-dependent amplitude of the Gaussian, wt, is given by

In eq , ΔT is an additional hyperparameter with units of temperature, which controls the strength of tempering. High values of ΔT correspond to weak tempering, where forces are allowed to accumulate, with ΔT → ∞ approaching conventional MetaD (constant w). Conversely, low values of ΔT correspond to systems that quickly taper their bias, with the ΔT → 0 limit corresponding to unbiased sampling. The value of A controls the initial bias-height increment. ΔT and A are set at the start of the simulation, alongside the other parameters (σ and τ) and the collective variables. With well-tempererd MetaD, at long times, the value of Bt(s⃗) provably converges to a fraction of the projection of the free energy onto the collective variable (up to an immaterial constant, c).[57]An estimate of ΔG(s⃗) can therefore be acquired from the converged bias.[58] Additionally, this equation illustrates the physical interpretation of ΔT. After convergence, the residual (i.e., uncorrected) free energy felt by the system is , implying that T + ΔT can be interpreted as the effective temperature experienced along a collective variable.[48]

While eq enables one to estimate ΔG, a preferred route is that of directly extracting the sought free energy from configuration histograms of simulation runs biased with the asymptotic Bt. This approach will be used to derive free energy landscapes in the remainder of this article, unless specified otherwise.

Supporting Information, Note 1 and Figure S1, demonstrate the implementation of MetaD with a basic one-dimensional example, while, in the reminder of this paper, we illustrate its applications to mapping deformation free energy landscapes for increasingly complex DNA nanosystems, simulated with MD and the coarse-grained oxDNA force field.

Information on the implementation of MetaD in oxDNA and specific simulation details for all case studies can be found in the Methods section and Tables S1 and S2.

While the convergence of (well-tempered) MetaD is very robust, the free parameters σ, A, and τ, alongside system-dependent features such as physical size, intrinsic diffusion times, and collective-variable dimensionality, have been shown to influence errors in free energy estimates and convergence time scales.[59,60] In the Methods we discuss these factors and other practical considerations that guided our parameter choice.

Bending and Buckling Free Energy of a DNA Duplex

In this section we demonstrate the application of MetaD to coarse-grained oxDNA simulations using a simple case study: the response of double-stranded (ds) DNA under strong bending. A similarly simple application is discussed in ref (61), which explores bubble formation in a basic bead-and-spring model of a DNA duplex.

dsDNA is often thought of as a Worm-Like Chain (WLC)—an elastic beam whose bending energy is quadratic in local curvature, much like a macroscopic beam. If the ends of such a duplex are compressed together, then the WLC model predicts that the curvature will increase everywhere. However, experimental evidence indicates that, under a sufficient compressive load, a short dsDNA duplex will not bend continuously. Instead it will buckle, and in this buckled state there will be a single point of high curvature—a kink.[62] Experimental observations of force-induced kinking have been identified for a DNA-based molecular vice in fluorimetry experiments,[63] in the vulnerability of dsDNA minicircles to single-stranded (ss) DNA-specific enzymatic degradation,[64] and also via atomic force microscopy (AFM) of said minicircles.[65] Similarly, kink formation under conditions of end-to-end compression has also been observed in an atomistic simulation[66] and with the oxDNA force field.[34,35] Both atomistic and coarse-grained simulations indicate that the origin of kinking is a local break in the continuity of coaxial stacking in the helix[34,35,66] and may be also associated with the loss of a Watson–Crick bond. Here we have used sampling of DNA kinking as a simple test application of MetaD in oxDNA.

We wish to apply MetaD to calculate how free energy varies with the end-to-end distance of a short duplex DNA, which, due to the complex buckling transition, is impossible to calculate analytically. Figure a shows snapshots of the unbuckled (left, A), and buckled (right, B) configurations of the duplex. The distance x between the centers of mass of the two collections of six cyan beads was used as a collective variable onto which the free energy is projected and the MetaD bias Bt applied. In Figure b, the time evolution of Bt(x) was plotted, along with reference free energy ΔG(x)—the true energetic cost to bend the duplex. The uncorrected potential Bt(x) + ΔG(x), that is, the residual potential felt by the system, progressively flattens as Bt(x) evolves according to eqs and 2 to counteract ΔG(x) (Figure c).

Figure 1

MetaD enables automated sampling of dsDNA buckling. (a) Snapshots of unbuckled (left, A) and buckled (right, B) configurations of a DNA duplex from a MetaD simulation. The buckled state features disrupted stacking roughly in the center of the duplex. The distance x between the centers of mass of the two collections of six cyan beads was used as the collective variable. (b) The time dependence of the bias Bt for a system with ΔT = 8T. Also plotted is ΔG, the unbiased potential experienced by the system (black dashed line). Bt is initially flat, and then it builds up according the history of visited configurations (eq ). (c) The simulation experiences a potential equal to ΔG + Bt – the uncorrected potential. As illustrated, the initial uncorrected potential is sharply varying, but then it progressively flattens as the bias grows, enabling access to a wider x-range. Different colors mark different numbers of MD time steps, as indicated in the (b, c) legend. (d) Implied free energy from the asymptotic Bt, for varying ΔT (eq ). ΔG(x) is plotted as a black dashed line. (e) Trajectories of the collective variable x for ΔT = 0 (ordinary MD), ΔT = 8T, and ΔT = 32T. (f) 2D free-energy landscape acquired from biased MD simulation. The y-axis indicates Ufourth least stack, which rises to 0 only if at least four nonterminal nucleotides lack stacks—i.e., a buckled state. Locations marked as A and B correspond to the snapshots in (a).

In the examples given in Figure b,c, the tempering parameter ΔT (eq ) was set to 8T, so that the bias converges to (eq ); that is, the asymptotic uncorrected potential is of the true value. Figure d shows free energy implied according to eq for different values of ΔT, compared with the reference free energy (see Figure S2 for proof of convergence of the biases). As expected, larger ΔT values produce accurate estimates of ΔG away from the minimum.

In Figure e, the time-varying values of x are given for three different values of ΔT. Under conventional MD (ΔT = 0, blue), only the unbuckled state is sampled. When metadynamics is turned on (ΔT = 8T or 32T, yellow and green, respectively), the bias repels the system from previously visited configurations, resulting in a wider exploration in the unbuckled free energy minimum. From ∼5 × 107 MD time steps, both biased systems begin exploring the buckled state at smaller x-values, only briefly for ΔT = 8T and more persistently for ΔT = 32T. The latter simulation then experiences frequent transitions between buckled and unbuckled states.

In the one-dimensional free energy profile projected along x, the configurations corresponding to the buckled and unbuckled states do not appear separated by a free-energy barrier. However, such a potential barrier exists and can be visualized along alternative coordinates, as shown with the two-dimensional free energy landscape in Figure f. Here, we introduce a second collective variable, Ufourth least stack, defined as the value of the fourth weakest stacking interaction, which we expect to increase as the duplex buckles and a kink forms. Indeed, in Figure f we observe two distinct states: a broad minimum at large x and finite (negative) Ufourth least stack, associated with the unbuckled duplex, and a second minimum centered at smaller x and with Ufourth least stack = 0, corresponding to the buckled duplex. Transitions between the two minima are not effortless even with ΔT = 32T, but good sampling is possible with many replicas that are run simultaneously, sharing and contributing to the same bias (see Methods).

ΔT can be used to control which parts of the free energy landscape should be explored and the trade-off between sampling a large region of collective variable sparsely or a small region well. It can also be used to eliminate the sampling of states that may be undesirable. For example, a low value of ΔT could be used to prevent sampling of a kink formation if the objective were to identify only bending close to the free energy minimum.

Two-Dimensional Isomerization Landscape of Bistable Motifs

While in our first case study a single collective variable was sufficient to bias the simulation and extract the sought information, it is often the case for (relatively) more complex DNA architectures that multidimentional free energy landscapes need to be explored. To this end, Holliday junction isomerization provides a useful case study. The immobile Holliday junction was the first nontrivial DNA motif to be intentionally constructed,[67] and it consists of four helices joined at a central four-way junction. Its configuration in the presence of divalent or high concentrations of monovalent cations is that of two quasi-continuous helices joined at a strand crossover location. This is referred to as the stacked-X configuration[55] and is shown in Figure a (left, right). In the absence of such cations, the construct acquires an unstacked planar configuration, where each of the four arms can move flexibly about the central junction (Figure a, center).[55]

Figure 2

MetaD enables sampling of the isomerization free energy landscape of bistable Holliday junctions. (a) A Holliday junction consists of two quasicontinuous duplexes joined by a crossover, as illustrated in the snapshot. There are two dominant conformers, one where the gray and red strands are fully stacked (left, isoI) and another where the orange and blue strands are fully stacked (right, isoII). An unstacked structure is believed to be the intermediate (center). (b) For MetaD simulations we used a two-dimensional collective variable, (x1, x2), where x1 is the distance between the centers of mass of the orange and blue sets of beads, while x2 is the distance between centers of mass of the red and gray beads. (c) A 5 × 107 time-step trajectory simulated under unbiased MD (left) and one of the same duration collected with MetaD (center), both overlaid with the (x1, x2) free energy profile. Simulating over 10 times this period in MetaD results in many transitions between conformers, enabling accurate sampling of the free energy landscape (right). Dots of different colors indicate sampled configurations. (d) Free energy surfaces corresponding to (i) the sequence-averaged model at 500 mM ionic strength and the sequence-specific model[68] at (ii) 500 and (iii) 100 mM ionic strength. See Figure S3 for data on the sequence-specific model at 200 mM ionic strength. In each case, there are two minima, corresponding to the two stacked-X conformers, and a saddle-point region associated with the intermediate. (e) (top) Probabilities for the isoI and isoII states for the four studied systems (bars), compared with experimental values (dashed (isoI) and solid (isoII) black lines).[69] (bottom) Free energies of the isoI, isoII (color-coded as in (d)) and intermediate states (black). The free energy of the intermediate falls by ∼7 kBT between the systems with 500 and 100 mM ionic strength, consistently the experimentally observed phenomenon of faster isomerization at low salt concentrations. Error bars are the standard error based on six replicas (too small to see for isoI and isoII).

Stacked-X Holliday junctions can exist in two conformers, distinguished based on which of the four helices are stacked at the junction (Figure a, left and right). These conformers, previously referred to as isoI and isoII,[55] are structurally equivalent if the base sequence is ignored, while asymmetry of base pairs at the junction results in one conformer being favored. In the presence of MgCl2, each conformer is long-lived—single-molecule Förster Resonance Energy Transfer (FRET) experiments indicate lifetimes of milliseconds to seconds.[69] Consequently, sampling transitions between the two conformers is intractable with typical molecular-simulation approaches. For alternative representations of DNA, transition sampling has required running the simulations at a vastly increased temperature[70] or using a coarse-grained force field that overestimates the stability of the transition state.[71] The properties of the oxDNA representation of a Holliday junction have been explored previously.[46] However, transitions and sequence-dependent conformer probability have remained unexplored due to the nonergodicity of this system under conventional MD sampling. Here we show that, using a two-dimensional reaction coordinate, MetaD can successfully sample conformer transitions and determine the relative conformer stability.

The particular structure investigated here is similar to the J3 junction, previously characterized experimentally,[69] with the only difference being that the dsDNA “arms” have been truncated to 11 base pairs (bp) to enable faster simulations (see Table S3 for sequences). To favor the formation of an unstacked intermediate state, thus enabling transitions between conformers, we used a two-dimensional collective variable, corresponding to the two diagonal distances across the Holliday junction (x1 and x2 in Figure b). In the stacked-X state, one of these distances takes a high value, corresponding to the width of the junction, while the other takes a low value, corresponding approximately to the axial rise of two base pairs. Meanwhile, the planar transition state corresponds to high values of both collective variables.

To demonstrate the enhanced sampling made possible in MetaD versus conventional MD, in Figure c we plotted small sections of trajectories for both techniques. While a trajectory simulated under MD remains stuck in a single minimum (Figure c, left), using MetaD it is able to escape and sample several transitions (Figure c, center and right).

Free energies projected onto the collective variables are plotted in Figure d, as acquired from MD using an asymptotic bias from MetaD. Illustrations of the similarity between the converged MetaD bias and the free energy from biased MD simulation are given in Figure S3. Two different constructs were tested, one which ignores the base identity by using sequence-averaged parameters and one utilizing the sequence-dependent force field.[68] The former construct was simulated at 500 mM ionic strength, while the latter was simulated at three different ionic strengths (500, 200, and 100 mM). As expected, while the sequence-averaged calculations produce a symmetric landscape, sequence dependence results in asymmetry, with one conformer being favored over the other. Additionally, we observe that the intermediate region between the two conformer minima flattens at lower ionic strengths. This region corresponds to the unstacked intermediate, which thus appears to be favored by a reduction in salt concentration.

Experiments indicate that, in conditions of 50 mM MgCl2 (ionic strength 150 mM), the J3 junction will display the isoII conformer 77.4% of the time.[69] In Figure e we show the simulated conformer probability for the four studied systems. For the sequence-specific model we find that conformer probability is independent of ionic strength and in quantitative agreement with experimental observations. This agreement with experimental results is intriguing, as stacking interactions in oxDNA have not been parametrized to reproduce Holliday junction conformer prevalence but, instead, the melting transitions of duplexes and hairpins based on the Santa-Lucia parameters.[32,72] The reproduction of conformer probability is further validation of the oxDNA model of stacking. Definitions of the stacked and transition states are discussed in the Methods.

Figure e (bottom) shows how the free energies of the two conformers, as well as the unstacked intermediate state, depend on ionic strength. While the values for the stacked-X configurations remain constant, ion concentration is critical in controlling the free energy of the unstacked intermediate, as previously noted. Indeed, as ionic strength falls from 500 to 100 mM, the relative free energy of the intermediate decreases by ∼7kBT, making it ∼3 orders of magnitude more likely. This effect is due to increased electrostatic repulsion associated with a less-concentrated electrolyte; the stacked Holliday junction has a high density of negative charge and is therefore disfavored when electrostatic screening is reduced.

An association can be made between the stabilization of the intermediate at lower ionic strengths and the increase in conformer interconversion rate, defined as sum of the rates of isoI → isoII and isoII → isoI transitions.[69] The latter, as determined experimentally for a junction of slightly different sequence, rises from 20 s–1 at 2 M Na+ to 800 s–1 at 400 mM Na+. A similar increase is observed in systems with magnesium counterions if their concentration is dropped from 100 mM (interconversion rate 10 s–1) to 7 mM (interconversion rate 500 s–1). By assuming a direct proportionality between the interconversion rate and the probability of the unstacked intermediate, oxDNA would predict that reducing the ionic strength from 500 to 200 mM would result in a sevenfold increase in the isomerization rate, while reducing the ionic strength further, to 100 mM, would accelerate isomerization by a factor of 1000. However, we note that these considerations are purely qualitative and that rare-event sampling techniques, which do not create fictitious dynamics[73] are typically required to make definite claims about transition rates and paths.

As an additional example, in Supporting Information Note 2 and Figures S4–S6 we test MetaD on a second bistable unit where the transition between two conformers requires a breaking of stacking interactions. This tile has been utilized as the elementary unit of reconfigurable origami that can spatially relay information through the propagation of conformational transitions along an array of units.[27] Similar to the case of the Holliday junction, we are able to efficiently reconstruct the transition free-energy landscape utilizing both one- and two-dimensional collective variables for biasing, which would not be feasible with unbiased MD. We are also able to gather information on the transition pathway between conformers; however, this needs to be interpreted with care owing to potential artifacts introduced by the biasing potential.

Bending Free Energy of a Compliant Origami Joint

The principle of compliant mechanism design is to control mobility and mechanical properties via local thinning of material, rather than through rigid body linkages.[74] This is a popular approach when designing DNA origami with an intended pattern of motion, where the number of helices is reduced in regions of the structure where compliance is desired.[12,49] Here, we consider a DNA origami-compliant joint as a useful case study for the mechanical predictions of the MetaD approach. The joint has been previously characterized experimentally[49] and computationally with oxDNA using unbiased MD.[37] The latter study demonstrates that oxDNA can accurately capture the shape of compliant DNA structures, although it underpredicts the width of conformational distributions.[37] We simulated a truncated version of the experimentally realized joint, illustrated in Figure a, where truncation improves computational efficiency. The joint is composed of two 18-helix bundles, connected by a thinner 6-helix layer. Consequently, bending will preferentially occur in a plane, localized to the thinned layer. See Figure S7 for the caDNAno routing of the device.

Figure 3

Metadynamics enables sampling of mechanically stressed states in DNA origami. (a) A mechanically compliant DNA origami joint, truncated here from its experimental realization.[49] The cross section of the 18-helix bundle is also shown. The yellow beads were used as references in the bending angle ϕ (Figure ). The collective variable x is defined from the distance between the centers of mass of the top collection of six cyan particles and the six at the bottom. Individual staples and scaffold, whose routing is depicted in Figure S7, are color-coded in the right-hand-side image. (b) MD simulation (top) only samples around the free energy minimum of the compliant joint, yielding little information about the force required to actuate it. MetaD simulation (bottom) learns to sample a wider range of configurations. The illustrated trajectories correspond to approximately half the total time sampled in MetaD simulations. Different colors indicate parallel replicas that, for MetaD, contribute to, and experience, the same bias potential. (c) Snapshots are illustrated from MD simulation (top) and MetaD simulation (bottom). (d) Time evolution of the MetaD bias Bt (solid lines). Letters refer to biases at simulation times corresponding to those illustrated in (b). The long time limit bias (final) is also shown. The true free energy ΔG is shown as a dashed gray line with black 1σ error bars (often too small to see). (e) The sum of ΔG and Bt—the uncorrected potential—is plotted for different simulation times. (f) Free energy profiles (continuous lines) as implied from Bt are plotted alongside ΔG (gray dashed, 1 σ error bars), demonstrating convergence. The thicker red curve in the bottom-right subpanel represents the free energy profile as determined from direct sampling of unbiased trajectories in (b) (top).

Figure 4

Metadynamics allows one to predict the mechanical response of a DNA origami to an external force. (a) Definition of the angle ϕ. Straight lines defining the angle are those passing through the centers of mass of the groups of yellow beads nearest and furthest from the joint. See further details in the Methods. (b) The bending free energy profile against ϕ as estimated with MD using the converged MetaD bias. (c) Renders of the ssDNA connections across the joint, as implemented experimentally.[49] There are three short ssDNA sections (yellow) and three long sections (blue). (d) Free energy profiles of WLCs for each of the ssDNA distributions under study. (e) Predictions of mean compliant joint angle ⟨ϕ⟩ compared to those from TEM experiments.[49] (f) Predictions of the standard deviation of angle width σ(ϕ) compared to those from TEM. Data points are color-coded as for the corresponding WLC free energy curves in (b). All error bars represent the standard errors evaluated from different replicas, as discussed in the Methods; those in (e) are smaller than the symbols.

Through MetaD simulations, we can explore the bending free energy of the joint, sampling highly deformed configurations inaccessible to conventional MD. We bias the simulations using the collective variable x, defined as the average distance between the centers of mass of top and bottom collections of cyan beads, illustrated in Figure a. As demonstrated in Figure b (top), MD explores states only close to the free energy minimum, physically corresponding to an unstressed six-helix section. Snapshots corresponding to these trajectories are illustrated in Figure c (top). By contrast, MetaD initially explores the free-energy minimum and then is pushed by the bias to explore other regions of configuration space (Figure b, bottom). Trajectories in MetaD widen with time, not just because of diffusion but because the free energy landscape felt by the system is progressively flattened. Snapshots illustrated in Figure c (bottom) indicate the sampling of high free energy states that would never have been reached in unbiased MD.

The time dependency of the bias is illustrated in Figure d, with letters corresponding to the times marked in Figure b. Notice the shape similarity of the reference free energy to the final bias reached by the simulation. Similarly, Figure e shows the time evolution of the uncorrected potential, Bt(x) + ΔG(x), which progressively flattens as previously noted in Figure , while Figure f shows how the implied potential converges to the reference curve. In Figure f (bottom right) we also show the free energy as determined from a direct sampling of unbiased MD simulations (panel (b), top), which expectedly are only able to reconstruct the profile for thermally accessible configurations.

Experimental investigations of the compliant joint have relied on a bending angle, rather than a distance, to classify the deformation state of the nanomachines. For a direct comparison and, thus, to demonstrate the predictive power of the oxDNA MetaD approach, we defined the bending angle ϕ as illustrated in Figure a, closely matching the definition used in Transmission Electron Microscopy (TEM) experiments.[49] The associated bending free energy profile is plotted in Figure b—note once more how MetaD enables sampling high free energy states associated with extreme bending, with free energies reaching ∼60 kBT above the ground state.

In experiments, controlled bending of the joint has been induced through the addition of ssDNA segments, bridging the 18-helix bundles across the flexible section of the joint at the locations illustrated in Figure c. Three strands join yellow beads (each containing Nshort nucleotides), and another three, with possibly different lengths, join blue beads (each with Nlong nucleotides). The segments act as entropic springs, bending the six-helix bundle and determining the configuration (and flexibility) of the joint. The bending state can thus be controlled by changing the length and number of the springs.[49]

Besides assessing the flexibility of the unconfined joint, a useful role for simulations would be that of predicting the mean bending angle that results from a given set of ssDNA springs, so to inform experimental design. To this end, one approach would be to perform separate simulations for many possible lengths of ssDNA springs[37] and then manufacture the system whose behavior is closest to the desired outcome. While this is computationally costly, it is the only possible approach when states far from the location of minimum free energy cannot be sampled.

MetaD, instead, enables a much more efficient approach thanks to its ability to sample with a single simulation the entire distribution of angles, as we have shown. Once this free energy profile is known in the absence of any ssDNA, one can indeed analytically account for the constraints imposed by ssDNA springs. Specifically we can predict the bending angle distribution, by reweighting the distributions from biased MD to account for the energetic contribution of the springs, as described in the Methods. Each ssDNA section is modeled as a separate WLC between attachment points, and free energies from each contribute to the reweighting. Figure d illustrates the free energy contribution for each combination of long and short chains used here, as a function of extension. This strategy offers an efficient alternative to determine the bending-angle distribution of the joint for any choice of ssDNA springs, ensuring that the inverse problem of designing ssDNA sections to produce a given angle is approachable. Similarly, it offers a way to estimate the flexibility of the joint under applied force, useful if it were later used in a load-bearing application.

Figure e compares our predictions for the mean bending angle ⟨ϕ⟩ with experimental data of the corresponding systems, finding good agreement. Good correspondence is also observed between simulated and experimental standard deviation σ(ϕ). It is feasible that the small discrepancies between simulation and experiments emerge from inaccuracies in the WLC model of the springs. Indeed, such a model may be inappropriate for some of the shorter sections used here (down to 11 nucleotides), especially given that sequence has been ignored. Additionally, the use of WLC springs neglects possible stacking effects at the attachment points of the ssDNA springs with the 16 helix bundles, on either side of the joint. Nevertheless, despite small discrepancies, the automated reconstruction of accurate profiles of mean bending angle (to within 10°) confirms the applicability of this method to the rapid prediction of the mechanical and structural properties of DNA origami before manufacture.

Conclusions

Molecular simulation is essential in the design and interpretation of systems that use DNA to build mechanical structures. However, unbiased MD simulation gives little information about the mechanical response of these structures to an applied force. Additionally, for multistable systems with nonergodic dynamics, unbiased simulation may entirely miss certain states, which may be critical to the function of the construct. To address these limitations, here we have combined well-tempered metadynamics and the popular oxDNA force field, thus introducing a tool for the fast and automated reconstruction of one- and two-dimensional free energy landscapes of deformable DNA nanostructures, including a sampling of multiple minima and transition states in multistable devices.

To exemplify the utility of our metadynamics implementation in DNA nanotechnology, we have applied it to four case studies, associated with systems of different scales and conformational complexity. First, we have demonstrated automated sampling of the reversible kinking of a short DNA duplex under compression, replicating experimental and computational observations on the process.[34,35,63−66] We have then reconstructed the free energy landscape of bistable DNA systems whose dynamics would be nonergodic under conventional MD, even using coarse-grained models. In particular, we have analyzed a bistable Holliday junction exhibiting two possible conformers and found quantitative agreement between simulated and experimental conformer occupancy.[69] The obtained free energy profiles also helped to quantify the effect of ionic strength on the accessibility of the transition state, which qualitatively correlates with experimental trends in switching rates.[69] Additionally, we have reconstructed the free energy landscape of a bistable motif previously used for information relaying in DNA origami,[27] for which we have identified plausible reaction intermediates—a useful insight for integrating these units into signal transduction architectures.

To demonstrate the applicability of our oxDNA MetaD implementation to larger constructs, we have predicted the mechanical response of a compliant DNA origami joint to varying force. We have further shown how, thanks to its ability to map out thermally unaccessible conformational landscapes, MetaD provides an efficient route for the computer-assisted design of joints with prescribed equilibrium angles and stiffness, which we have benchmarked against experimental data.[49]

By combining oxDNA with metadynamics, we have enabled faster prediction of free energy profiles without compromising the detail of the underlying DNA model. The process of landscape acquisition can be fully automated, as it does not require a manual tuning of biasing weights and uses a single thermodynamic window, contrary to umbrella sampling, and it can therefore be accessed by users lacking advanced computational expertise. Furthermore, our approach efficiently exploits parallelization between multiple CPUs or GPUs.

Our simulation strategy offers a much-needed design and characterization tool for the growing community interested in applying DNA nanotechnology to engineer dynamic, reconfigurable devices[27,75] and nanorobots,[13] both of which would benefit from a rapid in silico prediction of free energy landscapes. This is especially the case for large origami structures, composed of multiple DNA scaffolds,[13] where conventional MD simulation is even more costly. Our technique would also be particularly suited for the better and faster calibration of nanoscopic mechanical probes,[24,51,76] especially in cases where simple analytical models may yield inaccurate results.[77] Finally, MetaD is not only relevant when exploring deformation in fully hydrogen-bonded motifs but could be also applied to free energy landscapes associated with hybridization/dehybridization by defining suitable collective variables, for example, in terms of the number of hydrogen-bonded nucleotides in the system.[78,79] In general, our approach will enable a faster and more detailed acquisition of information related to the mechanical behavior of nanostructures, improving the feasibility of simulation-informed design and facilitating a direct comparison of molecular modeling to experimental measurements.

Methods

oxDNA Implementation

The oxDNA stand-alone executable was extended to enable support for tabulated potentials and corresponding forces between the centers of mass of collections of particles on a one- or two-dimensional grid (CPU implementation) or a one-dimensional grid (CUDA implementation).[80] The source code was otherwise unchanged.

A Python interface was then used to launch multiple MD oxDNA simulations (replicas) in parallel, analyze distributions of collective variables, and update the bias. Each of the N replicas was initialized from a different location in collective variable space and simulated under the effect of the time-evolving bias Bt, shared between all replicas. After each MetaD cycle, corresponding to τ MD time steps, the bias was updated with N Gaussians placed at the instantaneous locations of each of the replicas in configuration space, as discussed above. The parallel replicas therefore share their history to construct an optimal potential, which leads to a more efficient exploration of the configurational space and to an N-fold speed up in bias convergence.[52] We also note that, at early times, the replicas repel each other, encouraging them to search different regions. However, this effect diminishes at later times. Parallel simulations were run on CPUs for dsDNA buckling, Holliday junction isomerization, and bistable unit isomerization, while GPUs were used for the origami-compliant joint. The number of replicas used in each case study is reported in Table S2.

Following convention,[48] the bias was defined on a grid, necessary to avoid a slowdown as the number of forces involved increases. The grid spacing δx has a value chosen to be at maximum one-fifth of the MetaD σ—values are given in Table S1. Our implementation is compatible both with Monte Carlo, where the potential felt by the particle is calculated from a bilinear interpolation of the gridded bias, and MD, where the force is calculated from the numerical derivative.

Simulations for the compliant origami joint required 36 h over four GPUs (total 144 GPU-hours, Nvidia P100 GPU 16GiB), while the other three case studies required between 36 and 60 h over 32 CPUs (total 1152–1920 CUP-hours, 2× Intel Xeon Skylake 6142 processors, 2.6 GHz 16-core). These time scales represent substantial improvements from unbiased simulations, which may require tens of thousands of GPU hours for characterizing the mechanical behavior of origami nanomachines.[38]

Choice of MetaD Parameters

The analysis by Laio et al.[59] and Bussi et al.[60] highlighted the influence of MetaD free parameters σ, A, and τ on the errors of inferred free energies and convergence time scales. These studies recommend optimal choices for the Gaussian width σ at a fraction of the system’s size in the collective variable space.[59] Because in this work we extract free energy surfaces from configuration sampling of simulations biased with the asymptotic Bt, rather than directly from the bias, we only followed the heuristic consideration that σ should be smaller than the length scale of the free energy features one wished to map, to prevent overbiasing. The ratio A/τ determines the (initial) rate of growth of the bias and, therefore, the convergence time, with larger A/τ implying faster convergence.[59] For well-tempered MetaD, A/τ is not critically important, as the amplitude of corrections decays exponentially. Using small values of τ (while appropriately rescaling A) reduces “discreteness” in potential deposition and errors in free energy estimates.[59] Because in our implementation the MetaD bias is computed and updated by a Python script, which then relaunches the stand-alone oxDNA executable after each cycle, computational inefficiencies emerge when reducing τ. These were considered in our choices of τ. The MetaD parameters for each of the systems simulated are summarized in Table S1.

Choice of Collective Variables

The choice of collective variables in MetaD should follow key criteria, detailed in Bussi et al.[48] First, the collective variables should be designed to force the system to explore the high free energy transition states one wishes to sample, which is done by ensuring that these states correspond to unique values of the collective variables, which are not accessible when the system occupies low free energy configurations. The application of this criterion is well-exemplified by the definitions of the two-dimensional collective variables for our Holliday junction and bistable unit case studies, where the two isomers are clearly separated from the intermediate transition states in the (x1,x2) planes. Second, and critical when mapping deformation free energy of large DNA nanostructures, one must ensure that the collective variables are properly coupled to the deformation modes one wishes to characterize. For example, if one would like to study bending of helices or bundles, the collective variables should be defined based on the coordinates of multiple nucleotides on different strands, to avoid that bias buildup leads to rupture of hydrogen bonds and nanostructure disassembly rather than bending.

Case Studies

Unless otherwise stated, simulations used the oxDNA2 force field with 0.5 M ionic strength. The sequence-averaged version of the force field was used in all cases except when mapping the free energy landscape of Holliday junction isomerization, where the sequence-dependent force field was instead adopted.[68] Molecular dynamics was used to sample configurations. A time step of 0.004 simulation units was used for the bucking dsDNA and Holliday junction case studies, while a time step of 0.005 simulation units was used for the compliant origami joint and bistable tile. To maintain a temperature of T = 300 K, an Andersen-like thermostat was used—time evolution is Newtonian, but every 103 time steps, a fraction of particles has their velocities drawn from a Maxwell–Boltzmann distribution. The fraction of particles subject to velocity changes correspond to a diffusion coefficient of 2.5 oxDNA units for all examples except the DNA origami joint, where the diffusion constant was set to 5 oxDNA units. Configurations were saved every 1 × 105 time steps for all systems except the origami, where they were saved every 1 × 103 time steps. Well-tempered metadynamics simulations were run with multiple walkers with parameters listed in Tables S1 and S2. Subsequently, the converged bias from those simulations was used in MD to verify correct convergence. Details concerning replicas and time scales are given in Table S2.

Kink-Induced Buckling in dsDNA

A DNA duplex of length 30, with sequence 5′-ATG CAC AGA TTA GGA CCA ACC AGG ATA GTA-3′, was initialized using the generate-sa.py script in the oxDNA software package. MetaD was run with a bias on the collective variable x, the distance between virtual particles at the centers of mass of the six nucleotides on one end of the duplex and the corresponding six at the other end. This choice was made to guarantee that the applied bias induces duplex bending, rather than dehybridization. Details of simulations are given in Tables S1 and S2.

To evaluate a reference free energy—ΔG(x)—the bias from the ΔT = 16T system was used in biased MD to acquire a large number of states (Table S2). Convergence of the free energy implied by the bias to the reference free energy is demonstrated in Figure S2. Here, an equilibration period of 1 × 108 time steps was used to decorrelate initial states. To demonstrate convergence, the ΔG(x) values were constructed from either the first half or the second half of the simulation, see Figure S2a. Differences between the two are substantially smaller that the width of lines used to make the plot.

Figure f features a two-dimensional free energy landscape. The quantity on the y-axis, Ufourth least stacked, was chosen to distinguish the buckled from the unbuckled state. This energy is defined by first acquiring the 5′ and 3′ stacking energies associated with each nonterminal nucleotide. Subsequently, the lesser of these two values was stored for each nucleotide. The fourth greatest (i.e., least negative) value in the list then defined Ufourth least stacked. Since the buckled state breaks two internal base-pair stacking interactions (where each is between a pair of nucleotides), this value rises to 0 if the duplex is buckled.

To illustrate the two distinct buckled and unbuckled states we plotted a kernel density estimator (KDE) with bandwidth 0.05 units—either nanometers or kBT in Figure f. This should not be overinterpreted other than to imply bistability when x is constrained to a value below ∼6 nm. For example, the buckled state has Ufourth least stacked exactly zero, so the density here is very high, and the exact free energy values will depend strongly on the KDE bandwidth.

Holliday Junction Isomerization

Holliday junctions were based on the J3 junction, as previously studied using single-molecule FRET measurements,[69] and an alternative coarse-grained force field.[71] Here, the junction is truncated so that arms are each 12 bps or ∼4 nm long, slightly over four D lengths for 100 mM ionic strength; sequences are given in Table S3. Truncation was necessary for faster simulation, and it is unlikely that nucleotides so far from the junction contribute to configuration probabilities. Structures were initialized using the MrDNA[31] software, then subsequently refined in the oxDNA-viewer software.[81] Four sets of simulations were run: either with a sequence-averaged force field at 500 mM ionic strength or with a sequence-specific force[68] field at either 500, 200, or 100 mM ionic strength. In each case, the counterion is modeled implicitly through control of the Debye length over which electrostatic screening operates. As simulations at reduced electrostatic screening are slower, runs at 100 mM ionic strength necessarily have fewer steps.

A two-dimensional collective variable was used in MetaD simulation, (x1, x2), as illustrated in Figure a. These variables were designed to clearly distinguish the two stacked isomers, where x1/2 takes small values and x2/1 takes large values, from the unstacked transition state where x1 and x2 both have high values, incompatible with the stacked isomers. Simulations were performed with the parameters from Table S1. After MetaD runs, a biased MD simulation was initialized from the terminal states of each of the six metadynamics walkers, each in eight replicas. The first 1 × 107 steps were discarded to allow for decorrelation. Contour plots in Figure c are acquired from histograms of biased MD. Convergence of the free energy implied by the MetaD bias and comparison to that acquired by histograms of biased MD is illustrated in Figure S3.

For identification of the states isoI, isoII, and the intermediate, the following criteria were used. For both the stacked-X conformers and the intermediate we required that all hydrogen bonds in the eight nucleotides adjacent to the junction were formed (internal energy < −1 kBT). For the stacked-X conformers, we further required that stacks were formed between pairs of neighboring arms at the junction, with a stack being said to occur if its internal energy is less than −5 kBT. We thus identified the isoI and isoII conformers based on which stacks were formed. As expected, stacked-X states, with two formed and two unformed stacks as illustrated in Figure a, are dominant in all explored conditions. The intermediate was defined as the state with no stacks formed but all hydrogen bonds present. The aforementioned state definitions were used to acquire the probabilities and free energies in Figure f. The isoI and isoII states for the sequence-averaged force field should be equal by symmetry, so the ∼3% difference in state probability is a reasonable indication of simulation convergence. Errors of estimates are given as one standard error, using six repeats initialized from different positions in the collective variable space.

Bistable Unit Isomerization

The bistable unit studied in the Supporting Information, Discussion 2, was designed in caDNAno.[11] The strand routing is given in Figure S4. A two-dimensional order parameter was used to bias MetaD, based on distances x1 and x2 as illustrated in Figure S5a. These distances were defined between groups of four nucleotides adjacent to each of the four nicks in the structure, making sure that the two isomers are clearly separated form the intermediate transition state on the (x1,x2) plane, as was done for the case of the Holliday junction.

Parameters for simulations are listed in Table S1 with durations listed in Table S2. Convergence is illustrated by the plots in Figure S6. The terminal states of six walkers were then used as initial states in MD simulations, biased with the converged Bt from MetaD. The biased MD simulations were used to construct the free energy distribution in Figure S5b (left) (3 × 107 steps discarded prior to collection for decorrelation).

Additionally, MetaD was run with a 1D collective variable, . To use this order parameter, analytical derivatives with respect to position were calculated. After convergence of the 1D bias, multiple MD simulations were then run with said bias. These biased MD runs were used to reconstruct the 2D free energy landscape given in Figure S5b (right).

DNA Origami Compliant Joint Bending

The DNA origami compliant joint studied here was based on an experimentally realized structure.[49] However, for reasons of speed, it was truncated, reducing the length of the helix bundles on the two sides of the joint. The experimental structure had six ssDNA scaffold sections routed across the compliant joint to apply a bending moment, whose magnitude could be controlled by the ssDNA length. Here we have removed these sections, relying instead on the MetaD bias to bend the joint. The caDNAno routing is given in Figure S7.

A collective variable was defined as detailed in Figure a and discussed in the main text. After they were designed in caDNAno, the structures were relaxed via a gradient descent to prevent large forces. Simulations were then run with the GPU-accelerated version of oxDNA.[80] For metadynamics, four walkers were run in parallel on separate GPUs associated with the same compute node. MetaD parameters were used as detailed in Table S1.

To establish a reference free energy (ΔG) to validate the convergence of MetaD predictions and later evaluate the ϕ distribution, the MetaD bias was frozen, and biased MD simulations were run. As detailed in Table S2, four different initial x configurations were used to generate samples, with six replicas run from each of those four initial configurations. These were run for 4× 106 steps to decorrelate replicas, followed by a production run (see Table S2 for details). To evaluate uncertainties for all estimates, the standard error from simulation runs initialized from different initial configurations was used.

To evaluate the distribution of ϕ, the reference beads illustrated in Figure a (yellow beads) were used. The top and bottom sections of the bundle each have 12 reference nucleotides selected. These are composed of two groups of six, one further and one nearer to the joint. Each of those groups of six corresponds to three base pairs, chosen to be adjacent to crossovers, guaranteeing that the applied bias does not induce unwanted structure disassembly. The distance along the bundle between the near six and far six was chosen to be 21 nucleotides (two helical turns), so that base pairs used as references have the same orientation. The center of mass of each of the four groups was acquired. For convenience, we use the notation x⃗neartop, x⃗fartop, x⃗nearbottom, and x⃗farbottom to denote these centers of mass. Every 20 000 steps, the locations of the centers of mass were saved. Subsequently, two vectors were defined.

The angle between these two vectors was used to define ϕ(0,180). This angle is not the ϕ that is then used in free energy calculations. It is important to then consider a definition of ϕ on (0°,360°), rather than (0°,180°) (so that, in Figure c, location E corresponds to a ϕ > 180°, while Figure , location F corresponds to a ϕ < 180°). Therefore, an additional vector was defined, corresponding to the direction into the page in Figure a (far left). We label this and define it using the two sets of beads furthest from the location of the junction. Looking at Figure a (far left), v⃗ortho corresponds to the average vector from the center of mass of the yellow beads nearest the reader to those into the page. Subsequently, the value ofwas acquired. This takes a value that is negative if ϕ is less than 180° and positive otherwise. Hence ϕ was acquired as

To evaluate the effect of ssDNA springs on the bending angle, the following approximations were used. There are two sets of three ssDNA that bridge the compliant DNA origami joint gap in the experimental system. These correspond to one set of three ssDNA segments that bridge the short gap, and one set that bridges the long gap (where the short and long gaps are illustrated in Figure c). The set of three ssDNA sections that bridge the short gap each have Nshort ssDNA nucleotides; the others have Nlong ssDNA nucleotides.

We then evaluated the free energy contribution from each of the ssDNA springs, ΔGWLC(L), using an analytical approximation for the free energy of a WLC.[82]Here Lp is the persistence length of ssDNA, which we took as 2 nm,[83] while L0 is the contour length. This was acquired from assuming that the contour length of ssDNA was 0.676 nm/nt.[84] One subtlety is that the number of nucleotides Nnts refers to is 1 greater than the number in the actual chain. This may seem surprising, but consider the case where there are 0 nucleotides in the ssDNA spring; there would still be one nucleotide of separation between the two sides of the joint. To evaluate the total free energy, we summed the contributions from the six chains, three of which contain Nshort nucleotides and three of which contain Nlong nucleotides.

There are six springs in total, so the total statistical weight used to compute averages isHere Bt(x) is the MetaD bias; i ∈ {0,1,2} indexes short springs, while i ∈ {3,4,5} indexes long springs. For each, L is the separation distance measured in simulation between attachment points. This weighted distribution was used to acquire both ⟨ϕ⟩ and σ(ϕ), as plotted in Figure e,f. Uncertainties here were acquired from the standard error over repeating this procedure for simulations run with four different initial conditions uniformly spaced in the range of x studied.