Coarse-Grained Molecular Dynamics Model of Double-Stranded DNA

Apr 27, 2017 - the nonbond, stack, angle bending, and electrostatic interaction. The average twisted angle and the persistence length of the model wit...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCB

Coarse-Grained Molecular Dynamics Model of Double-Stranded DNA for DNA Nanostructure Design Hiromasa Yagyu,*,† Jae-Young Lee,‡ Do-Nyun Kim,‡ and Osamu Tabata§ †

Department of Mechanical Engineering, Kanto Gakuin University, 1-50-1 Mutsuura-higashi, Kanazawa-ku, Yokohama 236-8501, Japan ‡ Department of Mechanical and Aerospace Engineering, Seoul National University, 1 Gwanak-ro, Daehak-dong, Gwanak-gu, Seoul 08826, Republic of Korea § Department of Micro Engineering, Kyoto University, Kyoto Daigaku-Katsura C3, Nishikyo-ku, Kyoto, 606-8501, Japan ABSTRACT: A new coarse-grained molecular dynamics double-stranded DNA model (nCG-dsDNA model) using an improved beads−spring model was proposed. In this model, nucleotide comprising phosphate, sugar, and base group were replaced by a single bead. The double stranded model with 202 base pairs was created to tune the parameters of the bond, the nonbond, stack, angle bending, and electrostatic interaction. The average twisted angle and the persistence length of the model without electrostatic interaction were calculated at 35.3° and 120.3 bp, confirming that the proposed model successfully realized the experimentally observed double-stranded DNA structure. Moreover, the model with electrostatic interaction was discussed. From calculation results, we confirmed that the dependency of the salt concentration on the persistence length of the nCG-dsDNA model at the 30% charge is in good agreement with the Poisson−Boltzmann theoretical model.

■

INTRODUCTION DNA origami is the most promising method for developing nanoscale structures (nanostructures), which can be used selfassemble functions. This technology, which is applied to combine a large number of single stranded DNA chains using complementarity of bases in the DNA chain, enables the production of two-dimensional1 and three-dimensional nanostructures.2 A DNA molecule has repeating units called nucleotides, which comprise phosphate, sugar, and a base group. The chain backbone is constructed from phosphate and the sugar group. The base groups (adenine, thymine, cytosine, and guanine) are connected to sugar group. A base pairing (adenine with thymine or cytosine with guanine) with a planar stack interaction between the bases and bonds within the backbone leads to the double helix structure. Recently, nanoscale mechanical parts were constructed with DNA chains.3 However, an inclusive knowledge of their mechanical properties, especially of their strong dependency on the nick (backbone discontinuity) and crossover (bridge between neighboring chains) in the sequence is required to use DNA origami structures effectively as larger systems.4 To address this, we develop a new simulation tool to analyze the relation between structure and mechanical properties in complex DNA origami structures with a few thousand base pairs. The molecular dynamics method is expected to be able to calculate the mechanical properties of double-stranded DNA molecules. However, a full atomistic model cannot simulate DNA chains with more than a thousand base pairs due to the © 2017 American Chemical Society

long calculation time required. Coarse-grained molecular dynamics simulations of polymer melts have been performed using the Kremer−Grest model (bead−spring model).5 The Kremer−Grest model was replaced several monomer units by a single bead, and the used potential energies were only two types (bond and nonbond potential). By this coarse graining, the model cannot account for the chemical effects of the chain. However, chain dynamics with long-term scale can be simulated.6,7 Several coarse-grained DNA models have been proposed.8−16 Ding16 reported the interwound structures using the Janus chain model which builds on a system of semiflexible excluded volume multibead chains. Doi8 and Korolev11 reported a coarse-grained model of a double-stranded DNA chain, where the nucleotides (phosphate, sugar, and base groups) were treated as a single bead. To construct the helix structure, virtual beads connected to backbone beads and a harmonic spring were used to produce a stack interaction, i.e., an attractive force between the bases (π−π interaction). However, this model was not suitable for a DNA chain model because the stack interaction in the actual chain is a nonbonding interaction, and the harmonic potential represented by connecting spring was used as stacking interactions of nonbonding interaction. Savelyev10 used a stack interaction as a nonbond potential. However, this model included several nonbonding pairs; consequently, the required calculation time Received: April 26, 2017 Published: April 27, 2017 5033

DOI: 10.1021/acs.jpcb.7b03931 J. Phys. Chem. B 2017, 121, 5033−5039

Article

The Journal of Physical Chemistry B was long because a pair of beads had to be identified for all beads in the molecular dynamics calculation. The oxDNA model reported for DNA origami12 also requires a long calculation time since this model comprises a backbone bead and a base bead, and phosphate and sugar groups are represented in a single backbone bead. The stack interaction is defined as a harmonic bond and a Morse bond between the base beads. The new coarse-grained molecular dynamics double-stranded DNA model (nCG-dsDNA model) proposed in this study represents a B-DNA chain using a simple bead−spring model. In this model, each nucleotide comprising phosphate, sugar, and a base group is replaced by a single bead; the different sequences in the molecule correspond to the beads that are classified into four types: adenine (A), thymine (T), cytosine (C), and guanine (G). In the molecular dynamics simulation, the calculation time is proportional to N2(O(N2)) for using all electrostatic and nonbond interactions, and to N(O(N)) for using the Cutoff method with neighbor lists. Where N is the number of beads. Therefore, it is expected that a large chain can be simulated within a short calculation time because the number of beads in the nCG-dsDNA model is half that in the conventional DNA model (oxDNA model). Therefore, the proposed approach is suitable for simulating DNA origami structures requiring a large number of beads. In this article, parameters, such as the bond and nonbond potentials for the nCG-dsDNA model without electrostatic interaction, were optimized through comparison with a double helix structure of a realistic DNA molecule. The persistence length of the optimized nCG-dsDNA model was compared with that of the oxDNA model. Moreover, the nCG-dsDNA model with electrostatic interaction was created using conversion units, and the persistence length of the model was compared with the Poisson−Boltzmann theoretical model.

Figure 1. Molecular structure of the double-stranded DNA chain and the nCG-dsDNA model. Chains 1 and 2 were stranded and bonded with a hydrogen bond.

beads (adenine (A), thymine (T), cytosine (C), and guanine (G)) were used. Six types of potential energy were applied to the model, as shown in Figure 2.

■

MODELING AND SIMULATION METHOD We created for the double-stranded DNA model with 202 base pairs, as shown in Table 1. Figure 1 shows the molecular structure of the double-stranded DNA chain and the nCGdsDNA model. In this study, nucleotides (phosphate, sugar, and base groups) were treated as a single bead. Four types of

Figure 2. Potential energies of the model. In the stack interaction between base pairs, each base pair was utilized as a single nonbonding site, and the interaction was defined as the Lennard−Jones (LJ) potential. The left beads at outermost (k = 0) indicate to the chain end.

Table 1. Sequence of the DNA Model chain

sequence

chain 1

5′-GCCGGCCGCCGGCGCGGCGGAGCG GCTGCCAGTTCCCAGCGCGGGCGCGG CTGCCGCGCCGGAGCTCGTCGCCCGC AGGAGACCCCCGTGCCGCGCGGTCTC GGGCCCACGGGCCCGGACCGCGGCCA CACCCTCATGGCCACGCGGGCCGGTG TTGCCTCCCCGACGCGCCGCCAAATC GCAGTCGGGGGACCCCGGCGCT-3′

The bead−spring model was used for creating coarse-grained double-stranded DNA chain. The basic bead−spring model utilizes a polymer chain that comprises beads connected with bond and nonbond potential energies.5 In this work, the improvement of the basic bead−spring model by including the potentials for angle bending, stack interaction, base pairs bond, and electric interaction was conducted to simulate the doublestranded DNA. In this model, the following potential energy U of the model were used: U = U bond + U angle + U nbond + U stack + U bp + U elec (1)

chain 2

3′-CGGCCGGCGGCCGCGCCGCCTCGC CGACGGTCAAGGGTCGCGCCCGCGCC GACGGCGCGGCCTCGAGCAGCGGGCG TCCTCTGGGGGCACGGCGCGCCAGAG CCCGGGTGCCCGGGCCTGGCGCCGGT GTGGGAGTACCGGTGCGCCCGGCCAC AACGGAGGGGCTGCGCGGCGGTTTAG CGTCAGCCCCCTGGGGCCGCGA-5′

where the individual contributions are 1 U bond(rij) = k b(rij − r0)2 2 U angle(rijk) =

1 kθ(θijk − θ0)2 2

(2) (3)

where 5034

DOI: 10.1021/acs.jpcb.7b03931 J. Phys. Chem. B 2017, 121, 5033−5039

Article

The Journal of Physical Chemistry B cos θijk =

(rj − ri) ·(rj − rk)

U elec(rij) =

rijrjk

(4) bond

εk rij

e−rij / λD (7)

where

These potentials are the bond potential U (rij) and the angle-bending potential Uangle(rijk). The bond between the two beads in the backbone was defined as the harmonic potential. The angle bending between three beads (j = i+1, k = i+2) was controlled by the harmonic bending potential. Here, rij is the distance between beads i and j. θijk is the inner angle. kb and r0 are spring constants and the equilibrium bond length. kθ and θ0 are angle spring constant and the bond angle in the equilibrium state. The effect of the bond and the angle-bending potential is discussed in the next section. The short-range interaction between two beads separated by a distance rij is defined as the nonbond potential Unbond(rij) using a repulsive Lennard-Jones (LJ) potential. Here, the stack potential Ustack(rij) of a DNA chain, corresponding to a π−π interaction, is realized by treating each base pair as a single nonbond site using the LJ potential with attractive force.

⎛ ε ε k T ⎞1/2 λD = ⎜ 0 k 2B ⎟ ⎝ 2NAe I0 ⎠

(8)

In eq 7, nominal charges qi and qj are located at a distance rij. ε0 and εk are the permittivity of the vacuum and the relative permittivity of the aqueous solution, respectively. In eq 8, λD, NA, I0, and T are the Debye length, Avogadro’s number, the ionic strength of the salt, and temperature of the system, respectively. The charge in the coarse-grained DNA model is chosen to be 40% of the nominal charge density,17 and the nominal charge qi is calculated by the following equation. qi = Q i

NAe 2 4πε0σε

(9)

where, Qi is a realistic charge of each bead, and the charges are fixed on Qi = −1.0 C. σ and ε are length unit and energy unit. In this study, the effective charges qieff on each bead was found by changing the negative charge in the range of 20−50% of the nominal charge density (100%). The model without electrostatic interaction uses only the short-ranged excluded volume by the repulsive LJ potential. For this reason, we fit to experimental data at [Na+] = 500 mM. At this ionic concentration, the Debye length is approximately 0.43 nm. In contrast, the diameter of the excluded volume of backbone−backbone interactions in our model is below 0.5 nm. The Debye length at [Na+] = 500 mM is smaller than the diameter of excluded volume. The time evolution of the bead at a position ri was calculated using the following Langevin equation:

where σij and εij are the distance and energy parameters between beads i and j, and rc is the cutoff distance, respectively. For the stack potential Ustack(rij), nonbond site of two beads of base pair was treated as single site, and rij was defined as the distance between the geometrical centers of the base pair. The LJ potential (stack) is acting between the geometrical centers of base pairs. The chains could be twisted like a circular arrow in Figure 2 by balancing the forces between the beads in the backbone (bond) and those between the stack sites (stack) as shown in Figure 2. The 1−2 and 1−3 nonbond interactions, which is defined as the connected beads separated by two and three bonds, were not calculated. Base pair bonds between beads in a different backbone (chain) are defined as hydrogen bonds, as indicated by the Morse function: U bp(rij) = D[exp[−α(rij − l0)] − 1]2 − D

qiqj

m

d 2 ri dt

2

=−

dr ∂U − Γ i + W( i t) dt ∂ri

(10)

where m is the mass of a bead of the system, Γ is the friction constant. W i (t) is Gaussian white noise vector with independent components, which is defined as the following equation:

(6)

where l0 is the bond length in the equilibrium state. D is the binding energy of the hydrogen bonds. Hydrogen bonds in a GC pair were set to 1.5 times stronger in comparison with those in an AT pair because the number of hydrogen bonds is 2 bonds for AT pair, and 3 bonds for GC pair, as shown in Figure 1. α is a parameter depending on the force constants 2α2D. In this model, the force constants of a GC and an AT pair were considered constant. The electric interactions are represented by the cutoff Coulomb potential from the Debye−Hückel theory Uelec(rij). To consider the electrostatic interaction in the model, Debye length should be defined by reduced units. The coarse-grained molecular dynamics simulation uses reduced units, for which units of length, energy, mass, time, and temperature are defined as σ, ε, m, τ (=σ(m/ε)1/2), and ε/kB, respectively. Here, kB is Boltzmann constant. In this article, first, the model without electrostatic interaction was created, and reduced units was converted to SI unit. Subsequently, the model with electrostatic interaction was discussed.

⟨Wi (t )Wj(t ′)⟩ = 2kBTm ΓδijIδ(t − t ′)

(11)

where δ is the delta function, and I is the identity matrix. In this study, the friction constant Γ and the temperature T were set to 0.5τ−1and 1.0ε/kB, respectively.5 The velocity Verlet algorithm with a time step Δt = 0.005τ was used for solving the equation of motion of the beads. All walls of the cell were defined as the periodic boundary conditions. The cell size was set to 20σ × 20σ × 250σ for the model without electrostatic interaction, and 500σ × 500σ × 500σ for the model with electrostatic interaction. The system was calculated using the NVT ensemble, which used constant temperature and density condition. All the calculations were conducted using the coarse-grained molecular dynamics program COGNAC in the OCTA package distributed as an open source code.18

■

RESULTS AND DISCUSSION nCG-dsDNA Model without Electrostatic Interaction. Simulation Parameters. In this section, the parameters of the model without electrostatic interaction are defined. In the 5035

DOI: 10.1021/acs.jpcb.7b03931 J. Phys. Chem. B 2017, 121, 5033−5039

Article

The Journal of Physical Chemistry B

realistic double-stranded DNA chain. A strand of the DNA chain usually circles the axis of helix once every 10 base pairs.19 Figure 5 shows the mean twisted angle as a function of the binding energy of the Morse potential function of an AT base

definition of the simulation parameters, it was important that the model maintained the helix structure of the chain during the calculation. In this study, the simulation was run for 5 000 000 calculation steps (25 000τ), and the mean twisted angle of the model was evaluated for selecting the simulation parameters. In this model, the equilibrium length r0 of the bond potential and the equilibrium angle θ0 of the angle-bending potential were set to 0.7σ and 0°, respectively. αAT and l0 in the Morse potential were defined as 3.25σ−1 and 1.1σ, respectively. σij, εij, and c of the nonbond interaction in the LJ potential were set to 1.0σ, 1.0ε, and 2.5σij, respectively. σij, εij, and c for the stack interaction in the LJ potential were set to 0.26σ, 20ε, and 5.0σij, respectively. In the electric static interaction, the permittivity ε0εk was set to 1.0. The mean twisted angle of the model was calculated using the following equation. θ=

1 N−1

N

⎛ v ·ivi + 1 ⎞ ⎟ ⎝ |vi||vi + 1| ⎠

∑ cos−1⎜ i=1

(12)

where vi is the vector between base pairs i and N is the number of base pairs in the chain. The mean twisted angle was calculated by averaging the angles between the nearest neighbor vectors vi and vi+1 of the base pair, as shown in Figure 3.

Figure 5. Mean twisted angle as a function of the binding energy of the Morse potential function.

pair with kθ = 20.0ε. From the bonding energy of 192ε, the mean twisted angle was almost at an equilibrium state. Therefore, we selected DAT = 192ε in the parameters for the nCG-dsDNA model. Table 2 summarizes the simulation Table 2. Parameters of the nCG-dsDNA Model potential Ubond

Figure 3. Schematic of the twisted angle of the nCG-dsDNA model. θ is the twisted angle of the nearest neighbor base pairs.

Uangle Ustack Unbond Ubp

In this study, the attractive energy εij between stacks was fixed to 20ε and an analysis of the relation between the spring constant kθ of the angle-bending potential and spring constant kb of the bond potential was conducted. Figure 4 shows the

model 1 kb = 120εσ−2 r0 = 0.7σ kθ = 10.0ε θ0 = 0

model 2

model 3

model 4

kb = 120εσ−2 kb = 150εσ−2 kb = 150εσ−2 r0 = 0.7σ r0 = 0.7σ r0 = 0.7σ kθ = 12.5ε kθ = 17.5ε kθ = 20.0ε θ0 = 0 θ0 = 0 θ0 = 0 σij = 0.26σ, εij = 20ε, rc = 1.3σ σij = 1.0σ, εij = 1.0ε, rc = 21/6σ l0 = 1.1σ, DAT = 192ε, αAT = 3.25σ−1

parameters. In a coarse-grained molecular dynamics simulation, a model of a polymer chain such as that in rubber materials uses only two potentials, i.e., the bond and nonbond potentials. The bond potential is defined as a harmonic bond potential, and the spring constant is set to kb = 200εσ−2.6,7 In this work, the simulation was run for 5 000 000 calculation steps (25 000τ), and the persistence length of the model was analyzed using the trajectory data. The persistence lengths of the three models in Table 2 were evaluated for realizing the stiffness of the DNA chain. Model Structure. Figure 6 shows a created nCG-dsDNA model with 202 base pairs. The twisted angle of the model, which is the angle between vectors of neighboring base pairs,

Figure 4. Mean twisted angle as a function of the spring constant of the angle-bending potential.

mean twisted angle as a function of the spring constant of the angle-bending potential. The mean twisted angle linearly decreased with the angle-bending potential kθ. When the angle-bending potential kθ was 10.0ε and 12.5ε for kb = 120εσ−2, and 17.5ε and 20.0ε for kb = 150εσ−2, we confirmed that the mean twisted angle approached 35°−36° of the

Figure 6. Snap shot of the nCG-dsDNA model 4. Chain 1 (red) and Chain 2 (green) were twisted by stack interaction. 5036

DOI: 10.1021/acs.jpcb.7b03931 J. Phys. Chem. B 2017, 121, 5033−5039

Article

The Journal of Physical Chemistry B

GHz CPU. Although the calculation time in molecular dynamics simulation depends on a programing code, our model was attributed to minimizing a calculation time to a quarter of the oxDNA model, which is shorter time comparing with that according to the rule of O(N). This means that our model has the validity of the nCG-dsDNA model as discussed in the Introduction Section. Persistence Length. The persistence length Lps of the nCGdsDNA model, indicating the stiffness of the chain, is compared with that of the oxDNA model using the same sequence. Here, the oxDNA of the version incorporates electrostatics only through the short-ranged excluded volume, and can be simulated double stranded DNA at high salt concentration of around [Na+] = 500 mM. The persistence length Lps is a basic mechanical property reflecting the stiffness of a polymer. When Lps is small, the chain shows low stiffness. Lps is defined by calculating the average dot production of the unit vectors between the midpoints of a base pair, and by fitting the results to the exponential equation ⟨nk·n0⟩ = exp (−0.34k/ Lps).20 k is defined as the maximum base pair number for using calculations. nk is the unit vector between a base pair at k and k + 1. The analysis starts at the 10th base pair n10 and calculates the time-averaged dot production, indicated to the average of dot production during the calculation of correlation each pair from n10 with n11 up to n50. Figure 7 shows the time-averaged dot production of the unit vector between the midpoints of the stack site ⟨nk·n0⟩ in the

was 35.6° for model 1, 35.1° for model 2, 35.8° for model 3, and 35.3° for model 4. These values are in good agreement with the 36° in the realistic double-stranded DNA chain. In model 4, the average distance between the base pairs (stacks) connected by the LJ potential st, average distance between the beads in the backbone connected by the harmonic potential of model bc, and average distance between the backbone beads connected by the Morse potential hb are st = 0.30σ, bc = 0.45σ, and hb = 1.11σ, respectively. The length unit σ is defined as 0.34/st = 0.34/0.30 nm = 1.13 nm because the distance between the base chains in a realistic DNA molecule is 0.34 nm.19 The distance between the beads in the backbone l is l = bc × σ = 0.45 × 1.12 = 0.50 nm. From these calculations, the diameter of the helix was calculated as hb + σ = 2.36 nm based on the average distance between the backbone beads hb and diameter of a bead σ, indicated as the length unit. This value is close to the actual helix diameter of the double-stranded B-DNA chain. In this model, since a single bead is used to approximate phosphate, sugar, and base groups, the weight of the bead, indicated as the unit of mass m, is calculated at m = 308.95 amu (atomic mass unit) = 5.13 × 10−22 g using the sum of three groups. The molecular mass of the base group is defined as the average of adenine, thymine, cytosine, and guanine. Next, since the testing temperature (296 K) is set at 1.0ε/kB in the simulation, ε is calculated at ε = 4.09 × 10−21 J using kB = 1.38 × 10−23 J/K. From this scaling, the time τ = σ (m/ε)1/2 is calculated at τ = 1.25 × 10−11 s. Therefore, in this study, the required calculation time for 5 000 000 steps (25 000τ) is 3.14 × 10−7 s (0.314 ms). Table 3 shows the conversion results of the simulation units for the models. Table 3. Conversion Results of the Simulation Units for the Models model model model model model

1 2 3 4

length [m]

energy [J]

mass [g]

time [sec]

1.13 × 10−9

4.09 × 10−21

5.13 × 10−22

1.27 × 10−11

1.12 × 10−9

4.09 × 10−21

5.13 × 10−22

1.25 × 10−11

Table 4 shows the standard deviation of the distance between the connected beads in the backbone SR, the stack distance Sst, Figure 7. Time-averaged dot production of unit vector between the midpoints of the stack site ⟨nk·n0⟩ in the nCG-dsDNA model 4 as a function of the maximum base pairs number.

Table 4. Standard Deviation of the Distance between the Connected Beads in the Backbone SR, the Stack Distance Sst, Distance between Base Pairs Sl, and the Twisted Angle Sθ during the Calculation model model model model model

1 2 3 4

SR [nm]

Sst [nm]

Sl [nm]

Sθ [degree]

0.00221 0.00216 0.00215 0.00205

0.00107 0.00095 0.00128 0.00131

0.00109 0.00119 0.00113 0.00100

0.287 0.281 0.277 0.258

nCG-dsDNA model 4 as a function of the number of base pairs. Lps is 40.9 nm (120.3 bp) for the nCG-dsDNA model. In contrast, Lps is 41.7 nm (122.7 bp) for the oxDNA model, as shown in Figure 8. Table 5 shows the results of the persistence length Lps analysis. The experimentally calculated persistence length of the double-stranded DNA chain Lps is 50 nm.21 The result for the nCG-dsDNA model 4 are in good agreement with that of oxDNA model. nCG-dsDNA Model with Electrostatic Interaction. Simulation Parameters of Electrostatic Interaction. The electrostatic interaction was added to the model 4, and the persistence length at a salt concentration of [Na]+ = 1, 10, 100, and 500 mM was calculated. Table 6 shows Debye length λD at each salt concentration. The nominal charge density was defined as qi = −7.1 from eq 9 using the length unit and the energy unit as shown in Table 3, and the effective charge qieff

distance between base pairs Sl, and the twisted angle Sθ during the calculation. The distributions of SR, Ss, and Sl in all models showed small value. In contrast, Sθ was relatively large and its value of model 4 was smallest. Therefore, we concluded that model 4 is a suitable for the double stranded DNA model. Calculation Time. The calculation time for 5 000 000 calculation steps (25 000τ) was 5 h for the nCG-dsDNA model, and 24 h for the oxDNA model using Intel Xeon 2.5 5037

DOI: 10.1021/acs.jpcb.7b03931 J. Phys. Chem. B 2017, 121, 5033−5039

Article

The Journal of Physical Chemistry B

using a bare (nonelectrostatic) persistence length of 50 nm.22 The persistence length, resulting from our models decreased with increasing a salt concentration. The tendency of the result at the 30% charge, agrees with the Poisson−Boltzmann theoretical model. At the salt concentration of 1 mM, the persistence length of oxDNA model increased in comparison with that of the Poisson−Boltzmann theoretical model, our model. The nCG-dsDNA model has a wide range of salt concentration in comparison with the oxDNA model. Figure 10

Figure 8. Time-averaged dot production of unit vector between the midpoints of the stack site ⟨nk·n0⟩ in the oxDNA model as a function of the maximum base pairs number.

Table 5. Persistence Length Lps of the nCG-dsDNA Model and the oxDNA Model model nCG-dsDNA model nCG-dsDNA model nCG-dsDNA model nCG-dsDNA model oxDNA model

1 2 3 4

Lps [bp]

Lps [nm]

126.5 124.4 135.3 120.3 122.7

43.0 42.3 46.0 40.9 41.7

Figure 10. Twisted angle of the nCG-dsDNA model with electrostatic interaction as a function of salt concentration.

shows the twisted angle as a function of salt concentration. The twisted angle of the model at 40% and 50% charge decreased at low salt concentration. This indicated to breaking of the double helix structure. Therefore, these results confirmed that nCGdsDNA model at 30% charge explains the salt concentration dependency of behavior of double stranded DNA, and it was expected that our model have the validity of the nCG-dsDNA model for DNA origami.

Table 6. Debye Length λD of the nCG-dsDNA Model 4 salt concentration [mM]

λD [nm]

1 10 100 500

9.71 3.07 0.97 0.43

■

CONCLUSIONS The parameters of the nCG-dsDNA model without electrostatic interaction were optimized through comparison with the double helix structure of a realistic DNA molecule. The persistence length of the optimized nCG-dsDNA model was compared with that of the oxDNA model. The persistence length of the nCG-dsDNA model without electrostatic interaction was calculated as 120.3 bp (40.9 nm). The result was agreed well with the persistence lengths of the oxDNA model and the experiments. Moreover, the nCG-dsDNA model with electrostatic interaction was created using converted units. The persistence length of the model was compared with the oxDNA model and the Poisson−Boltzmann theoretical model. The result confirmed that the persistence length at the 30% charge, agrees with the Poisson−Boltzmann theoretical model. In a future study, mechanical properties, such as the stretch, twist, and bend, of the nCG-dsDNA model without nick will be evaluated. In addition, the dependencies of the nick and crossover in the sequence on the mechanical properties will be analyzed using the nCG-dsDNA model.

was changed in the range of 20−50% of the nominal charge density. Effect of Salt Concentration of Persistence Length. Figure 9 shows the persistence length of the nCG-dsDNA model as a function of salt concentration. A filled diamond symbol is the persistence length of the oxDNA model using the same sequence. Here, the oxDNA of the version incorporates electrostatics using the Debye−Hückel theory (oxDNA2). A dot line was predicted from the Poisson−Boltzmann theory

■

AUTHOR INFORMATION

Corresponding Author

*Tel: 81-45-786-7118; Fax: 81-45-786-7118; E-mail: yagyu@ kanto-gakuin.ac.jp.

Figure 9. Persistence length of the nCG-dsDNA model with electrostatic interaction. A filled diamond symbol is the persistence length of the oxDNA model. A dot line is predicted from the Poisson− Boltzmann theory.22

ORCID

Hiromasa Yagyu: 0000-0002-1256-2006 5038

DOI: 10.1021/acs.jpcb.7b03931 J. Phys. Chem. B 2017, 121, 5033−5039

Article

The Journal of Physical Chemistry B Notes

(20) Ouldridge, T. E.; Louis, A. A.; Doye, J. P. K. Structural, Mechanical and Thermodynamic Properties of a Coarse-grained DNA Model. J. Chem. Phys. 2011, 134, 085101. (21) Bustamante, C.; Bryant, Z.; Smith, S. B. Feature Ten Years of Tension: Single-molecule DNA Mechanics. Nature 2003, 421, 423− 427. (22) Markegard, C. B.; Fu, I. S.; Reddy, A.; Nguyen, H. D. Coarsegrained Simulation Study of Sequence Effects on DNA Hybridization in a Concentrated Environment. J. Phys. Chem. B 2015, 119, 1823− 1834.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS A part of this work was supported by a grant from Japan Society for the Promotion of Science (JSPS), Bilateral Joint Research Projects (Japan-Republic of Korea).

■

REFERENCES

(1) Rothemund, P. W. K. Folding DNA to Create Nanoscale Shapes and Patterns. Nature 2006, 440, 297−302. (2) Seeman, C. Feature DNA in a Material World. Nature 2003, 421, 427−431. (3) Dietz, H.; Douglas, S. M.; Shih, W. M. Folding DNA into Twisted and Curved Nanoscale Shapes. Science 2009, 325, 725−730. (4) Ma, Z.; Kim, Y. J.; Kim, D. N.; Tabata, O. Encyclopedia of Polymeric Nanomaterials; Springer-Verlag, 2015; pp 589−603. (5) Kremer, K.; Grest, G. S. Dynamics of Entangled Linear Polymer Melt: A Molecular-dynamics Simulation. J. Chem. Phys. 1990, 92, 5057−5086. (6) Yagyu, H.; Utsumi, T. Coarse-grained Molecular Dynamics Simulations of Nanofilled Crosslinked Rubber. Comput. Mater. Sci. 2009, 46, 286−292. (7) Yagyu, H. Coarse-grained Molecular Dynamics Simulation of the Effects of Strain Rate on Tensile Stress of Crosslinked Rubber. Soft Mater. 2015, 13, 263−270. (8) Doi, K.; Haga, T.; Shintaku, H.; Kawano, S. Development of Coarse-graining DNA Models for Single-nucleotide Resolution Analysis. Philos. Trans. R. Soc., A 2010, 368, 2615−2628. (9) Potoyan, D. A.; Savelyev, A.; Papoian, G. A. Recent Successes in Coarse-grained Modeling of DNA. Comput. Mol. Sci. 2013, 3, 69−83. (10) Savelyev, A.; Papoian, G. A. Chemically Accurate Coarse Graining of Double-stranded DNA. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 20340−20345. (11) Korolev, N.; Luo, D.; Lyubartsev, A. P.; Nordenskiold, L. A Coarse-grained DNA Model Parameterized from Atomistic Simulations by Inverse Montecarlo. Polymers 2014, 6, 1655−1675. (12) Doye, J. P. K.; Ouldridge, T. E.; Louis, A. A.; Romano, F.; Šulc, P.; Matek, C.; Snodin, B. K.; Rovigatti, L.; Schreck, J. S.; Harrison, R. M.; Smith, W. P. J. Coarse-graining DNA for Simulations of DNA Nanotechnology. Phys. Chem. Chem. Phys. 2013, 15, 20395−20414. (13) Hinckley, D. M.; Freeman, G. S.; Whitmer, J. K.; de Pablo, J. An Experimentally-informed Coarse-grained 3-site-per-nucleotide Model of DNA: Structure, Thermodynamics, and Dynamics of Hybridization. J. Chem. Phys. 2013, 139, 144903. (14) Naômé, A.; Laaksonen, A.; Vercauteren, D. P. A Coarse-grained Simulation Study of the Structures, Energetics, and Dynamics of Linear and Circular DNA with its Ions. J. Chem. Theory Comput. 2015, 11, 2813−2826. (15) Dans, P. D.; Darre, L.; Machado, M. R.; Zeida, A.; Brandner, A. F.; Pantano, S. Assessing the Accuracy of the SIRAH Force Field to Model DNA at Coarse Grain Level. In Advances in Bioinformatics and Computational Biology; Setubal, J. C., Almeida, N. F., Eds.; Springer International Publishing: Gewerbestrasse, Switzerland, 2013; pp 71− 81. (16) Ding, Y.; Kröger, M. Phase Behavior and Formation Dynamics of Helically Wound Networks: Generalized Janus Chain Model. Macromolecules 2009, 42, 576−579. (17) Ramachandran, A.; Guo, Q.; Iqbal, S. M.; Liu, Y. Coarse-grained Molecular Dynamics Simulation of DNA Translocation in Chemicallymodified Nanopores. J. Phys. Chem. B 2011, 115, 6138−6148. (18) Aoyagi, T.; Sawa, F.; Shoji, T.; Fukunaga, H.; Takimoto, J.; Doi, M. A General-purpose Coarse-grained Molecular Dynamics Program. Comput. Phys. Commun. 2002, 145, 267−279. (19) Watson, J. D.; Crick, F. H. C. A Structure for Deoxyribose Nucleic Acid. Nature 1953, 171, 737−738. 5039

DOI: 10.1021/acs.jpcb.7b03931 J. Phys. Chem. B 2017, 121, 5033−5039