De novo reconstruction of DNA origami structures through atomistic molecular dynamics simulation.

The DNA origami method has brought nanometer-precision fabrication to molecular biology labs, offering myriads of potential applications in the fields of synthetic biology, medicine, molecular computation, etc. Advancing the method further requires controlling self-assembly down to the atomic scale. Here we demonstrate a computational method that allows the equilibrium structure of a large, complex DNA origami object to be determined to atomic resolution. Through direct comparison with the results of cryo-electron microscopy, we demonstrate de novo reconstruction of a 4.7 megadalton pointer structure by means of fully atomistic molecular dynamics simulations. Furthermore, we show that elastic network-guided simulations performed without solvent can yield similar accuracy at a fraction of the computational cost, making this method an attractive approach for prototyping and validation of self-assembled DNA nanostructures.

INTRODUCTION

Self-assembly offers a route to arranging matter at the nanoscale with high accuracy and at low cost. One successful application of the self-assembly principle is DNA origami (1), a method of folding a long DNA strand into a pre-determined three-dimensional (3D) shape (2,3). Such self-assembled DNA nanostructures are already finding practical applications in single-molecule studies of enzyme activity (4,5), in vivo sensing of ion concentrations (6) biomimetic membrane channels (7,8) and in nanometer-precision placement of proteins, metal nanoparticles, and carbon nanotubes (9–11). Further development of the DNA origami method will likely lead to larger, more complex, and more useful self-assembled systems. For example, it is not hard to imagine a DNA-based precision drug delivery system (12–16), as there already exist DNA origami devices that can enter a living cell (17), detect lung-cancer-specific microRNAs (18) and perform logical operations to decide on the release of a payload (19,20).

Structural characterization of DNA origami objects is an essential part of their design and optimization process. The leading DNA origami design tool, caDNAno (21), represents DNA helices as perfectly rigid rods, Figure 1A. In solution, however, deviations from idealized design can be expected. Transmission electron microscopy (3,22), super-resolution optical imaging (23), fluorescence resonance energy transfer (9) and magnetic tweezers (24) have been applied to infer information about the in situ structure and dynamics of DNA origami objects. The only experimentally derived 3D structure of a DNA origami object was obtained using the cryo-electron microscopy (cryo-EM) method (25), revealing considerable deviation of the in situ structure from the idealized design, Figure 1B. A pseudo-atomic model of the pointer object, Figure 1C, was derived from the cryo-EM data by fitting an atomic-scale model into the electron density map (25).

Figure 1.

Structural models of DNA origami. (A) A three-dimensional (3D) model of the pointer object (25) built according to its idealized design. Each DNA helix is represented by a cylinder. (B) Cryo-electron microscopy reconstruction of the DNA pointer structure (25). The object's structure is characterized by a 3D electron density map visualized in the figure as a surface of constant electron density. (C) All-atom model of the DNA origami object. Atoms comprising the scaffold and staple strands are shown as blue and white van der Waals (vdW) spheres.

Recently, we demonstrated the utility of molecular dynamics (MD) simulations for providing realistic, atomic-scale models of simple DNA origami constructs (26,27). Here, we assess the accuracy of the MD method for de novo prediction of the atomic structure of large DNA origami objects through direct comparison with the results of cryo-EM reconstructions. We find the MD method capable of reproducing the structural features of the pointer object with accuracy approaching that of state-of-the-art cryo-EM microscopy. Furthermore, similarly accurate structures were obtained when solvent was treated implicitly.

MATERIALS AND METHODS

General simulation protocol

All MD simulations were performed using the program NAMD (28), the CHARMM36 (29–31) parameters for nucleic acids, the TIP3P model of water (32), a custom parameterization for ions (26,33), periodic boundary conditions, smooth truncation of Lennard-Jones and short-ranged Coulomb interactions at 10 Å with the switching function active above 8 Å, and, except where specified, particle-mesh Ewald (PME) for long-ranged electrostatics with a PME grid density of about 1 Å per grid point.

The temperature was held constant at 300 K by applying Langevin forces (34) to all non-hydrogen atoms; the Langevin damping constant was set to 1 ps−1, except where specified. All production simulations were performed in the NVT ensemble (constant number of atoms N, volume V, and temperature T).

All-atom explicit solvent simulations

To build an all-atom model of the pointer object, the caDNAno design of the object (provided by H. Dietz) was converted to an atomic model using our custom conversion script cadnano2pdb as previously described (26). Magnesium–hexahydrate complexes were randomly placed in proximity of DNA using an in-house visual molecular dynamics (VMD) script (33). Water was added using the solvate plugin of VMD (35). A pre-defined number of randomly chosen water molecules was replaced by sodium and chloride ions using the autoionize plugin of VMD. Each system was energy minimized for 2400 steps, and then equilibrated in the NVT ensemble having all non-hydrogen atoms of DNA harmonically restraining to their idealized structure coordinates, which allowed water and ions to attain equilibrium distribution. The spring constant of each restraint k = 1 kcal/mol Å2.

Prior to the production simulation of the pointer object, we performed a set of exploratory MD simulations to determine the solution conditions that correspond to the experimental situation (25). The number of ions in the bulk (>5 nm from DNA atoms) and the volume of the simulation system were iteratively adjusted until the system's pressure reached ∼1 bar and the bulk concentrations of Na+, Cl− and Mg2 + ions were 5, 45 and 20 mM, respectively, which were the conditions realized in experiment (25). The final system contained 7954 magnesium–hexahydrate complexes, 162 Na+ and 994 Cl− ions and 1 709 144 water molecules and measured 307.03 × 350.05 × 521.9 Å3.

To build the system for the production simulation, we merged the idealized model of the pointer structure with water and ions in amounts determined by the exploratory simulations. The resulting system was energy minimized and equilibrated in the NVT ensemble; the system's dimensions were set to the above values. The system was first simulated for 3 ns having all non-hydrogen atoms of DNA restrained to their initial coordinates; the spring constant of each restraint k = 1 kcal/mol Å2. Next, the harmonic restraints to Cartesian coordinates were replaced by an elastic network of pairwise-defined intra-helical restraints. Specifically, harmonic springs (k = 0.1 kcal/mol Å2) were placed between all intra-helical pairs of non-hydrogen DNA atoms located within 0.5 nm of one another; the rest length of each harmonic spring restrain was taken to be the distance between the atoms in the initial (ideal) atomic structure. The use of elastic network restraints allowed the global conformation of the DNA origami object to relax while preserving the local base pairing and base stacking patterns. The structure was simulated with the elastic network restraints for ∼15 ns. The production simulation of the final system was carried out in the absence of any restraints for 210 ns.

Elastic network-guided simulations

As a starting configuration to elastic network-guided simulations, we used the idealized coordinates of the pointer object. Water and ions were not included in the simulation system. During the elastic network-guided simulation, the intra-helical elastic network restraints used in the explicit solvent simulation were retained. No attempt was made to optimize the rest lengths or spring constants of the elastic network restraints. Electrostatic interactions were completely truncated at 10 Å; the PME calculation for long-ranged electrostatics interactions was not performed. To mimic inter-helical electrostatic repulsion, harmonic bonds (k = 1 kcal/mol Å2) with a 31 Å rest length were placed between atoms in adjacent helices. Specifically, these bonds were placed between pairs of phosphorous atoms for every pair of nucleotides having the same 3′-to-5′ orientation (so that a nucleotide of a staple strand was connected to a nucleotide of the scaffold strand and vice versa) located in the same plane of the structure, at least 10 bp away from any crossover between the helices.

Explicit solvent simulation restarted from elastic network-guided structure

The all-atom structure obtained after 1.7 ns elastic network-guided simulation was submerged with solvent as described above. Following a brief minimization, the resulting structure was equilibrated for 43 ns in the absence of the elastic network restraints; weak restraints (k = 0.1 kcal/mol Å2) were applied to maintain the Watson–Crick pairing. Finally, the system was equilibrated for another 74 ns in the absence of any restraints.

RESULTS

Following a previously described protocol (26), we built an idealized model of the pointer structure from its caDNAno design and submerged the resulting all-atom model in water containing magnesium ions at experimental concentration. The resulting model of 5 767 572 atoms (shown in the inset of Figure 2B) was equilibrated for 15 ns in the presence of restraints that preserved the base pairing and base stacking pattern of the structure, followed by ∼200 ns unrestrained simulation. At the beginning of the simulation, the DNA origami structure quickly departed from its initial idealized configuration of parallel helices toward a conformation characterized by a global right-hand twist, consistent with the cryo-EM reconstruction, Figure 2A. The helices, which were tightly packed at the beginning of the simulation, spread into the characteristic chickenwire configuration. The root mean square deviation (RMSD) of the simulated pointer structure from the cryo-EM derived pseudo-atomic model gradually reduced to a value similar to the reported resolution of the cryo-EM reconstruction, Figure 2B, indicating a good agreement between the simulated and experimentally derived structures.

Figure 2.

Explicit solvent all-atom MD simulation of the pointer object. (A) Instantaneous conformations of the pointer object during a 200-ns MD trajectory. The DNA atoms are depicted as blue and white vdW spheres. Single-stranded DNA at the ends of the helices are not shown for clarity. This simulation trajectory is illustrated by Supplementary Movies 1–3. (B) Root mean square displacement (RMSD) of the pointer object's coordinates with respect to the pseudo-atomic model obtained from cryo-EM reconstruction (25). The RMSD was calculated using non-hydrogen atoms of the scaffold strand. The gray region indicates the part of the MD trajectory where harmonic or elastic restraints were applied to DNA. The inset illustrates the full extend of the explicit solvent model. (C) The dimensions of the pointer object relative to the cryo-EM derived model. The dimensions of the pointer structure along the x- and y-axes were defined as the center of mass (CoM) distances between the three center helices at the opposite faces of the pointer structure. The size along the z-axis was defined as the CoM distance between two terminal 10-bp slices of the same nine center helices (at the opposite faces of the structure). The background images define the x, y and z dimensions using a pointer structure. (D) The inter-helical distance map of the pointer structure. The color of each helix indicates its average separation from the neighboring helices. The top and bottom rows illustrate the inter-helical distances computed along the x and y dimensions of the structure, respectively. The left and right columns illustrate the inter-helical distances computed for the last frame of the explicit solvent MD trajectory and the cryo-EM derived structure, respectively. The inter-helical distances were calculated by averaging the CoM distances between adjacent 10-bp segments within the center 40 bp of the pointer object.

For further quantitative comparison, we plot the dimensions of the simulated pointer structure versus simulation time, Figure 2C. The pointer's dimension along the DNA helical axis (the z axis) quickly relaxed to the value of the cryo-EM derived model. In the x − y plane, the structure expanded gradually throughout the simulation reaching 95% (x axis) and 91% (y axis) of the dimensions of the cryo-EM derived model. To elucidate the origin of the asymmetric expansion, we calculated the average inter-helical spacing of neighboring DNA helices along the x and y axes independently, for both simulated and experimentally derived models, Figure 2D. Whereas the simulated model exhibits a fairly uniform inter-helical spacing along the x and y axes, the separation of the helices in the cryo-EM derived model is a few percent larger, on average, along the y axis.

One possible explanation for the asymmetric inter-helical spacing in the experimentally derived structure is the surface interaction of the structure with the carbon substrate used to obtain the cryo-EM images (25). Because of its asymmetric shape, the contact area of the pointer object with the substrate surface depends on the orientation of the object. When the y − z plane of the pointer object, rather than the x − z plane, lays against the surface, less of the structure extends into the solution away from the surface. Therefore, if surface interactions were to contribute to spreading of the helices, one would expect the object to appear more elongated along the y axis than along the x axis. Alternatively, the asymmetric inter-helical spacing could reflect the actual asymmetry of the pointer object that is not captured by the MD method. Indeed, the pointer structure is shorter along the z axis on the −y face than on the +y face. Furthermore, there are only three locations (except those proximal to the ends of the helices) where consecutive crossovers span more than three helices along the y axis whereas eight such motifs exist along the x axis, including two that span six helices. However, a Pearson correlation analysis of the center 40 bp of the average cryo-EM derived model revealed a rather modest correlation (correlation coefficient of −0.111) between the inter-helical distance and the presence of consecutive crossovers between adjacent 10-bp fragments. The same analysis performed over the conformations sampled during the last 48 ns of the MD simulation yielded a virtually identical correlation coefficient (−0.114). Thus, although consecutive crossovers can indeed locally compress the structure, they do not appear to be responsible for the compression of the pointer structure along the x axis in the cryo-EM derived model.

The good overall agreement between the simulated and experimentally-derived structures of the pointer object suggests that all-atom MD simulations could be used for testing and prototyping of DNA origami designs. However, explicit solvent simulations of a typical DNA origami object are still too demanding, computationally, to be practical for routine characterization of DNA origami objects. Hence we sought an alternative approach that would be less computationally costly but still could provide a reasonably accurate atomistic model.

Elimination of explicit solvent combined with a custom elastic network of restraints reduced the computational cost of DNA origami simulation by a factor of 10 000 without compromising accuracy of the average structure characterization. In these elastic network-guided simulations, all non-bonded interactions between DNA atoms were limited in range to 10 Å using the standard 8–10 Å smooth cutoff scheme (28). In the absence of the long-range electrostatic forces that results in inter-DNA repulsion, custom elastic restraints were necessary to produce the characteristic chickenwire-like organization of DNA helices observed in our previous MD simulations (26) and experiment (1,25). Specifically, harmonic bonds with a spring constant of 1 kcal/(mol Å2) and the rest length of 31 Å were placed between the phosphorous atoms of the neighboring helices for each pair of base pairs located within the same plane of the structure and at least 10 bp away from any crossover between the helices, see inset of Figure 3A for a schematic representation. In addition to restraints that enforced the chickenwire pattern, the local base pairing and base stacking patterns were preserved using a network of harmonic bonds (36), same as the network employed at the beginning of each explicit solvent simulation (26).

Figure 3.

Elastic network-guided simulation of the pointer object. (A) Chickenwire representation of the DNA origami construct before (left) and after (right) a 1.7-ns elastic network-guided simulation. During the simulation, intra-helical elastic network restraints maintained basepairing and basestacking (represented with atomic detail in the lower magnified region), while inter-helical restraints (represented schematically in the upper magnified region) maintained the prescribed inter-helical distances. This simulation trajectory is illustrated by Supplementary Movies 4–6. (B) RMSD of the pointer object's coordinates with respect to the cryo-EM derived model during the elastic network-guided simulation. The inset shows the RMSD of the pointer object during explicit solvent simulations started using the configuration obtained at the end of the elastic network-guided run (blue) and the idealized conformation (red, same as in Figure 2B). (C) The dimensions of the simulated DNA origami object relative to the cryo-EM derived model. The object's dimensions along x, y and z axes were calculated as in Figure 2C. Each plot depicts the dimensions along one axis during explicit solvent simulations started using the configuration obtained at the end of the elastic network-guided run (thick lines) and the idealized conformation (thin lines, same as in Figure 2C). (D) The inter-helical distance map of the pointer structure. The color of each helix indicates its average separation from the neighboring helices. The maps were computed as described in the caption to Figure 2D. The simulated map characterizes the state obtained at the end of the elastic network-guided simulation. (E–G) Comparison between the simulated (blue) and cryo-EM derived (red) structures of the pointer object. The simulated structure characterizes the state obtained at the end of the elastic network-guided simulation. A movie comparing the structures is available in Supplementary Data (Movie 7). (E) The entire pointer object, viewed from the y axis. (F) A slice of the object showing the characteristic chickenwire pattern. (G) Spreading of the helices due to crossover omission in the design of the pointer structure. In panels E and F, the simulated and cryo-EM derived models were aligned according to the coordinates of the scaffold strands. (H) A scheme for characterizing a crossover geometry (25). Each crossover consists of four ‘legs’ that diverge from a central connection. The line connecting the CoM of basepairs within each helix at the junction defines the crossover vector x, depicted in dark gray. The four legs are associated with vectors a–d that are each given by the line connecting the CoM of one of the base pairs at the junction with the base pair located two base pairs away from the junction within the same duplex. These vectors are used to calculate the away-from-plane bend β. The vectors are then projected into the plane defined by x before the leg-to-leg inclination angle γ and in-plane helix bend, α, are obtained. (I) Distributions of the crossover angles. The crossover angles, defined in panel H, were calculated for the cryo-EM derived model (pink) and for the model derived by averaging the MD trajectory (blue). The distributions of the in-plane and away-from-plane bends computed using instantaneous snapshots of the MD trajectory (gray) are significantly broader because of the thermal fluctuations.

Starting from the idealized all-atom model, 2 ns of elastic network-guided MD simulation brought the pointer structure into a stable conformation characterized by RMSD from the cryo-EM derived structure of <1 nm, Figure 3B, which is below the resolution of the cryo-EM density map (25). To test if the structure produced by the elastic network-guided simulation is compatible with the in situ environment of DNA origami, we merged the elastic network-guided structure with solvent and simulated the resulting explicit solvent model for ∼130 ns. During the explicit solvent simulation, the RMSD of the structure increased slightly due to larger thermal fluctuations, see inset of Figure 3B, but remained considerably lower than the final RMSD of the structure simulated in explicit solvent starting from the idealized conformation, indicating that the latter simulation did not reach a steady state.

The dimensions of the pointer object in the elastic network-guided simulation rapidly reached the dimensions of the cryo-EM derived structure along the x and z axes, but remained ∼5% shorter along the y axis, Figure 3C; the pointer structure maintained its dimensions during the subsequent explicit solvent simulation. The average inter-helical spacing of the elastic network-guided model was very close to 2.5 nm along both x and y directions, Figure 3D, in contrast to the cryo-EM derived model that showed larger inter-helical distance along the y axis only.

Overall, the local features of the pointer structure produced by our elastic network-guided method are in remarkable agreement with the features of the cryo-EM derived model, see Figure 3E–G. For quantitative analysis of the local conformations of crossovers, atomic coordinates of the pointer structure were averaged over the last 19.2 ns of the explicit solvent simulation that began from the structure determined by the elastic network-guided method; the RMSD of the ‘average structure’ from the cryo-EM derived model was only 0.92 nm. Overall, the distributions of the crossover angles in the average MD and cryo-EM derived structures are similar, Figure 3I, except that the in-plane helix bend α was more open in the simulation structure whereas the away-from-plane bend was more closed. When the crossover angle distributions were calculated for the ensemble of conformations observed during the explicit solvent MD simulation, much broader distributions were obtained for the in-plane and away-from-plane bends due to thermal fluctuations whereas the leg-to-leg inclination was not affected by them, Figure 3I. The more compact distribution of the crossover angles in the cryo-EM derived model is not surprising because thermal fluctuations were implicitly averaged during the electron density reconstruction process.

The pointer object also included a number of unusual structural motifs. Consecutive crossovers span up to six helices causing local compaction. These crossovers were well-resolved in the cryo-EM reconstruction, indicating overall good structural stability. The overall configuration of each set of consecutive crossovers was similar in models derived from simulation and cryo-EM reconstruction. However, rupture of several basepairs was observed during the MD simulations, indicating the presence of substantial stress in these motifs. Similarly, three basepairs ruptured during explicit solvent MD simulation in another unusual motif—a left-handed ‘psuedohelix’ formed by a strand weaving through the pointer. Another notable motif was the omission of a single basepair at the edge of the pointer object. In the cryo-EM derived model, the helix of the missing basepair was bent strongly toward a neighboring helix and its respective crossover was distorted. In the simulated models, the two helices were bent by similar amounts with a left-handed twist. The basepair forming the crossover in the shorter helix ruptured after 40 ns of explicit solvent simulation that began from the idealized coordinates but was intact throughout the explicit solvent simulation that began from the state obtained at the end of the elastic network-guided simulation.

DISCUSSION

We have shown that the explicit solvent MD method can produce realistic atomic-scale structural models of large DNA origami objects with accuracy of the cryo-EM reconstruction technique. We have developed an elastic network-guided simulation protocol that dramatically reduces the computational cost of the structure prediction, making it suitable for routine use in design and development of DNA origami structures. Our structure prediction method is complementary to continuum mechanics (37) and coarse-grained MD (38) approaches in providing microscopically accurate information about the local structure of self-assembled DNA systems. A web service that provides all configuration files required to run our elastic-network guided simulation of an arbitrary DNA origami object is available at http://bionano.physics.illinois.edu/origami-structure.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online.