End-to-end attraction of duplex DNA.

Recent experiments [Nakata, M. et al., End-to-end stacking and liquid crystal condensation of 6 to 20 basepair DNA duplexes. Science 2007; 318:1276-1279] have demonstrated spontaneous end-to-end association of short duplex DNA fragments into long rod-like structures. By means of extensive all-atom molecular dynamic simulations, we characterized end-to-end interactions of duplex DNA, quantitatively describing the forces, free energy and kinetics of the end-to-end association process. We found short DNA duplexes to spontaneously aggregate end-to-end when axially aligned in a small volume of monovalent electrolyte. It was observed that electrostatic repulsion of 5'-phosphoryl groups promoted the formation of aggregates in a conformation similar to the B-form DNA double helix. Application of an external force revealed that rupture of the end-to-end assembly occurs by the shearing of the terminal base pairs. The standard binding free energy and the kinetic rates of end-to-end association and dissociation processes were estimated using two complementary methods: umbrella sampling simulations of two DNA fragments and direct observation of the aggregation process in a system containing 458 DNA fragments. We found the end-to-end force to be short range, attractive, hydrophobic and only weakly dependent on the ion concentration. The relation between the stacking free energy and end-to-end attraction is discussed as well as possible roles of the end-to-end interaction in biological and nanotechnological systems.

INTRODUCTION

Self-assembly properties of nucleic acids are vital to the basic functions of a biological cell and have been extensively exploited in biotechnology. DNA hybridization—self-assembly of complementary sequence single-stranded DNA (ssDNA) into a double helix—is a central biotechnological process (1), used, among others, in platforms for DNA detection (2), programmable assembly of DNA nanostructures (3,4), directional transport of cargo (5), molecular computing (6) and nanofabrication (7,8). Another process of outstanding importance is DNA condensation, where counterions transform electrostatic repulsion between naked DNA molecules into attraction, facilitating packaging of double-stranded DNA (dsDNA) in cell nuclei and viral capsids (9,10).

Recently, an entirely different type of DNA self-assembly was discovered: spontaneous end-to-end aggregation of short duplex DNA fragments into rod-like structures (11). When water was evaporated from solution containing a high concentration of short (6–20 bp) DNA fragments, liquid crystal phases were observed. Since the DNA fragments were nearly as wide as they were long, the observation of axial ordering could only be explained if the fragments formed rod-like supramolecules, suggesting end-to-end aggregation. Further experimental evidence of end-to-end association was obtained from the analysis of small angle X-ray scattering data from a system containing short DNA fragments and a divalent electrolyte (12,13). The second virial coefficient extracted from these data was shown to be positive for DNA fragments capped with a short hairpin (indicating overall repulsion) and negative for DNA fragments without such caps (indicating overall attraction). It was concluded that end-to-end attraction was large enough to overcome electrostatic repulsion in a divalent electrolyte. While side-by-side forces between long DNA molecules has been the subject of many experimental (14–16) and theoretical (17,18) studies, little is known about the conditions and microscopic mechanism of DNA association end-to-end. Furthermore, the effects of end-to-end attraction of duplex DNA in biological and technological processes are entirely unexplored.

Whereas traditional single molecule experiments have provided extensive information about DNA hybridization and side-by-side interactions (14,19), applying these tools to study end-to-end assembly is extremely difficult, as a DNA duplex cross section is just 5 nm2 and the effective concentration of DNA ends in a solution amenable to single-molecule studies is typically small. All-atom molecular dynamics (MD) method is well suited for the study of systems that share the length scale of short dsDNA molecules and can be used to probe the atomic origin of intermolecular forces (20). Here, we use the MD method to characterize end-to-end association of duplex DNA in unprecedented detail, elucidating the microscopic mechanism of spontaneous association, its free energy costs and the kinetic rates. At the end of the article, we discuss the relationship between our findings and the pertinent experimental observations.

MATERIALS AND METHODS

General simulation methods

All MD simulations were performed using the program NAMD (21), the parmbsc0 refinement of the AMBER-parm99 force field (22,23), the TIP3P model of water (24), standard parameters for ions (25), periodic boundary conditions, and particle-mesh Ewald (PME) full electrostatics with a PME grid density of ∼1 Å per grid point. Except where specified, van der Waals and short-range electrostatic energies were calculated using a smooth (10–12 Å) cutoff, and integration was performed using 1–2–4 fs multiple timestepping (21). The temperature was held constant using a Langevin thermostat (21) applied to all non-hydrogen atoms; the Langevin damping constant was set to 0.1 /ps. For simulations in the NPT ensemble, constant pressure was maintained at 1 bar using the Nosé-Hoover Langevin piston pressure control (26).

Each simulation reported in this study used one of the following three system types: elongated along the z-axis to minimize the amount of solvent around two DNA fragments (∼24 000 atoms, Figure 1a); isotropic to allow two DNA fragments to tumble freely (∼56 000 atoms, Figure 2a); and large and isotropic to allow unbiased interaction between 458 DNA fragments (∼1.4 M atoms, Figure 4a). The DNA sequence was poly(dA·dT) in all systems. Counterions were added to each system to neutralize the DNA charge prior to the addition of a number of ions corresponding to the reported molarity (100 mM, except where specified) of NaCl electrolyte. Steric clashes that were introduced during the assembly of each system were removed from each system through minimization using a conjugate gradient method (27). Equilibration was performed in the NPT ensemble, and subsequent production simulations were performed in the NVT ensemble, except where specified.

Figure 1.

Collapse of aligned dsDNA. (a) Simulation system containing axially aligned duplex DNA. Each DNA fragment is free to move along and rotate about the common axis. The DNA duplexes (blue and green) are shown in van der Waals representation; sodium and chloride ions are shown as yellow and cyan spheres; water is not depicted. An animation illustrating spontaneous end-to-end collapse of duplex DNA is available in Supplementary Data. (b and c) The end-to-end distance (b) and the relative azimuthal angle ϕ (c) of two duplex DNA in representative simulations of the end-to-end collapse. Data from the same pair of simulations are plotted in (b) and (c). (d and e) Scatter plot showing the relative azimuthal angle ϕ at the time of collapse (d), and at the end of simulation (e). One data point is shown for each of 36 simulations of blunt-ended (black circles) and 5′-phosphorylated (red squares) dsDNA fragments in 100 mM NaCl electrolyte.

Figure 2.

Stability of the end-to-end DNA assembly. (a) Simulation system containing two spatially unrestrained dsDNA fragments in 100 mM NaCl electrolyte. Both fragments are free to rotate and move about the simulation box. The system is drawn as in Figure 1a. An animation illustrating a typical MD trajectory is available in Supplementary Data. (b and c) End-to-end junction of 5′-phosphorylated dsDNA with ϕ=−20° (b) and blunt-ended dsDNA with ϕ=180° (c). (d and e) End-to-end distance (d) and relative azimuthal angle, ϕ, (e) of the DNA fragments in three unrestrained MD trajectories.

Figure 4.

Spontaneous aggregation of duplex DNA into long rod-like structures. (a) A system containing 458 duplex DNA fragments placed at random. NaCl solution is shown as a semi-transparent molecular surface. (b and c) Initial (b) and final (c) conformations of the DNA fragments that composed the longest 10 aggregates at the end of a 260-ns MD simulation. DNA fragments forming each aggregate are shown in a different color. (d) The instantaneous number of aggregates, N(s), formed by s DNA fragments in the MD simulation. The lines show the equilibrium distribution of the aggregate length according to the reverse step-growth polymerization model (37) for specified values of Gbind. A movie of the simulation trajectory is provided in Supplementary Data.

Collapse of aligned dsDNA

Thirty-six systems were built using the anisotropic unit cell. The dimensions of each system were chosen to provide a minimum distance of 2 nm between the surfaces of the DNA fragments across the periodic boundary, which should accommodate a majority of screening counterions for NaCl solutions of 100 mM or greater concentration (Debye length ≤1 nm). Steric clashes were removed through 3000 minimization steps. Each system was subsequently equilibrated for 65 ps with the DNA backbone atoms harmonically restrained to their initial positions. Axial alignment of the DNA fragments was enforced by harmonically restraining each phosphorous atom of the DNA to the surface of an 11-Å radius cylinder (with spring constants of 139 pN/nm per atom). The DNA fragments could translate along the axis of the cylinder and rotate azimuthally. The starting conformation was characterized by a 20.5 Å end-to-end separation, which we define as the distance between the centers of mass of the nearest terminal base pairs, taking the periodic boundary condition into account. The relative azimuthal angle ϕ of the terminal base pairs was defined as the angle between the projections of the vectors connecting the O5′ and O3′ atoms of the terminal base pairs into the plane normal to the common DNA axis. For two consecutive base pairs in a B-DNA helix, ϕ ≈ 36°, depending on the sequence.

Stability of the end-to-end complex

Systems were built by placing collapsed end-to-end DNA assemblies in an isotropic volume of 100 mM NaCl electrolyte (Figure 2a). Steric clashes with the solvent were removed through minimization having the DNA backbone restrained. Subsequent simulation was performed in the NPT ensemble.

Mechanics of end-to-end dissociation

The coordinates of the system containing 5′-phosphorylated DNA fragments in a 100 mM electrolyte were taken after 140 ns of simulation described in the section ‘Stability of the end-to-end complex’ to provide the initial conformation for simulations of rupture of the DNA fragments. The NVT ensemble was used during these simulations; Langevin thermostat was applied to water oxygens and ions. A harmonic spring of ks=4000 pN/nm was used to produce the rupture, by increasing its rest length at a rate of 0.4 or 0.2 Å/ns. The work performed at both rates were in good agreement. Simulations of the end-to-end fragments having complementary single-stranded overhangs are described in Section 1.4 of Supplementary Data.

Potential of mean force of axially aligned DNA duplexes

Umbrella sampling simulations were performed using the anisotropic systems (Figure 1a) and two simulation protocols different by the method used to set up initial systems and the alignment restraints. Both protocols enforced the end-to-end distance r using a harmonic spring of ks = 4000 pN/nm for 3.5 < r < 12 Å in 0.5-Å intervals and ks = 1000 pN/nm for 13 < r < 19 Å in 1.0-Å intervals. The first protocol was used to provide the estimate of the potential of mean force (PMF) (Figure 3). The initial conformation for each simulation was obtained by placing the DNA fragments a specified distance r apart at one of the four ϕ = 0, 90, 180 and 270° (four simulations for each r). The systems were equilibrated for 2 ns before data accumulation during production simulations lasting ∼16 ns. In the second protocol, which we used to compute the relative binding free energies (Table 1), the initial conformations were generated iteratively by shifting the minimum of the restraining potential in steps and followed by 0.5-ns equilibration, starting from the final frames obtained in the simulations described in the section ‘Collapse of aligned dsDNA’. Subsequently, each system was equilibrated for at least 2.5 ns before accumulation of data during production simulations lasting 7.5–15 ns. Axial alignment was maintained as described in the section “Collapse of aligned dsDNA”, using ks = 13.9 and 139 pN/nm for the first and second protocols, respectively. In the second protocol, a torque pointing along the common DNA axis was distributed among the phosphorous atoms of each DNA molecule to restrain ϕ about −20°, 36° or 180°, with a spring constant of 219.4 pN nm/rad2, which roughly corresponds to an 8° root mean squared fluctuation.

Figure 3.

Representative dependence of the effective force (red) and the free energy (black) of two axially aligned DNA fragments on the end-to-end distance. The data result from 100 independent simulations of two 5′-phosphorylated fragments, the end-to-end distance of which was maintained at a specified value by a harmonic spring. Additional restraints were applied to maintain the axial alignment. The strength of the restraints was found to affect the values but not the general shape of both curves (see text). The image in the background illustrates the simulation method. The DNA fragments were immersed in 100 mM NaCl electrolyte, and were free to rotate about their helical axes.

Table 1.

Relative free energy change, ΔΔG, upon formation of the end-to-end complex

5′ chemistry	ion concentration (M)	ϕ	ΔΔG (kcal/mol)
Phosphoryl	0.1	36°	0
Phosphoryl	1.0	36°	0.4
Phosphoryl	1.0	−20°	−1.8
Phosphoryl	1.0	180°	0.1
Hydroxyl	1.0	36°	2.3
Hydroxyl	1.0	180°	1.2

Here, the free energy change ΔG is approximated by the minimum of the end-to-end PMF obtained from umbrella sampling simulations (Supplementary Figure S4). The values of ΔG are given relative to the value measured for the system containing 5′-phosphorylated fragments end-joined in the 36° orientation in 100 mM NaCl. In each simulation, the relative azimuthal orientation of the DNA fragments was enforced using harmonic restraints (see ‘Material and Methods’ section). The application of such restraints introduced a bias to the estimates of ΔG. As all simulations employed the same restraints, this bias was assumed to cancel out in the calculation of ΔΔG.

Spontaneous assembly of long end-to-end aggregates

The system (depicted in Figure 4a) was assembled through the sequential placement of 458 DNA fragments in a cubic volume (250 Å on each side) of 100 mM NaCl electrolyte. To place a DNA fragment, trial positions and orientations were randomly selected until the DNA coordinates did not clash with any previously placed fragments. During the first 50 ns of equilibration, the system shrank to its equilibrium size of 238 Å on each side. To improve computational efficiency, a 7–8 Å cutoff was used along with 2–2–6 fs time stepping scheme.

RESULTS

Spontaneous end-to-end association of duplex DNA was observed in the simulations of two (dA·dT)10 fragments constrained to diffuse along a common axis in a volume of 100 mM NaCl electrolyte. Figure 1a illustrates the initial state of a typical simulation system. Figure 1b plots the distance between the DNA fragments versus time for two simulation systems differed by the termination of the DNA's 5′-ends. The DNA fragments were observed to freely diffuse along the common axis until the end-to-end distance fell below ≈8 Å, whereupon the fragments rapidly collapsed into an end-to-end bound complex. The relative azimuthal orientation of the DNA fragments ϕ continued to change after the collapse, switching between several metastable orientations (Figure 1c). We define the relative azimuthal angle ϕ as the angle between the projections of the vectors connecting the O5′ and O3′ atoms of the terminal base pairs into the plane normal to the common DNA axis (see ‘Materials and Methods’ section).

In the final conformation adopted by the blunt-ended fragments, the 5′ to 3′ direction of the backbone was discontinuous at the end-to-end junction. In the case of the 5′-phosphorylated fragments the 5′ to 3′ direction of the backbone was continuous at the end-to-end junction as though the backbone of a continuous 20 bp B-DNA had been cut.

To determine the statistical significance of the above observations, we performed 17 additional simulations for each system type, different by the relative azimuthal orientation of the DNA fragments at the onset of the simulation: ϕ = i × 20°, where i = 1, … , 17. In all simulations, we observed collapse of the DNA fragments into an end-to-end bound complex. Figure 1d plots the relative azimuthal orientation at the time of collapse against the time to collapse. The collapse of blunt-ended fragments occurs irrespective of their azimuthal orientation, whereas formation of a 5′-phosphorylated end-to-end assembly did not occur with ϕ between 90 and 230°.

After collapse, ϕ continued to evolve, reaching the states characterized in Figure 1e. The clustering of ϕ values around −20, 36 and 180° indicates the three preferred binding states. At ϕ ≈ 36°, the conformation of the bound complex is similar to that of a continuous B-form DNA. At ϕ ≈ −20°, the backbones of the terminal base pairs overlap slightly (Figure 2b). At ϕ ≈ 180°, the 5′-ends of the fragments are in direct contact (Figure 2c). The conformations of the end-to-end assemblies with ϕ = −20 and 36° are similar, differing primarily in the order in which the 5′- and 3′-termini are lapped and in their base stacking geometry. Thus, relative to the coordinates of two consecutive base pairs in a canonical B-DNA helix, the terminal base pairs forming an end-to-end junction have a time-averaged root mean squared deviation (RMSD) of 2.1 and 0.9 Å for the ϕ = −20 and 36° conformations, respectively. For reference, base pairs in the middle of one of the DNA fragments had an RMSD of 0.7 Å.

The preference for these three orientations suggests a hydrophobic origin of the attractive force, as such conformations reduce exposure of the hydrophobic bases to water. The relative orientations of the blunt-ended DNA were nearly equally split between the −20 and 180° states. More than 50% of the 5′-phosphorylated fragments formed the −20° state, 35% formed either the 36° or 180° state (three systems each) and 12% formed the state with ϕ ≈ 100°. We attribute such preferential alignment of the 5′-phosphorylated fragments to the electrostatic repulsion between the terminal phosphate groups, which apparently renders the 180° orientation energetically less favorable than the −20° one. The free energy difference between these bound states is discussed below.

DNA fragments initially forming a bound state were simulated in the absence of any restraints using the isotropic system shown in Figure 2a. Three systems were constructed: one containing 5′-phosphorylated DNA bound with ϕ = − 20° (Figure 2b) and two containing blunt-ended DNA fragments bound at ϕ = 180° (Figure 2c) and ϕ = −20°. All three systems contained 100 mM NaCl electrolyte.

The plot of the end-to-end distance, Figure 2d, indicates that all three assemblies remained bound for the entire duration of the simulations (>200 ns). The standard deviation of the end-to-end distance in the 5′-phosphorylated system was 0.68 Å, twice less than that of the blunt-ended systems. The greater stability of the 5′-phosphorylated complex may be due to hydrogen bonds that were observed between the 5′-terminal phosphate of one DNA fragment and the 3′-terminal hydroxyl of the other fragment. The plot of the relative azimuthal angle reveals that the 5′-phosphorylated system maintained the same stable conformation at ϕ ≈ −20°, depicted in Figure 2b. Starting from a similar conformation, the blunt-ended complex underwent two sudden rotations at 70 and 140 ns that brought ϕ from −20° to −110° and to 180°; the relative orientation continued to evolve after that.

In the simulation of the 180° blunt-ended complex, the terminal base pair of one of the fragments ruptured after 33 ns. During the next few nanoseconds, Watson–Crick pairs within that fragment stochastically broke and re-annealed, propagating the unpaired base toward the opposite end of the DNA fragment, and slipping the entire DNA strand with respect to the other by 1 bp. Despite the slippage, the DNA fragments remained stably bound.

The lifetime of a bound complex sharply decreases when an external force is applied to disrupt it (28). Thus, under a constant force of 150 pN directed along the common axis of the DNA fragments, the assembly remained intact during a 50-ns simulation but dissociated within 5 ns under a 200 pN force (see Section 1.1 in Supplementary Data for simulation details).

To determine the dissociation pathway, the rupture force and the mechanical work required to dissociate the assemblies, the DNA fragments were subjected to the force of a harmonic spring whose equilibrium-extension length was increased at a constant rate. Spring-driven rupture of this sort has been used extensively in the study of proteins (29,30). Rupture was induced by pulling apart the fragments either along (axial stretching) or perpendicular to (transverse shearing) their common symmetry axis, applying the force either to the nearest ends or to the centers of mass (CoM) of the fragments. At least four simulations were performed for each protocol (14 in total). Animations in Supplementary Data illustrate typical simulation trajectories. In all cases, rupture was observed to occur by sliding of one terminal base pair relative to the other and was preceded by stretching of the duplex in the case of the CoM axial pulling. Although the three rupture protocols yielded different dependencies of the force on the separation distance (Supplementary Figures S1 and S2), the average work performed was 9.4±1.5 kcal/mol, independent of the rupture protocol. The typical rupture forces varied between 100 and 200 pN and were considerably lower for CoM pulling. Inclusion of short overhangs of complimentary sequence ssDNA at the ends of the fragments increased the work required to rupture the end-to-end assemblies (for details, see Section 1.4 and Figure S3 in Supplementary Data).

PMF of axially aligned DNA duplexes

To improve our estimates of the force and free energy of the end-to-end interaction, 100 variants of the system containing two-DNA fragments in 100 mM NaCl electrolyte, Figure 1a, were simulated with a constant end-to-end distance enforced by a harmonic spring potential. Each simulation explored a unique combination of the end-to-end distance and the azimuthal angle and lasted 18 ns (aggregate simulation time was 1.86 μs). The DNA fragments were kept aligned by weak harmonic restraints (see ‘Materials and Methods’ section) that allowed the terminal base pairs to shear.

Figure 3 shows the dependence of the effective end-to-end force on the end-to-end distance and the PMF reconstructed from this set of simulations by the weighted histogram analysis method (31,32). The PMF can be thought of as the change of free energy along a chosen coordinate. The force sharply increases with the end-to-end distance between 3.5 Å—the distance between consecutive base pairs in a DNA helix—and 6.5 Å, the separation allowing water molecules to penetrate the volume between the ends of the fragments. The force rapidly decreases as the end-to-end distance exceeds 6.5 Å and becomes slightly repulsive after ∼13 Å (2–10 pN). Thus, the end-to-end force has a large but very short-range attractive component caused by the hydrophobic effect and a much smaller long-range repulsive component that originates from screened electrostatic interactions between the DNA fragments (19).

A variation of the above protocol (described in Section 2 of Supplementary Data) was used to calculate the PMF for DNA fragments different by their terminal chemistry and relative azimuthal orientation and for several concentrations of the surrounding electrolyte. Table 1 lists the change in the depth of the PMF minima (Δ(min[PMF] – PMF(∞)) relative to its value for the 5′-phosphorylated fragments with ϕ = 36° in 100 mM NaCl. These calculations (detailed in Supplementary Data) demonstrate that increasing the electrolyte concentration from 0.1 to 1 M has negligible effect on the PMF. Among the three most likely orientations that 5′-phosphorylated fragments form, the ϕ = −20° state has the lowest free energy, in good agreement with the ϕ = −22° angle frequently observed in the crystal structures of poly(dA·dT) oligonucleotides (33). The 5′-phosphorylated fragments exhibit deeper minima than the blunt-ended fragments. For the latter, the conformation of ϕ = 180° is preferred over ϕ = 36°. These variations in the depth of the PMF are consistent with the occupancy of bound states observed in our simulations of spontaneous collapse of aligned DNA fragments (Figure 1e).

Standard binding free energy of end-to-end association

We computed the standard free energy of binding Gbind using a variation of the method described previously (34). For a system of two DNA fragments, Gbind determines the fraction of time the fragments form a bound state, which can be, in principle, observed directly in an all-atom MD simulation. However, the dissociation rate of the end-to-end assembly was found to exceed 200 ns (Figure 2d) and therefore using such a brute force approach was not possible. Equivalently, Gbind can be obtained from the logarithm of the equilibrium binding constant, which is the ratio of the kinetic rates of end joining and rupture (kon and koff, respectively).

For our calculations of Gbind, we considered a process consisting of the following four steps (Supplementary Figure S5). First, we evaluated the free energy cost of enforcing axial alignment restraints on a pair of infinitely separated DNA fragments. Second, we computed the cost of bringing a pair of DNA fragments from infinity to the maximum-separation state considered in our simulations (CoM–CoM distance of 52 Å). Third, we determined the free energy cost of forming the end-to-end complex of two axially aligned DNA fragments. Finally, we evaluated the cost of releasing the axial alignment restraints from the end-to-end assembly. The sum of these terms yielded the free energy change upon formation of the end-to-end DNA complex Gbind = −6.3 ± 1 kcal/mol for a DNA concentration of 1 M in 120 mM NaCl. Since a pair of DNA ends has only one binding configuration, we can express the standard binding free energy in terms of DNA ends, so that  kcal/mol. A complete description of the methods used in this section is provided in Section 3.1 in Supplementary Data.

The above value for Gbind represents our best effort to quantify the strength of the end-to-end interaction. In general, Gbind may depend on the DNA sequence. Although the question of sequence dependence is enticing, we defer investigations of the sequence dependence of the end-to-end interaction to future studies.

The rate of dissociation of an end-to-end assembly, koff, can be estimated from the PMF and diffusion coefficient, D, by computing the mean first passage time under the assumption that reannealing does not occur (35,36) (see Section 3.2 in Supplementary Data for details). From the simulations performed in the section ‘Collapse of aligned dsDNA’, we estimate D ∼25 Å2/ns. Due to uncertainty in the location of the barrier peak, we obtain the range 170 000–860 000 ns, with a best estimate of 480 000 ns. The use of axial alignment restraints in our calculations of the PMF limits pathways ordinarily available to rupturing DNA, so these values likely represent an upper bound for .

Assuming that end-to-end binding is limited by translational diffusion, the upper bound estimate for kon is 4πDR0 ≈ 7 (ns M)−1, where R0 = 37 Å is the CoM distance between a pair of DNA fragments at which binding is expected to occur. The true value of kon should be smaller because the DNA fragments must be axially aligned for binding to occur. Furthermore, long-range electrostatic repulsion may reduce the value of kon. Our estimates of Gbind and koff suggest a range for kon of 0.03–0.16 ns−1 M.

Our results until this point described the interaction of two DNA fragments in isolation. To simulate multifragment aggregation, 458 DNA fragments, each 10 bp in length, were placed in a cube of 100 mM NaCl solution (23.8 nm on each side) to form the system shown in Figure 4a. During a 260-ns MD simulation, the DNA fragments diffused about their initial positions and interacted with their neighbors to form aggregates up to 11 DNA fragments (110 bp) in length. The DNA fragments that formed the longest 10 aggregates are shown at the beginning, Figure 4b, and the end, Figure 4c, of the simulation. The number of aggregates of a given length at three different instances of the MD trajectory is shown in Figure 4d. The plot reveals rapid growth of end-to-end aggregates and roughly exponential distribution of the lengths of the aggregates. The number of aggregates did not reach a steady state by the end of the simulation (Supplementary Figure S6).

We model the process of end-to-end aggregation using a simplified reversible step-growth polymerization model, which has an analytical time-dependent solution that relates the mean aggregation number, 〈N〉 to Gbind as a function of DNA concentration, c (37) (see Section 4 in Supplementary Data). The DNA concentration in our multifragment system was ∼56 mM. Using our prediction of Gbind = −6.3 kcal/mol, we obtain 〈N〉 = 39 for the mean aggregation number at equilibrium for this system.

The above simulation partially mimics the experimental assay of Nakata et al. (11), which involved very dense fluids of short DNA fragments. For 10 bp DNA fragments, Nakata et al. found a transition from an isotropic phase to a nematic liquid crystal phase at c = 875 mg/ml (∼150 mM). Through a combination of calibration experiments and theory, the authors estimated 〈N〉 = 9 at the isotropic–nematic phase transition. Under the model of reversible step growth polymerization, 〈N〉 = 9 at c = 150 mM corresponds to Gbind = −3.87 kcal/mol. In contrast, using our prediction of Gbind = −6.3 kcal/mol yields 〈N〉 = 64 at a DNA concentration of 150 mM.

In the reversible step growth polymerization model, the total number of associated molecules is described by a second-order kinetic equation. Since unbinding of DNA fragments is negligible in our system, kon can be extracted by fitting to the data, where N(t) is the total number of end-to-end bound DNA fragments at time t, N0 = N(0) and c0 is the initial DNA concentration. In our simulation of spontaneous aggregation, the association rate reduced from kon = 0.37 ns/M observed within the first 50 ns to 0.069 ns/M for the 60–250-ns interval (Supplementary Figure S6). The change in association rate occurred as longer end-to-end aggregates formed (>3 duplex fragments) and the rotational and translational diffusive motions of shorter aggregates slowed.

The aggregation simulation provides a lower bound estimate for the dissociation rate koff. During the course of the simulation, 307 DNA ends were bound for 217 ns on average and unambiguous unbinding was observed for only one end-to-end associated complex. A few events were observed in which partial unbinding occurred via rupture of the Watson–Crick base pair of terminal nucleotides, as well as one instance where a bound DNA fragment was transferred to an unbound fragment; we neglect these events for the subsequent analysis. A Poisson distribution yields the probability that exactly one end-to-end bound DNA pair would rupture during the simulation given a value of koff. For  <1900 ns, it is likely that more than one rupture would have occurred. For >59600 ns, it is likely that no rupture would have occurred. Thus, the probability that exactly 1 rupture occurred is >10% for –596 000 ns. The greatest probability of exactly 1 rupture occurring is 37% for  ns. The ranges estimated for kon and koff suggest −Gbind in between 4.4 and 7.6 kcal/mol, with a most probable range of 5.2–6.2 kcal/mol [kon = 0.069–0.37 ns/M and  ns]. Alternative statistical analysis of the results yields a similar range for Gbind and is described in Section 4 of Supplementary Data.

DISCUSSION AND CONCLUSIONS

Stability of a dsDNA molecule can be conveniently described as a sum of base stacking and base pairing energies contributed by individual nucleotides that constitute the molecule. It is tempting to conceptually equate the base stacking interactions within a continuous molecule with the base stacking interactions that drive the assembly of two disjoint DNA duplexes. Thus, the unified nearest neighbor parameters, which can predict the energy for DNA hybridization based on DNA melting data, suggest Gbind = −16.94 + 2 × 6.94 = −3.06 kcal/mol (38) for the association of two 10 bp poly(dA·dT) DNA fragments into a continuous 20-bp molecule. Such a simple calculation may be, however, flawed as additional conformational flexibility afforded by the lack of phosphodiester bonds at the end-to-end interface should allow the base stacking geometry to be optimized, magnifying the base stacking contribution to the free energy. Accordingly, the average interaction energy between adjacent base pairs with ϕ = −20° was measured to be 2 kcal/mol lower in bases forming the end-to-end junction than in bases in the middle of one of the DNA fragments (Supplementary Figure S7). This finding is in good agreement with a survey of crystallographic structures that found DNA fragments formed of AT base pairs to stack with ϕ = −22° (33). We note that the stacking geometry may depend on the sequence of the DNA termini.

The base stacking free energy has been experimentally quantified by introducing a dangling nucleotide or a nick (a cut in the backbone of one strand) to a DNA molecule and observing the change in melting temperature (39–41), and by observing the mobility of a nicked DNA molecule relative to intact DNA and DNA with a gap (42,43). These experiments provide estimates for the base stacking free energy between −0.65 and −2.0 kcal/mol for stacks formed by thymine and adenine. However, the extraction of the free energy values is indirect with these methods.

Here, we provide the first direct estimate of the standard binding free energy of end-to-end association of DNA fragments, obtaining Gbind = −6.3 ± 1 kcal/mol. The simulations reported here are similar to the liquid crystal condensation experiments of the Clark and Bellini groups (11), which estimated −3.8 kcal/mol for the end-to-end free energy, in reasonable agreement with our estimate. The small-angle X-ray scattering experiments of the Pollack group demonstrated that the end-to-end interaction dominates over electrostatic repulsion in a divalent electrolyte (12,13), which, unfortunately, is not sufficient to estimate the standard free energy of end-to-end binding. We note that the value of Gbind obtained in this study is larger than values reported in experiments. Having employed multiple methods, we are confident that the range of values obtained for Gbind accurately reflects the standard binding free energy within the limits of the molecular force field used in our study. Nevertheless, we cannot rule out the possibility that the present MD force field somewhat exaggerates the interactions driving end-to-end self-assembly of duplex DNA.

It is interesting to note that hydrophobic interactions between DNA bases and inorganic materials can be significantly stronger than the stacking energies observed in biochemical assays. AFM experiments indicate that ssDNA adheres to graphite with free energies of −4.9 and −6.8 kcal/mol per nucleotide for cytosine and guanine, respectively (44). For comparison, a computational investigation of hydrophobic interactions revealed a free energy of −55 kcal/mol for the adhesion of a pair of 11 × 12 Å2 graphene sheets, or about 10 times the energy per unit surface area (45). Interactions of similar strength were found to promote DNA–fullerene and ssDNA–carbon nanotube association (46,47).

Given the relatively large free energy of the end-to-end interaction, we pose the following question: why has end-to-end association only recently been observed? In most biological and nanotechnological settings, the concentration of DNA ends is too low for end-to-end association to be statistically significant, and the lifetime of end-to-end interactions (∼70 μs) is somewhat shorter than the temporal resolution of many experimental techniques. To illustrate this point, we plot in Figure 5 the fraction of bound DNA ends versus the concentration of free DNA ends in chemical equilibrium. The concentration of DNA ends is relatively high in DNA cyclization assays that measure the fraction of DNA molecules bent into a circle. The J-factor, which represents the concentration of one end in the proximity of the other, has a maximum value of ∼10−4 mM (48), which is two orders of magnitude too small to promote cyclization of a significant fraction of blunt-ended DNA. The introduction of sticky ends increases the interaction energy by a few kcal/mol so that the DNA cyclization can be observed.

Figure 5.

The effect of end-to-end attraction in different DNA systems. The binding free energy (black) and the fraction of bound DNA ends (red) are plotted against the reference concentration of DNA ends. Background images schematically illustrate four DNA systems in which the end-to-end attraction may or may not play a role. From top left to bottom right: the maximum local concentration of DNA ends [i.e. the J-factor (51)] is about two orders of magnitude too low to induce an observable fraction of blunt-ended DNA circles of any length (orange); translation and rotational confinement of DNA ends at dsDNA breakage [e.g. by the protein Ku (PDB:1JEY)] will promote binding of the DNA ends, which likely aids repair of DNA during non-homologous end joining (49) (purple); the structure factor obtained from small-angle X-ray scattering experiments of short DNA duplexes in a divalent electrolyte reveals end-to-end attraction (12) (blue); at very high DNA concentrations, long DNA aggregates form and align in liquid crystal phases (11) (green).

The large standard binding energy of end-to-end association implies that end-to-end interactions will be important in systems containing a high concentration of DNA ends. For instance, a remarkable new method of creating patterned, self-assembled structures out of DNA—termed DNA origami (3)—introduces many nicks along a path of DNA, which may enhance the stability of the resultant pattern. We speculate that broadened awareness of end-to-end association will influence the development of nanotechnologies where the end-to-end interaction can be used advantageously, such as with DNA origami, or can pose a limitation, such as for DNA microarrays.

In cells, double-stranded DNA breakage, which poses a mortal threat, results in nearby blunt or sticky DNA ends. During non-homologous end-joining—the dominant repair pathway for such breaks in multicellular eukaryotes (49,50)—the ruptured DNA ends are held together by proteins such as the Ku heterodimer or DNA-PK, or by nucleosome interactions until damaged DNA can be removed and ligation of the DNA backbone occurs (49). Since DNA attracts end-to-end, it is not necessary that this complex, whose microscopic structure is not yet known, hold the DNA ends in strict alignment. It is sufficient for the ends to be held proximally so that the effective concentration of DNA ends is large enough to promote end-to-end association (Figure 5, purple), whereupon ligation may occur. The free energy of dsDNA end-to-end association found in this study suggests that placing DNA ends in a sphere of 3-nm radius will produce an end-to-end associated state ∼95% of the time.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online: Supplementary Methods, Results and Figures 1–7, Supplementary Animations 1–6 and Supplementary References [52-61].

FUNDING

Center for the Physics of Living Cells through a grant from the National Science Foundation (PHY-0822613). The supercomputer time was provided at TeraGrid (MCA05S028). Funding for open access charge: NSF.

Conflict of interest statement. None declared.