A coarse-grained model for DNA origami.

Modeling tools provide a valuable support for DNA origami design. However, current solutions have limited application for conformational analysis of the designs. In this work we present a tool for a thorough study of DNA origami structure and dynamics. The tool is based on a novel coarse-grained model dedicated to geometry optimization and conformational analysis of DNA origami. We explored the ability of the model to predict dynamic behavior, global shapes, and fine details of two single-layer systems designed in hexagonal and square lattices using atomic force microscopy, Förster resonance energy transfer spectroscopy, and all-atom molecular dynamic simulations for validation of the results. We also examined the performance of the model for multilayer systems by simulation of DNA origami with published cryo-electron microscopy and atomic force microscopy structures. A good agreement between the simulated and experimental data makes the model suitable for conformational analysis of DNA origami objects. The tool is available at http://vsb.fbb.msu.ru/cosm as a web-service and as a standalone version.

INTRODUCTION

DNA origami technique allows to construct 2D and 3D objects from DNA strands (1). This is achieved by folding a long ‘scaffold’ chain into a designed shape by site-specific complementary interactions with a custom set of short ‘staple’ strands (2). The scaffold chain has a typical length ranging from hundreds to thousands of nucleotides and often has a natural origin. Up to a few hundreds of staple strands with a typical length of 40 nucleotides are required to fold the scaffold chain properly. Every staple oligonucleotide binds several distant parts of the scaffold DNA, forming double-helical tracts and bringing them together. The individual tracts stitched with interhelical connections arrange in space, filling a dense lattice. The connections, or so-called crossovers, have a form of antiparallel Holliday four-way junction (3).

There are two major types of lattice: honeycomb and square. In the honeycomb lattice the DNA helices are located in the nodes of a hexagonal network (Figure 1A). Consequently, every helix has up to three potential neighbors separated by an angle of 120°. The lattice geometry results from the crossovers between the helix and its neighbors occuring with an interval of seven base pairs (Figure 1B, Supplementary Figure S1A), which corresponds to a B-DNA twist angle of 240° (Figure 1A and C). In the square lattice, a helix has up to four potential neighbors separated by an angle of 90°, and the crossover periodicity is eight base pairs, or 270° (Supplementary Figure S1B).

Figure 1.

The DNA origami fundamentals: lattice and crossovers. (A) The helices (black numbered circles) are located in the nodes of the honeycomb lattice (gray). Colored lines between the circles represent crossovers between the helices and their colors correspond to individual staples from the panel B. The colored sectors within the circles refer to the twist of the staples on their path between the helices shown with a dashed line in the panel B; (B) a side view of the design from the panel A, numbers of the helices correspond between the panels. The helices are formed by a scaffold chain (black) and a number of staple chains (colored). The dashed line shows the staple-formed path from helix 1 to helix 5 corresponding to the colored sectors and crossovers from the panel A; (C) a full-atom view of the area highlighted in gray on the panel B, colors of the DNA chains correspond between the panels.

There are computer-aided design tools for the DNA origami engineering (4,5); the most popular tool is caDNAno (6). It is well suited for 2D objects drafting, but designing of 3D systems with it is challenging. This issue was addressed in a few recent studies (7,8). The main difficulty is the complexity of interconnections between the individual helices in a non-planar DNA origami, which often looks confusing on a 2D plot. There are dedicated software tools for verification of 3D DNA origami designs, such as CanDo (9–11). CanDo utilizes a simplified model of DNA optimized for the DNA origami simulation intended to predict an equilibrium shape of the origami and estimate its mobility. CanDo is fast and provides a valuable structural feedback. Unfortunately, the DNA model used in CanDo does not consider non-covalent interactions (electrostatic and Lennard-Jones terms), while it was shown by Arbona and coauthors that electrostatic interactions are necessary to reproduce experimental structural data of DNA origami (12). It makes the tool not applicable for simulation of dynamic structures (Supplementary Figure S2) and also could lead to the incorrect shape prediction for some complex designs (Supplementary Figure S3).

These problems could be solved with an employment of a more sophisticated physical model of DNA. Methods for generation of a representative equilibrium ensemble, such as molecular dynamics (MD) or Monte-Carlo simulations, provide a powerful platform for realization of such model as a specialized force field. The size of DNA origami objects (millions of atoms for explicit solvent systems) dictates the choice of coarse-grained (CG) force fields for the task. Modern CG DNA models, such as Martini (13), SIRAH (14), oxDNA (15), Knotts model (16) and others (12,17,18) can reliably reproduce the physical and geometric properties of B-DNA, and there are specialized tools which can handle caDNAno output (15,19). However, geometry optimization for systems that are far from an equilibrium state (for example, 3D objects designed in a 2D space) remains a challenge. Optimization procedure usually converges in a local minimum, which depends on the initial coordinates and therefore could significantly differ from the real shape. There are several approaches to overcome this issue (20,21), reviewed by Valsson and coauthors (22).

The existing force fields provide a wide range of CG mapping approaches, with detalization varying from several particles to one particle per nucleotide. Unfortunately, DNA origami systems remain computationally expensive. We suggest that in this case the extreme low-resolution approach is justified, taking into account the size of origami systems and a large number of interconnections that define their shape.

In this paper we suggest a new coarse-grained model of the DNA origami. Parameters of the model were verified with atomic force microscopy (AFM), Förster resonance energy transfer (FRET) spectroscopy, and all-atom MD simulation. At the expense of atomistic and single-nucleotide level details our model is specifically tuned for geometry optimization and efficient scanning of a conformational landscape of DNA origami designs. The approximations underlying the model design result in a number of limitations that should be considered assessing the accuracy of simulations: (i) single-stranded regions are the most challenging objects in DNA modeling (23); due to lack of structural experimental data we were not able to estimate the reliability of reproduction of single-stranded origami substructures; (ii) nucleotide sequence details are neglected; (iii) the model does not consider the influence of temperature and salt concentration on the DNA origami self-assembly efficiency; (iv) the correspondence between the physical and simulated time is not straightforward; our model significantly speeds up the large conformational changes.

We also present a tool called COSM (Coarse-grained Origami Structures Modeling) based on this model dedicated for the DNA origami shape prediction and conformational analysis. The tool is available at http://vsb.fbb.msu.ru/cosm as a web-service and as a standalone version.

MATERIALS AND METHODS

CG modeling

GROMACS 4.6.5 software package (24) was used for geometry optimization and MD simulation. Two successive steps of energy minimization were applied for caDNAno-derived coordinates of an origami design to remove strains. First, the quasi-Newtonian limited memory BFGS algorithm (25) was used with cutoff distance of 23 Å for non-covalent interactions. Second, energy minimization using steepest descent algorithm (all cutoff values equal 23 Å) was performed. After that, MD simulation was carried out at T = 300 K with a time constant for coupling of 1.0 ps under the control of a velocity rescaling thermostat (26), all cutoff values for non-covalent interactions were 15 Å. The time step for integration was 100 fs, 900 000 steps in total. Coordinates were written to an output trajectory every 2 ps. We used the oxDNA 2.2.2 package (15) for a complementary CG modeling of the Hc-system. The simulation parameters used are given in the Supplementary Note S2.3.

All-atom modeling

The starting coordinates of the component DNA helices were obtained with the 3DNA v.2.3 software package (27). The helices were arranged in 3D with the inter-strand distance of 24 Å and, in case of the Hc-system, angle of 240° using PyMOL Molecular Graphics System, v. 1.8.6.0 (28). After that the strands were rotated to comply the chain flow at the crossover sites using the PyMOL software; the nucleotides order in the resulting coordinates file was edited to provide sequence integrity of the individual strands. GROMACS 2016.3 software package was used for the all-atom simulations. Explicit solvent simulations employing SIRAH water model (29) in parmbsc0 (χ OL4) force field (30) were performed at T = 300 K under control of velocity rescaling thermostat (26). Cutoff values for non-covalent interactions were 12 Å. To minimize the influence of the neighboring structures through the periodic boundary conditions, a triclinic box was added around the DNA; a distance to the box walls was 200 Å. A 2 fs time step was used for the simulations with 5 × 107 steps in total. The simulations were performed on the Lomonosov-2 supercomputer at the Research Computing Center of Moscow State University.

Synthesis of the scaffold DNA

A coding sequence of human glial cell-derived neurotropic factor (UniProt ID P39905, isoform 1) was used as a source of a linear scaffold DNA. The sequence was synthetised as described elsewhere (31). The desired gene fragment was extracted with polymerase chain reaction (PCR) using oligonucleotides 5′-d(CCACCATGTCACCAGATAAACAA)-3′ and 5′-d( TGGATCCCAGATACATCCACACCTTTTAGCGG)-3′ as a forward and reverse primers, correspondingly. Here and elsewhere, PCR was carried out using the Taq system (Evrogen) which included 10 mM (2.5 mM each) deoxynucleoside triphosphates (dNTPs), 10× Taq Buffer and 5 U/μl Taq DNA polymerase. Thermal cycling consisted of 94°C for 1.5 min, followed by 30 cycles of 94°C for 15 s, 59.1°C for 20 s and 72°C for 40 s, with final extension step of 72°C for 10 min. PCR product was separated with 1.5% agarose gel electrophoresis using EtBr for visualization and extracted with QIAquick Gel Extraction Kit (Qiagen). Resulting DNA template has a total length of 418 bp. The sequence is provided in Supplementary Note S2.1.

A single-stranded DNA was obtained from the template DNA with asymmetric PCR (aPCR) using the same primers set with concentrations of 250 nM for the forward primer, 5 nM for the reverse primer and 2.5 nM for the template DNA. Amplification protocol consisted of 94°C for 5 min, followed by 35 cycles of 94°C for 30 s, 60°C for 20 s and 72°C for 25 s, with a final extension step of 72°C for 1 min. Gel mobility of a single-stranded DNA turned out to be highly dependent of environmental conditions, therefore aPCR products separation was performed in 1.5% agarose gel in an ice bath using Sybr Gold (Invitrogen) for visualization. The single-stranded DNA was extracted from gel using QIAquick Gel Extraction Kit (Qiagen).

DNA origami assembling

The systems consisted of 10 nM of scaffold DNA and 200 nM of each staple oligonucleotide in TE-buffer (5 mM Tris–HCl, 1 mM EDTA) supplemented with 5 mM NaCl and 10 mM MgCl2. The mixture was cooled from 95°C to 20°C with a step of 1°C/min using a T100 PCR thermal cycler (Bio-Rad) to anneal the strands. Folded systems were isolated by three rounds of centrifugal filtration using Amicon Ultra-0.5 ml 30 K centrifugal filters (Millipore). On each round 1× volume of folding buffer was added to the sample to bring a volume to 500 μl and it was then centrifuged for 10 min at 14 000 rpm (20 000 g) at 12°C. After each round, flow-through was discarded. List of staple oligonucleotides is provided in a Supplementary Table S3.

Sample preparation for AFM

All AFM experiments were performed on freshly cleaved mica. Before deposition on the substrate, the origami solution was diluted five times with a buffer containing 5 mM Tris–HCl and 6 mM MgCl2. A drop of the diluted solution was applied on the mica surface for 5 min and then removed from the surface with a flow of nitrogen. This sample preparation method eliminates salts and also minimizes both aggregation and denaturing artifacts during absorption and drying.

AFM measurements and image processing

The AFM experiments were performed using a multimode AFM with Ntegra Prima controller (NT-Mdt, Russia) in tapping mode with a typical scan rate of 1 Hz and typical free amplitude of several nm. All measurements were performed in air using supersharp silicon cantilevers (32) with a tip diameter ∼1 nm.

Standard algorithms of AFM images flattening were used (subtraction of quadric surface and averaging by lines), no algorithms of resolution improvement were used. FemtoScan Online software (ATC, Russia) was used to filter, analyze, and present the AFM data.

FRET measurements

Fluorescence-labeled oligonucleotide st4L: FAM-5′-d(ATCCCAGATACATCACACCTTTGTTTATCTGGTGACATGG)-3′-TAMRA was purchased from Synthol, Russia. To assemble labeled Sq-systems we used 10 nM of scaffol DNA and labeled oligonucleotide and 200 nM of each of the remaining staple strands.

Annealing was done as described above, no filtration was performed. For every system from 100 to 140 data points were recorded using Infinite M200 PRO microplate reader (Tecan). To prevent excitation light from the flash lamp reaching into the detector, a distance of 34 nm between the excitation and emission wavelengths was kept. For descriptions of the FRET data processing see Supplementary Note S2.4.

Statistical analysis

P-values were determined by a Student’s unpaired t-test using R 3.0.2 (33). The significance threshold was set at 0.05.

RESULTS AND DISCUSSION

The main purpose of our CG model was to reliably describe geometry and dynamics of origami designs using effectively simplified representation of DNA origami. In our model, DNA origami design is represented by a worm-like chain (34) of particles corresponding to double- and single-stranded DNA regions. It traces the path of the scaffold strand, i.e. there are no explicit staples. The CG particles are connected through bonds realized as harmonic potentials and influence each other by electrostatic and van der Waals interactions.

Types of particles used in the CG model

Conversion of a DNA origami design to the CG model starts with identification of the boundaries of double- and single-stranded regions (see Supporting Note S1.1, which covers input data processing). After that, the regions are filled with CG particles (here and below refer to Figure 2 illustrating the model and Table 1 listing parameters of the CG particles).

Figure 2.

Transition from initial 2D scheme to the COSM model. Blue-encircled numbers specify the index numbers of the corresponding strands. (A) 2D-plot of the DNA origami. Black: scaffold chain; gray-blue: staple chains. Each cell corresponds to one base pair. Sites of insertions are denoted with a tear-shaped sign with a plus symbol; (B) the design from the panel A schematically shown in its COSM representation. The scaffold chain is shown by a black line. Crossover-forming particles are outlined by gray; (C) correspondence between the COSM model for the strand 1 and a B-DNA geometry. COSM particles are in gray, staple crossover site is shown with arrows; (D) an actual geometry of the encircled area from the panel B.

Table 1.

Particles of the model

Particle	Sizea (σii), Å	Charge, e		Molecular mass, g/mol		Bond length, Å (base pairs)
		Honeycomb	Square	Honeycomb	Square	Bond type	Honeycomb	Square
S	4.4	−0.111	−0.111	340	340	S−S, S−T	6.3 (1 nt)	6.3 (1 nt)
H	23.0	−0.172	−0.196	4760	5440	H−H	23.8 (7)	27.2 (8)
						T−H, B−H, PT−H	23.8 (7)	27.2 (8)
T, B, PT	4.4	−0.172	−0.196	2380	2720	T*−T*b	22c	22c
						B*−B*b	3.4 (1)	3.4 (1)
T1, B1	4.4	−0.025	−0.025	680	680	T1−PT, B1−PT	3.4 (1)	3.4 (1)
T2, B2	4.4	−0.049	−0.049	1360	1360	T2−PT, B2−PT	6.8 (2)	6.8 (2)
T3, B3	4.4	−0.074	−0.074	2040	2040	T3−PT, B3−PT	10.2 (3)	10.2 (3)
T4, B4	4.4	−0.098	−0.098	2720	2720	T4−PT, B4−PT	13.6 (4)	13.6 (4)
T5, B5	4.4	−0.123	−0.123	3400	3400	T5−PT, B5−PT	17.0 (5)	17.0 (5)
T6, B6	4.4	−0.148	−0.148	4080	4080	T6−PT, B6−PT	20.4 (6)	20.4 (6)
T7, B7	4.4	—	−0.172	—	4760	T7−PT, B7−PT	—	23.8 (7)
N	4.4	−0.025	−0.025	680	680	N−N, N−PT	3.4 (1)	3.4 (1)

aSize of the particles in our model is equal to the B–DNA diameter or the distance between base pairs increased by 1 Å to avoid intersections and direct contacts between nucleotides of corresponding all-atom models. bThe asterisk substitutes empty symbol or digit from 1 to 7. cB-DNA diameter.

Single-stranded regions are filled with S particles. Each S particle describes an individual nucleotide and has corresponding size, mass, charge and distance to the neighboring S particle.

Every double-stranded region begins and ends with a terminal particle, the space between the terminal particles is filled with helical particles H. The H particle corresponds to seven base pairs in the honeycomb lattice and has a corresponding mass and charge (corrected for charge shielding as described in Supporting Note S1.3.2). A bond length between the double-stranded particles is equal to the length of seven base pairs. As a result, each particle locates in the node of the staple-crossover network. The H particle has a diameter (excluded volume) of a DNA duplex to model interhelical interactions.

Terminal particles are smaller than the H particles to neighbor other DNA regions without radii overlap. There are two types of terminal particles: when two helices are connected through a staple crossover, the type T is used (Figure 2, strands 1 and 2), otherwise, the adjacent termini have the B type particles (strands 2 and 3). The T and B particles have similar physical properties: each particle represents 3.5 bp by molecular weight (this is needed for an even mass distribution along the helix), but has a size (excluded volume) of only 1 nt. The difference between these particles is in the bond length: T–T bond length is 22 Å to avoid bending of the helices at the ends, while the B–B bond has a length of the helical rise, corresponding to the actual distance between the ends of the helices (Figure 2D).

Often the length of a double-stranded region is not a multiple of seven base pairs. For an accurate description of such helices’ length we use subtypes of the terminal particles. For example, if a helix has one ‘excess’ base pair (right side of the strand 1), we place a T1/B1 particle at a helical rise distance from the last ‘nodal’ particle. For five excess base pairs we use particles T5/B5 (left side of the strand 1) and so on. To avoid radii overlap, we substitute the last H particle in the helix with a small ‘pre-terminal’ particle PT. It has the same physical properties as the particles T and B, but behaves as the H particle in terms of covalent interactions (see Supplementary Note S1.2 describing the types and parameters of covalent interactions).

DNA origami technique allows to design twisted or curved structures. This is achieved by insertion or deletion of base pairs inside individual helices on specified positions (central part of the strand 3). To model that, two H particles flanking the insertion/deletion site are replaced by the 7 bp-long segment -PT-(N)6-PT-. Here, N is a particle corresponding to a single base pair; the bonds N-N and N-PT have a length of a helical rise. Then a desired amount of the N particles is inserted or deleted within the segment. As a result, an exact length of the segment could be specified using only one additional particle type.

We use MD simulation method to predict conformational landscape of a DNA origami design. Covalent interactions between the CG particles are described using harmonic potentials with force constants derived from physical properties of double- and single-stranded DNA (Supplementary Note S1.2). Non-bonded particles can influence each other via Coulomb and Lennard-Jones interactions (Supplementary Note S1.3). In the same manner as staple oligonucleotides guide folding of a scaffold DNA, final shape of the CG strand in our model is controlled by restraints of the distances between individual particles (Supplementary Note S1.4). The interactions between the particles could be modified in accordance to the ionic strength of the system (Supplementary Note S1.5).

Benchmark systems

Two single-layer benchmark systems were designed using caDNAno software (6) and employed to validate our model (Figure 3). One system was designed in a honeycomb lattice and we named it ‘Hc-system’. Another was designed in a square lattice and named ‘Sq-system’. Both systems were of small size (scaffold length of 418 nucleotides) in order to verify structures predicted by the model with all-atom MD simulations, AFM and FRET.

Figure 3.

Designs of the benchmark structures. (A) scheme of the Hc-system; (B) scheme of the Sq-system with putative arrangement of the arms.

Hc-system

The Hc-system consisted of a ‘body’ and two ‘legs’ which were parallel to each other in the initial caDNAno scheme (Figure 3A, Supplementary Figure S10A). However, AFM imaging has shown that the real structure deviated from this scheme significantly (Figure 4A). According to AFM, the Hc-system had a V-like shape with an angle of 64.7 ± 25.6 degrees (N = 30) between the legs (Figure 4B). The height of the Hc-system’s body was 18 ± 4 Å, and the height of the legs was 11 ± 3 Å. Both values were higher than a typical AFM-measured height of a double-stranded DNA on mica (35). The difference between the two values suggests that the system was not planar, which implies dynamic behavior of both the legs and the body.

Figure 4.

Analysis of experimental and simulated shapes of the Hc-system. (A) AFM images of individual structures; (B) conformations of the Hc-system predicted by our model in side and front views; (C) scheme of the characteristic movements of the Hc-system and corresponding torsional angles: γ for swinging and θ for twisting motions of the legs; (D) probability distribution of the γ angle values according to our model (red) and AFM imaging (blue); (E) Dynamics of the θ angle according to coarse-grained (red) and all-atom (black) models; (F) mobility of atoms/particles of all-atom (black; only phosphorus atoms of the scaffold chain were considered for the RMSF calculation) and coarse-grained (red) models during MD simulation.

Results of conformational analysis of the Hc-system with our COSM tool were in a good agreement with the experimental data. According to the model, the system had a V-like conformation (Figure 4B) formed as a result of swinging movements of the legs. The legs movements were accompanied by twisting of the body, highlighting dynamic behavior of all parts of the system. The almost two-fold difference in the body and legs thickness in the AFM data may correspond to the body twist, which in this case correlates with the modeling.

A torsional angle γ between the legs (Figure 4C) was 64.41 ± 15.39° (Figure 4D); there was no statistically significant difference between the simulated and experimentally measured angles (independent samples t test: NMD = 2000, NAFM = 30, t = 0.1012, df = 2028, P = 0.82. Here and elsewhere: NMD, NAFM are the sample sizes obtained from MD and AFM data, correspondingly). According to a 90 ns-long MD trajectory, the terminal values of the angle were 18.02° and 106.8°; a 900 ns-long trajectory provided a slightly better sampling with the values between 15.35° and 133.57°. Thus, the simulated flexibility of the system corresponds to that observed on the AFM images, where the angle between the legs varied from 17° to 123°.

To further validate the COSM model, we performed an all-atom simulation of the Hc-system. Since 10-ns-scale all-atom simulations are capable to overcome free energy barriers of 5–7 kcal/mol only (23), therefore we were not able to obtain an appropriate conformational sampling with the method. We observed moderate deviations around the starting conformation; however, all-atom MD simulations allowed us to identify the main motion vectors of the system and provided the atomistic view on it.

Because of lack of structural data provided by the all-atom modeling, we also analyzed the dynamical behavior of the Hc-system using the oxDNA CG model (15). oxDNA and all-atom MD simulations agree that the legs indeed performed swinging motions directed both in- and out-of-body-plane (Supplementary Notes S2.2, S2.3 and Supplementary Figures S11A and S13). The motions were accompanied by twisting along the legs’ axes (Figure 4C, angle θ, Supplementary Figure S13), which we also observed in the COSM modeling of the system. The mean value of the twist angle according to the all-atom trajectory was 32.61 ± 10.27° (Figure 4E). There was no statistically significant difference with the value of 29.72 ± 19.50° obtained with our computational model (independent samples t test: Nall-atom = 200, NCOSM = 200, t = 1.8545, df = 398, P = 0.064; Nall-atom, NCOSM are the sample sizes randomly obtained from the all-atom and COSM MD trajectories, correspondingly).

The relative mobility of different parts of the system, expressed in root mean square fluctuations (RMSF) of the system’s particles, correlated well between the coarse-grained and all-atom modeling (Figure 4F). The only exception was the 5′-terminal part of the Hc-system, which showed high mobility in our model, while being stationary in the all-atom simulation. We believe that this is due to the insufficient sampling provided by the all-atom trajectory. The Hc-system is obviously mobile according to the AFM and coarse-grained modeling results (both COSM and oxDNA), but the largest movement observed throughout the 100 ns-long all-atom trajectory had an amplitude of 32.9 Å, while the amplitude of the same process was 250.2 Å in our computational model and 307.0 Å in oxDNA. It might attribute to the fact that the coarse-grained models greatly speed up the conformational sampling (Supplementary Figure S12 and Supplementary Notes S2.2.1, S2.3). Significantly longer all-atom MD simulations are needed to observe the large-scale conformational changes, which is very time- and computationally consuming. Note that larger trajectory time scale may lead to accumulation of force field parameterization errors corrupting the B-DNA geometry (36).

Our model allowed to reduce the Hc-system size 245 times, from 25 524 atoms (287 842 atoms with solvent) in the all-atom representation to 104 particles in the coarse-grained form. Processing time for the MD modeling part of the conformational analysis workflow was 29 seconds on a single core of Intel Core i7-6700K processor, the file size of the resulting trajectory was 60MB. The file size of the 100 ns-long all-atom trajectory was 70GB, its calculation took 4800 min of processing time on 224 Intel Xeon E5-2697 CPU cores and 16 NVIDIA Tesla K40s GPUs. The oxDNA model of the Hc-system consisted from 1608 atoms; the 5 × 107 steps-long MD trajectory calculation took 668 minutes on a single core of Intel Core i7-6700K processor, the trajectory file size was 631MB. Despite its simplicity, our model provided relevant data on both the global shape and flexibility of the Hc-system, which were supported by the experimental and both oxDNA and all-atom simulation results.

Sq-system

The Sq-system consisted of three arms that concurred at one point (Figure 3B). On the original caDNAno scheme all the arms were parallel to each other (Supplementary Figure S10B). However, according to AFM data, the shape of the Sq-system was a composition of conformations ranging from Y-like to T-like structures (Figure 5A). All three arms had the same height of 13 ± 2 Å. The angles α and β between the arms (Figure 5B) were equal to 108 ± 22° and 121 ± 24°, correspondingly. Crossover periodicity was occasionally seen on the AFM images as regular dark straps separated by ∼100 Å.

Figure 5.

Analysis of experimental and simulated shapes of the Sq-system. (A) AFM images of individual structures; (B) scheme of the Sq-system showing the parameters used for characterization of the system’s geometry: angles α, β and a distance between the ends of the Sq-system; (C) conformations of the Sq-system predicted by our computational model; (D, E) probability distribution of the α and β angle values, correspondingly. Red: the values according to our model, blue: obtained from AFM imaging; (F) mobility of atoms/particles of all-atom (black; only phosphorus atoms of the scaffold chain were considered for the RMSF calculation) and coarse-grained (red) models during MD simulation.

The COSM model of the Sq-system had an Y-like shape with the α angle of 120.3 ± 11.8° and the β angle of 110.0 ± 14.6° (Figure 5C–E). The predicted angle β distribution satisfied the experimental data (NMD = 8000, NAFM = 30, P = 0.092), but the angles α were statistically significantly different from the angles measured with AFM (NMD = 8000, NAFM = 30, P = 0.041). From the angle histograms (Figure 5D, E) it could be seen that on AFM images the Sq-system was more mobile than in the simulation. We performed a 900 ns-long simulations of the Sq-system to improve the sampling. The new values of the angles α = 114.9° ± 13.3° and β = 111.4° ± 14.1° were in a better agreement with the experimental data (P-values were 0.243 and 0.136, correspondingly; NMD = 89 000, NAFM = 30).

According to the initial scheme, crossover periodicity was 32 bp, or 108.8 Å (see Figure 3B, Supplementary Figure S1B). From cryo-EM data (PDB entry 4v5x (37)), helices in a DNA origami with a square lattice are slightly curved, resembling a chickenwire-pattern. This leads to a small decrease in the distance between the crossovers to 107.4 ± 1.3 Å. Our model successfully reproduced this curving pattern (Figure 5C). The distance between crossovers in the model was 106.2 ± 2.4 Å.

The short arm of the Sq-system by design had a possibility of its ends to perform mutual movements (Figure 5B). On the initial scheme of the Sq-system the distance between the ends was 36.3 Å (Supplementary Figure S4), but in the COSM model the distance fluctuated around the mean value of 73.7 ± 16.4 Å. To test the ability of our model to predict such fine structural details of the DNA origami correctly, we performed FRET measurements on the Sq-system. The distance between the ends of the short arm obtained with FRET was very close to the modeled distance, but with a smaller variance: 73.6 ± 0.7 Å (independent samples t test: N = 100, t = 0.0609, df = 198, P = 0.951). We used the system without the scaffold strand as a negative control and in this case FRET signal corresponded to the distance of 55.5 ± 0.3 Å and was significantly different from COSM data (N = 100, t = 11.0957, df = 198, P < 0.0001). Thus, our model was entirely consistent with the FRET data.

For additional characterization of dynamic properties of the Sq-system we performed an all-atom MD simulation of the object. As in the case of the Hc-system, the all-atom model did not move far from initial coordinates within 100 ns of all-atom trajectory (Supplementary Figure S11B). Despite this, dynamic profiles of the coarse-grained and all-atom representations of the Sq-system closely resembled each other (Figure 5F), revealing that our computational model correctly simulates dynamic behavior of the system.

The all-atom model of the Sq-system consisted of 26 256 atoms (1 030 921 with solvent), which is almost 365 times greater than its COSM representation (72 particles). Computation of a 100 ns-long all-atom trajectory took 6152 min of processing time on 224 Intel Xeon E5-2697 CPU cores and 16 NVIDIA Tesla K40s GPUs, the file size of the final trajectory was 262GB in compressed form. Processing time for MD simulation of the COSM model, including two geometry optimization steps, was 14 s on a single core of Intel Core i7-6700K processor, the file size of the resulting trajectory was 10MB. Processing time for a 900-ns long trajectory (final file size: 100MB) was 132 s.

Thus, our model reliably described the global shapes, dynamics and fine details of the Sq-system. However, we had to increase the default trajectory length to obtain correct estimates of the geometric parameters.

Verification of the model based on literature data

In order to verify that our model correctly predicts structures of diverse real multilayer DNA origami we used published cryo-EM and AFM data from several recent papers (10,37–39). The only atomic structure available at the RSCB Protein Data Bank is the cryo-EM structure of a ‘pointer’ object (37). This system contains a number of unusual DNA topologies, such as a vertical stack of five Holliday junctions and custom crevices between the helices, which influence its final shape. The most interesting feature of the structure is a slight right-hand twist along the long axis of the object (Figure 6A). The feature was reproduced in a recent all-atom MD simulation (19) starting from idealized model of the structure. We have taken the approach of Maffeo and coauthors to test the ability of our model to predict the geometry of this object using caDNAno scheme as a source of initial coordinates (Figure 6B). The resulting structure was in a good agreement with the cryo-EM data (Figure 6C–E); it acquired the twisted conformation of the archetype and demonstrated the characteristic chickenwire pattern with a minimum and maximum distances between neighboring strands of 21.5 and 34.8 Å, correspondingly. The only two misrepresentations were the ‘thumb’ and ‘index finger’ geometries. Both domains have acquired a tilt deriving the ‘fingers’ from the pattern set by the square lattice (Figure 6D). While in the case of the thumb its splay was planned in design of the object using a special arrangement of crossovers, the direction of the splay in our modeling was opposite to that observed in cryo-EM. Despite this, our model reproduced the longitudinal twisting of the index finger domain. Note that in the all-atom MD study of the object the index finger domain also acquired a slight tilt in the same direction that was predicted by our model (19). For the remaining structure the coarse-grained representation traced positions and characteristic bending of individual strands with high fidelity. Calculation of a 90 ns-long trajectory for this object took 269 seconds of wall time on Intel Core i7-6700K CPU.

Figure 6.

Verification of the model using simulations of previously published DNA origami structures. (A) Different views of the cryo-EM ‘pointer’ object, PDB entry 4v5x (37); (B) idealized COSM model of the pointer object used as a starting structure for the simulation; (C–E) different views of the pointer shape after simulation (black framework), superposed with the cryo-EM coordinates (blue and pink); F: predicted structure of an ‘S’-shaped object (10). Left: side view; right: snapshots of a front view at different timeframes of the trajectory, showing out-of-bending-plane flexibility of the object; (G) a hexagonal nanocontainer opening (38); (H) an icosahedron assembly from initial system to a final shape (39). Crossover sites are not shown for clarity.

To validate the ability of our tool to predict shapes of curved origami designs, we modeled an ‘S’-shaped multilayer object designed by Kim et al. as a benchmark for the CanDo model (10). The object includes two domains of opposite curvature that determine its shape. Our model successfully reproduced the global shape of the object (Figure 6F) tracing the curvature geometry observed with AFM (10). We also captured the out-of-bending-plane twist deformation of the object shown by Kim et al. Since our tool models dynamic behavior of origami objects, we discovered that the twist is not static and fluctuates distorting the shape of the object in both directions out of the bending plane (Figure 6F). Curved designs appeared to be the most computationally intensive objects due to a number of particles required to represent their structure. Calculation of a 90 ns-long trajectory for this object took 60 minutes of wall time on the same CPU.

We were also interested in simulation of a system which is able to perform predetermined motions. An article of Douglas et al. (38) describes creation and functional testing of a nanocontainer featuring aptamer-controlled unlocking mechanism responding to the origami microenvironment. Using our tool, we modeled the opening of the nanocontainer (Figure 6G) by removing the fastening interactions between the two domains of the nanocontainer. According to our model, it is a repetitive decaying process with a period of ∼105 integration steps. While we did not calibrate our model for correct time estimates, we believe that the real motion occurs over a time interval of hundreds of ns. It is a very computationally expensive task to model the motion using existing tools, considering its estimated timescale and the size of the structure. With our model, calculation of 150 ns-long trajectory took 28 min of wall time on Intel Core i5-2500 CPU (see Supplementary Figure S9). The range of the container’s opening angles predicted by our model (135 ± 30°) was in a good agreement with the AFM data presented by Douglas et al. (38) (144 ± 24°, independent samples t test: NAFM = 10, NMD = 1250, t = 0.9309, df = 1258, P = 0.352). Note that the chickenwire pattern is clearly visible on this system (Figure 6G). From the simulated data, the co-directional crossovers were separated by 70.5 ± 1.6 Å, which is very close to an estimate of 71.4 Å, the length of a 21 bp-long B-DNA.

Another example of DNA origami functionality comes from a different work of Douglas et al., which unveiled a self-assembly of an icosahedron composed of three monomers (39). The monomers were interconnected through 15 unique complementary contacting sites. Our model provides a variety of tools for simulation of multicomponent systems. By placing weak restraints on distances between the subunits’ contacting sites, we were able to simulate hybridization between them and to observe an assembly and dynamics of the icosahedron (Figure 6H) during 30 ns of MD trajectory. The calculation took 45 minutes in terms of wall time on Intel Core i5-2500 CPU. The predicted shapes of a monomer subunit and the assembled icosahedron (Figure 6H, Supplementary Figure S3) were in a good agreement with AFM images of the object presented in the paper (39). An alternative approach could be an incorporation of additional particle types to the model, which feature increased non-covalent affinity to each other. In that case, it could be possible to simulate the system consisting of stoichiometric mixture of the subunits placed at random positions. Recently Snodin et al. simulated the self-assembly of a small DNA origami using the similar system setup (40). We plan to realize this approach in the near future.

Thus, our model was able to reliably reproduce both the equilibrium shapes and the dynamic behavior of multilayer DNA origami designs, including the curved ‘S’-shaped object. Moreover, it provides a toolset for a computational analysis of assembly of multicomponent systems.

CONCLUSION

In summary, we have developed a new coarse-grained model for DNA origami objects simulation. Our model is realized as a force field for the MD simulation program package GROMACS and uses its schemes for acceleration and parallelization of calculations. The model greatly reduces the size of simulated systems, which allows to obtain large time scale MD trajectories using scarce computational resources. Despite its simplicity, the COSM model provides reliable data on geometry and flexibility of single-layer and multilayer DNA origami objects designed in both hexagonal and square lattices and having various size and topologies. The important feature of the model is its capability to simulate the dynamic behavior of the DNA origami systems. The accuracy of the model’s prediction on geometry and fine details of the objects was in a good agreement with AFM, FRET and cryo-EM experimental data. The tool based on the model provides a valuable insight on the mechanical properties and predominant motions of the origami designs. Due to a wide selection of trajectory analysis tools provided by GROMACS, the COSM model significantly extends the potential of computational analysis of DNA origami.

AVAILABILITY

The tool is available at http://vsb.fbb.msu.ru/cosm as a web-service and as a standalone version.

Click here for additional data file.