Improved Force Fields for Peptide Nucleic Acids with Optimized

Article Cite This: J. Chem. Theory Comput. 2018, 14, 3603−3620

pubs.acs.org/JCTC

Improved Force Fields for Peptide Nucleic Acids with Optimized Backbone Torsion Parameters Maciej Jasiński,†,‡ Michael Feig,*,† and Joanna Trylska*,‡ †

Department of Biochemistry and Molecular Biology, Michigan State University, East Lansing, Michigan 48824, United States Centre of New Technologies, University of Warsaw, Warsaw, Poland

‡

J. Chem. Theory Comput. 2018.14:3603-3620. Downloaded from pubs.acs.org by UNIV OF SUNDERLAND on 11/01/18. For personal use only.

S Supporting Information *

ABSTRACT: Peptide nucleic acids are promising nucleic acid analogs for antisense therapies as they can form stable duplex and triplex structures with DNA and RNA. Computational studies of PNA-containing duplexes and triplexes are an important component for guiding their design, yet existing force fields have not been well validated and parametrized with modern computational capabilities. We present updated CHARMM and Amber force fields for PNA that greatly improve the stability of simulated PNA-containing duplexes and triplexes in comparison with experimental structures and allow such systems to be studied on microsecond time scales. The force field modifications focus on reparametrized PNA backbone torsion angles to match high-level quantum mechanics reference energies for a model compound. The microsecond simulations of PNA-PNA, PNA-DNA, PNA-RNA, and PNA-DNA-PNA complexes also allowed a comprehensive analysis of hydration and ion interactions with such systems.

■

INTRODUCTION

complexes can be more stable than DNA-DNA and DNARNA complexes so that DNA strand invasion in the presence of PNA with a complementary base sequence is likely.1,7 This has led to an interest in developing PNAs for a number of biotechnology applications8 including antisense and antigene actions, e.g., as inhibitors of bacterial translation by targeting mRNA9,10 or rRNA,11,12 as inhibitors of viral replication13 or in gene therapies targeting host DNA14−16 or microRNAs.17,18 The stronger binding affinity of PNA to DNA and RNA is believed to result primarily from reduced electrostatic repulsion between strands since the PNA backbone is neutral in contrast to the charged phosphates in the backbones of DNA and RNA. The artificial nature of the PNA backbone prevents its degradation by cellular enzymes,19 but efficient delivery of PNA remains one of the major challenges20−23 that has so far limited therapeutic applications. Therefore, the design of modified PNA8,24 or PNA-conjugates25−27 that increase the bioavailability of PNA without reducing its binding affinity to DNA and RNA continues to be important.28 High-resolution structural information about PNAPNA,29−36 PNA-DNA,5,37−39 or PNA-RNA6,40 complexes is available from crystallography 5 , 6 , 2 9 , 3 2 − 3 6 , 3 8 , 3 9 and NMR.30,31,37,40 Crystal structures of PNA-PNA and PNADNA duplexes have revealed parallel and antiparallel arrangements,5 and PNA-PNA duplexes may be right- or left-handed helices.29,32,33,35,36 Right-handed helices follow mostly canonical A-like and B-like base geometries in PNA-RNA and PNA-

Peptide nucleic acids (PNA) are synthetic nucleic acid analogs where the usual phosphate backbone in DNA or RNA is replaced by N-(2-aminoethyl)-glycine units (see Figure 1).1,2 PNAs can form duplexes with DNA, RNA, or other PNAs as well as triplexes with DNA, binding via complementary Watson−Crick or Hoogsteen base pair formation and assuming structures similar to regular nucleic acids.3−6 PNA-DNA

Received: March 23, 2018 Published: May 23, 2018

Figure 1. Peptide nucleic acid (left) vs DNA and RNA (right). © 2018 American Chemical Society

3603

DOI: 10.1021/acs.jctc.8b00291 J. Chem. Theory Comput. 2018, 14, 3603−3620

Article

Journal of Chemical Theory and Computation Table 1. Studied PNA Systems type

sequence

handedness

PDB ID resolution

box size

Na+ ions

PNA-PNA

duplex

right

0

duplex

72.4 Å

0

PNA-DNA

duplex

66.8 Å

7

PNA-RNA

duplex

69.9 Å

7

PNA-RNA

duplex

68.7 Å

7

PNA-DNA-PNA

triplex

3MBS 1.27 Å 3MBS 1.27 Å 1PDT NMR 5EME 1.14 Å 5EMF 1.14 Å 1PNN 2.5 Å

72.4 Å

PNA-PNA

p(GGCATGCC) p(GGCATGCC) p(GGCATGCC) p(GGCATGCC) p(GCTATGTC) d(GACATAGC) p(GCTGCTGC) r(GCAGCAGC) p(GCAGCAGC) r(GCUGCUGC) p(CTCTTCTTC)-HGSSGH-(CTTCTTCTC) d(GAAGAAGAG)

79.6 Å

8

system

left right right right right

calculations. Less effort has been made to reproduce the relative energetics of different conformations. Only the cyclopentanemodified PNA force field incorporates detailed conformational energy comparisons including one-dimensional torsional scans.46 Since the flexible backbone of PNA is essential for its overall structure and its adaptability to DNA and RNA when forming hybrid duplexes and triplexes a more careful parametrization remains desirable. The PNA force fields were primarily validated based on the structural stability of PNA duplexes and triplexes in comparison to experimental structures from NMR and X-ray crystallography. However, while a number of PNA-PNA structures have been available for comparison, unmodified PNA-DNA and PNA-RNA complexes were until recently only available from early NMR studies with limited resolution.37,40 Another means for validating force field parameters is via comparison to thermal melting data of PNA duplexes,48,63 but the simulation of melting and reannealing of nucleic acids requires much more conformational sampling than what has been reached in most PNA simulations to date. All of the previous simulations of PNA reported to date are in the submicroseconds regime, with most in the 1−10 ns time regime and only a few recent studies reaching 100 ns time scales.25,43,48,65 Two studies have employed enhanced sampling techniques via accelerated molecular dynamics (AMD)63 and metadynamics48 to enable melting and reannealing of PNA duplexes but with only partial success.48 Recently, new structures of PNA-RNA complexes were reported from crystallography,6 and new high-resolution structures of PNA-PNA duplexes have been available for a few years.35,38 The new structures and the ability to extend simulations to microsecond time scales creates an opportunity to reassess and further improve PNA force fields, especially with respect to the sampling of backbone torsion angles which has received less attention in previous studies. We are reporting here microsecond-scale simulations of PNA-PNA, PNA-DNA, and PNA-RNA duplexes as well as a PNA-DNA-PNA triple helix with CHARMM (Nilsson et al.49) and Amber-type (Rathinavelan and Yathindra71) force fields to validate their performance with respect to the recent experimental structures. We found that both force fields with their original parameters lead to structural distortions relative to the experimental structures when simulations are extended beyond 100 ns time scales. We thus generated updated variants of each force field, where the dihedral torsion potentials were optimized to reproduce two-dimensional torsional scans from high-level quantum mechanical calculations. In testing the updated force

DNA complexes, whereas PNA-PNA duplexes are seen to prefer P-type forms with a small twist angle, large x displacement, and major groove that is wide and deep.4,29,36 There are also PNA-PNA32 and PNA-DNA5 triplex structures, and there is one example of a unique β-turn fold for a PNA with a modified backbone.34 Computational studies of PNAs have added dynamic and energetic insights. A large number of molecular dynamics (MD) simulations have been carried out for single-stranded PNA,41−48 double-stranded PNA-PNA/PNA-DNA/PNA-RNA duplexes,31,45,48−65 and PNA-containing triplexes.66 Many of these studies consider backbone- or base-modified PNA variants44,46,47,52,55,56,58−63,65 and a conjugate of PNA with cobalamin.25 The general insight gained from these studies has been that PNA is more flexible and binds more strongly in hybrid complexes than DNA or RNA45,48,49,51,63 and that certain backbone modifications can reduce flexibility and further increase binding affinities.44,58,61,62,65 Other insights from simulations are an apparent preference of PNA for membrane surfaces43 and an increased propensity of PNA to serve as a conduit for charge transfer due to its greater flexibility.57 In addition to the MD studies, there is an early molecular mechanics analysis of PNA triple helices67 and more recent quantum chemical analyses68,69 of specific conformational features of PNA. A variety of force fields has been used to describe PNA and modified variants in published simulation studies. A CHARMM-based force field was developed by Nilsson et al.42,49 and adopted further by Weronski et al.43 Manukyan et al. rederived a CHARMM-based force field for cyclopentanemodified PNA.46,47 A PNA force field compatible with Amber force fields was initially developed by Orozco et al.51,66 and has been used in some subsequent studies.25,44,52,54 Other works use Amber-based parameters with rederived charges that are different from each other.48,53,56−59,61,63−65 Some of the parameters are for PNA with modified backbones, but there are also multiple sets of rederived parameters for unmodified PNA. One study reports optimized parameters based on the Gromos 45A4 force field45 and a very early study used PNA parameters compatible with the CVFF force field.50 Generally, all of the PNA force fields borrow bonded parameters including torsion potentials from existing force fields and focus on rederiving partial charges for the PNA backbone and in some cases also the bases.53 Most of the force field parametrization has relied on ESP (electrostatic potential)49 or RESP (restricted electrostatic potential) fitting70 of classical point charges to match electrostatic fields from quantum mechanical 3604


Article

Journal of Chemical Theory and Computation fields, the conformational sampling improved when compared to experiment and stable structures were maintained on microsecond time scales. We further took advantage of the extensive sampling of PNA duplex and triplex structures to provide a detailed comparison of water and ion distributions around PNA duplexes and triplexes. In the following, the simulation methods and force field optimization are described in detail before results are presented and discussed.

The restraints were then removed over eight rounds, each lasting 1 ns and with a force constant that was reduced in half compared to the previous round. Subsequently, 5 ns in the NPT ensemble was performed to adjust the box size with a target pressure set to 1 bar. Following equilibration, all systems were then simulated in the NVT ensemble for 1 μs. CHARMM and Amber-type force fields were used. The CHARMM force field initially consisted of the PNA parameters introduced by Nilsson et al.72 and CHARMM36 nucleic acid parameters for DNA73 and RNA.74 The CHARMM version of the TIP3P75 model was used to describe water, and Na+ parameters used recent NBFIX modifications.76 Amber force field simulations initially used the PNA parameters with partial charges computed by Rathinavelan and Yathindra71 and backbone atom types introduced by Orozco et al.66 The Rathinavelan and Yathindra set is the first complete set of published charges for PNA within the Amber family and has been used previously for single-stranded PNA44 in addition to double-stranded PNA.71 The Amber ff99OL3 parameters were used for RNA (including parmbsc0 α/γ77 and χOL378 corrections to ff9979) and Amber OL15 (with parmbsc0 α/ γ77 and χ/ζOL1,80 χOL4,81 and βOL182 modifications to ff9979) for DNA. In the Amber simulations, the standard TIP3P model75 and recommended Na+ parameters83,84 were used. Additional simulations were carried out with newly optimized force fields (see below). CHARMM and Amber simulations both applied a cutoff of 10 Å to the Lennard−Jones and short-range electrostatic interactions with a switching function becoming effective at 9 Å. Long-range electrostatic interactions were calculated using the particle mesh Ewald (PME) summation method85 under periodic boundary conditions. SHAKE was applied to allow for a 2 fs time step, and a temperature of 300 K was maintained via a Langevin thermostat with a friction constant set to 5 ps−1. The simulations with the CHARMM force field were carried out with CHARMM 86 (version 42a1) together with OpenMM87 (version 7.0.1). Amber simulations were performed with the same version of the OpenMM package but with the electrostatic energy scaling factor k = 1/(4πϵ0) adjusted slightly to match results obtained with the sander program from AmberTools16.88 One simulation each was run with the original force fields. With the optimized force fields, five replicates were carried out for 1PDT and three replicates for all other systems. Force Field Optimization. We optimized the CHARMM and Amber force fields by refitting backbone dihedral torsion parameters to match energies from quantum mechanics (QM) torsion scans. The charges, Lennard−Jones parameters, and bond and angle parameters were not modified in order to maintain compatibility with the Amber and CHARMM force field families. The model compound shown in Figure 3 was used to calculate QM torsion maps and adjust the molecular mechanics (MM) torsion parameters subsequently. This molecule covers the entire PNA backbone between two bases but does not include the bases that were replaced with methyl groups. Twodimensional (2D) maps were generated in vacuum by fixing all but two torsion angles to typical values while scanning the remaining two torsion angles in 15° increments. Two sets of conformations were considered according to the torsion values given in Table 2. For each set of dihedrals, the model compound was optimized with fixed torsions at the MP2 631+G(d) level. Energies were then evaluated at the MP2 aug-

■

METHODS Systems. We focus here on PNA-PNA, PNA-DNA, and PNA-RNA duplexes as well as a PNA-DNA-PNA triplex (see Table 1). For all of the systems, experimental structures are available from NMR (1PDT37) and X-ray crystallography (5EME/5EMF,6 3MBS,35 1PNN5) (see Figure 2). The

Figure 2. Experimental structures for the PNA-PNA (3MBS), PNADNA (1PDT), PNA-RNA (5EME, 5EMF), and PNA-DNA-PNA (1PNN) systems studied here (see Table 1). PNA strands are shown in yellow and green, DNA and RNA in gray, and amino acids in magenta (there is a peptide linker in the 1PNN triplex used in crystallization).

structures of PNA-PNA (3MBS) and PNA-RNA (5EME/ 5EMF) are very recent and resolved at resolutions below 1.3 Å, significantly better than structures used previously in comparisons between simulations and experiment. MD Simulations. All of the systems listed in Table 1 were studied via all-atom molecular dynamics (MD) simulations. The crystallographic structures were solvated in cubic boxes large enough for each atom of the nucleic acid structures to be at least 15 Å from the edge of the box. The resulting box sizes are given in Table 1. PNA strands were capped by adding an acetyl (COCH3) group to the N-terminus and an amide (NH2) group to the C-terminus. The systems were subsequently neutralized via the addition of Na+ counterions that were placed randomly by replacing water molecules. The initial systems were subjected to energy minimization; 100 steps with the steepest descent algorithm were followed by 500 steps with the adopted-basis Newton−Raphson algorithm. The systems were then gradually heated from 30 to 300 K in 10 K steps, each step lasting 20 ps. During heating, positional restraints were applied to solute heavy atoms with a force constant of 50 kcal/mol/Å2. 3605


Article

Journal of Chemical Theory and Computation

the maps below 8 kcal/mol relative to the overall minimum were considered in the RMSD calculation. The resulting RMSD was then used as the target function in a Monte Carlo optimization procedure where the force constants kn in eq 1 were varied to improve the agreement between the MM and the QM maps under the adiabatic assumption. Because torsion angles overlap in different 2D maps, all of the parameters for all of the torsion angles were fit simultaneously. The optimization run involved 10 000 Monte Carlo cycles, which was sufficient to achieve full convergence. The resulting parameters were subjected to the same umbrella sampling protocol as described above, and the entire optimization cycle was repeated once. For the CHARMM force field, an initial analysis revealed significant differences in the sampling of the α torsion (C-NC6′-C5′) between crystallographic values and the distribution obtained with the original Amber force field (data not shown). Therefore, the first attempt to improve the force field involved simply replacing the original CHARMM torsion term for α with the original Amber torsion values. This resulted in some improvement but was overall not satisfactory (data not shown). However, the resulting torsion maps from this initial modification were used as the starting point for the optimization against the QM data as described above since they were closer to the QM maps. Analysis. The analysis of the MD trajectories was performed with the MMTSB Tool Set,91 the MDAnalysis library,92,93 and the VMD visualization package.94 MINT95 and 3DNA96 were used to analyze nucleic acid structures. Hydrogen bonding was identified based on a maximum distance of 3.5 Å between the donor and the acceptor atoms and an acceptor−proton−donor angle between 150° and 180°. Stacking interactions were quantified by calculating Coulomb and van der Waals (VdW) interactions between all atoms in interacting nucleobases. Results from the CHARMM simulations were evaluated using CHARMM force field parameters, and Amber simulation results were evaluated using Amber parameters. Two nucleobases were assumed to be stacked if their vdW energy was lower than −0.5 kcal/mol. Snapshots from the 1PDT trajectories were clustered to separate the dominant helical state from partially distorted structures. There were 10 000 frames randomly selected from a total of 2 500 000 snapshots saved during the simulation with each force field. The selected frames were subsequently clustered using the dbscan algorithm from the cpptraj program in AmberTools16.88 Clusters were required to have at least 4 members, and the distance cutoff between cluster elements was set to 1.5 Å (based on heavy atoms and excluding terminal base pairs). Custom scripts and programs were used to calculate water and ion densities. Solvent densities were obtained by counting occupancies for water oxygens (with a radius of 1.5 Å) and sodium ions (with a radius of 1.4 Å) on a 0.5 Å volumetric grid after reorienting each simulation snapshot so that the central part of the nucleic acid optimally superimposes to a reference structure.

Figure 3. Model compound used for optimizing backbone torsion angles. Torsion angles optimized here are indicated with Greek letters. Their definitions are given in Table 2. Atom names are given according to the naming in the CHARMM force field. Alternate atom names in the Amber force field are also given in Table 2. Atoms are colored according to atom type (oxygen, red; nitrogen, blue; carbon, cyan; hydrogen, white). Dotted magenta line is shown to indicate the bond separating two adjacent PNA residues.

Table 2. Standard Values for PNA Backbone Torsion Angles Fixed during Dihedral Scans for Two Sets of Conformationsa

a

torsion

set 1

set 2

CHARMM

Amber

α β γ δ ε

−120 60 60 75 −20

70 60 60 75 150

C-N-C6′-C5′ N-C6′-C5′-N2′ C6′-C5′-N2′-C2′ C5′-N2′-C2′-C N2′-C2′-C-N

C*-N1*-C2*-C3* N1*-C2*-C3*-N4* C2*-C3*-N4*-C5* C3*-N4*-C5*-C* N4*-C5*-C*-N1*

For atom names see Figure 3.

cc-pVQZ level. QM calculations were carried out with Psi4 (version 1.1).89 2D torsion maps were also calculated from the MM force fields for the model compound shown in Figure 3 via umbrella sampling in vacuum. For each map, two torsion angles were varied simultaneously in 5° increments while other torsions were restrained according to the values given in Table 2 as in the QM torsion scans. Torsion angles were restrained harmonically with a force constant of 10 kcal/mol/degree2. Each umbrella simulation was run in vacuum over 2 ns using a step size of 1 fs. The WHAM program90 was used to reconstruct 2D potentials of mean force from the individual umbrellas. The differences between the initial MM torsion maps and the QM-derived maps were then minimized by sampling parameters for modified torsional terms. For each torsion, Fourier series terms according to the standard torsion potential given in eq 1 were considered

■

5

V (φ ) =

∑ kn(1 + cos(nφ − φ0,n)) n=1

RESULTS AND DISCUSSION Simulations with Established PNA Force Fields. Previous simulations of PNA systems reported stable simulations on 1−100 ns time scales.31,45,48−66 We also found stable structures up to about 200 ns, but significant structural distortions emerged on the microsecond scale based on the RMSD analysis shown in Figure 4 (see also final structures in

(1)

with multiplicities n, offsets φ0,n, and force constants kn. To compare the MM and QM maps, the QM maps were interpolated at 5° increments, and root-mean-square deviations (RMSD) were calculated after shifting maps with a constant offset that resulted in minimal RMSD values. Only regions of 3606


Article


Figure 4. RMSD time series with the original CHARMM (A) and Amber (B) force fields with respect to experimental structures after superposition, calculated for heavy atoms of PNA, DNA, and RNA nonterminal residues. For 1PDT, the first RMSD trace is with respect to the NMR structure from the PDB; the second trace (indicated by an asterisk (*)) is with respect to a refined structure obtained via simulation (see text).

Figure S1). With Amber, the PNA-PNA duplex (3MBS) lost its helical structure, including complete strand separation. With CHARMM, the PNA-PNA duplexes were maintained stably but the PNA-DNA (1PDT) and the PNA-DNA-PNA triplex (1PNN) were distorted significantly and more minor structural perturbations were found for one of the PNA-RNA complexes (5EMF). The conformational sampling with the established CHARMM and Amber PNA force fields will be described in more detail below, but it is already clear, simply based on the RMSD analysis, that both of the existing force fields have difficulties in maintaining the experimentally known PNAcontaining structures on microsecond time scales. Force Field Optimization. The apparent issues with the existing PNA force fields prompted us to consider force field modifications. Rather than changing the charges or Lennard− Jones parameters, we focused on improving the sampling of backbone torsions. Altered charges may affect force field compatibility when PNAs interact noncovalently with other molecules such as DNA or RNA, whereas torsional parameters have been the critical parameters for fine tuning the backbone conformational sampling in nucleic acids77,78,80−82,97−100 and peptides and proteins.101−104 As described in the Methods, we calculated 2D torsion maps for a model compound representing the PNA backbone in vacuum using high-level QM calculations. The torsion pairs α/ β, β/γ, and ε/δ along the backbone (see Figure 3) were considered in two conformations (see Table 2). Only one set of torsions applies for the α/ε pair with not strictly adjacent torsions (see Figure 3). The pair γ/δ was not scanned and optimized because the planar amide around N2′ results in values around either 75°/60° for right-handed PNA or −75°/− 90° for left-handed PNA.35 The MM torsion terms of the classical CHARMM and Amber PNA force fields were subsequently optimized to match the QM data. The torsion parameters before and after optimization are given in Table 3. There are significant differences for most torsion angles between the original and the optimized parameters, and there are also significant differences between the optimized CHARMM and Amber variants. These differences mainly reflect variations in the underlying force fields as the torsion terms are used here primarily as an overall correction to match the MM and QM surfaces following a similar strategy commonly applied to peptide backbones.103−105

Table 3. Original and Optimized PNA Backbone Torsion Parameters for CHARMM and Amber Force Fields CHARMM torsion

n

φ0

koriga

α

1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5

0 0 0 0 0 180 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

1.8

β

γ

δ

ε

a

0.5 0.5

kopta 0.96 0.17 0.41 0.27 1.59 0.72 0.03

Amber koriga 2.0 2.0 0.4

kopta 0.29 0.12 0.10 1.54

0.1556

0.46 0.18 0.21

0.73 0.45 1.52 0.48 0.31 0.3 0.5 0.62 1.6

0.68 0.71 0.33

0.61 0.71 0.16 0.24 0.70 0.42 0.05

0.24 0.12

Force constants in kcal/mol.

Figures 5 and 6 show the QM torsion maps in comparison with the original and optimized MM maps for pairwise torsions along the PNA backbone. It can be seen that while the overall features are similar, there are significant differences between the original MM and QM maps. There are also differences between the original CHARMM and Amber maps. One major difference is the sampling of the α torsion. The original CHARMM force field has a pronounced minimum around 180°, while Amber has the major minimum at 75° when β is near 50°. The QM surface, however, has two minima, at 75° and 200°. Other significant differences are too favorable β torsion angles near 3607


Article


Figure 5. Comparison of torsional maps between QM and MM before and after optimization for the CHARMM force field for the α/β (A), β/γ (B), ε/δ (C), and α/ε (D) torsion angles. In each panel, the top, middle, and bottom rows show the QM, initial MM, and optimized MM maps. For panels A-C, the left column refers to set 1, the right column to set 2 (see Figure 3 and Table 2 for details). Coloring indicates relative free energies in kcal/mol according to the given color bar. Energies higher than 10 kcal/mol are not shown.

180° with both CHARMM and Amber force fields compared to the QM surfaces. After optimizing the torsion potential, the CHARMM and Amber maps became more similar to the QM reference maps. The optimized torsion potential addresses in particular the discrepancies for the α and β torsion but also improves other aspects of the 2D maps. Previously, cross-correlation maps (CMAP) have been used to reproduce target torsion maps

from QM more precisely.102,105 While CMAP could be used here as well for any one map, a two-torsion CMAP cannot capture additional dependencies on other variables. For example, the deviations seen in the ε/δ torsion maps stem from simultaneously fitting to set 1 and set 2 conformations. Therefore, an exact reproduction of these maps would require a three-torsion term. However, the underlying cause for the remaining differences may be fundamentally related to 3608


Article


Figure 6. Comparison of torsional maps between QM and MM before and after optimization for the Amber force field as in Figure 5.

limitations in the fixed-charge, classical nature of the MM force field; therefore, we did not pursue CMAP or other higher order correlation terms to further match the QM energies. Interestingly, the optimized CHARMM and Amber maps are very similar to each other despite other differences in the force fields. Therefore, a comparison of the optimized force field variants is expected to provide insights into how the backbone torsion energetics drives the overall helical sampling of PNA duplex and triplex structures. Simulations with New Force Fields. Additional simulations of the systems in Table 1 were carried out with the new

force field parameters. The stability with respect to the experimental structures was improved throughout (see Figures 7 and 8 and Table 4), and the large structural distortions seen with the original force field were mostly avoided (see Figure S2). Canonical base pairing was largely maintained (see Figure 7) apart from the terminal base pairs of most structures. With Amber, the 2:15 base pair in the 1PDT (PNA-DNA) structure and the 3:14 base pair in the 5EMF (PNA-RNA) structure also lost base pairing during part of the simulations. This observation is consistent with previous reports of slightly less favorable base pairing with the Amber force field compared to 3609


Article


Figure 7. Fraction of Watson−Crick hydrogen bonding for duplexes and Watson−Crick and Hoogsteen hydrogen bonding for the triplex (left and right, respectively) along the base sequence during the simulations (A and C) and contributions of van der Waals energies to base stacking interactions (B and D) for the new CHARMM (A and B) and Amber (C and D) force fields. Combined electrostatic and van der Waals contributions are shown in Figure S3.

Figure 8. RMSD time series with the new CHARMM (A) and Amber (B) force fields with respect to experimental structures after superposition, calculated for heavy atoms of PNA, DNA, and RNA nonterminal residues. For 1PDT, the first RMSD trace is with respect to the NMR structure from the PDB; the second trace (indicated by an asterisk (*)) is with respect to a refined structure obtained via simulation. Different colors distinguish simulation replicates started from different initial velocities (rep1, black; rep2, red; rep3, green; rep4, blue; rep5, yellow).

Table 4. Average RMSD Values with Original and Modified Force Fieldsa system PNA-PNA (r) PNA-PNA (l) PNA-DNA PNA-RNA PNA-RNA PNA-DNA-PNA

PDB 3MBS 3MBS 1PDT 1PDTb 5EME 5EMF 1PNN

RMSD old ff, CHARMM [Å] 1.34 1.36 7.91 7.64 1.15 1.09 3.75

RMSD new ff, CHARMM [Å]

(

Improved Force Fields for Peptide Nucleic Acids with Optimized

Recommend Documents