Improved Force Fields for Peptide Nucleic Acids with Optimized

Subscriber access provided by UNIVERSITY OF TOLEDO LIBRARIES

Molecular Mechanics

Improved Force Fields for Peptide Nucleic Acids with Optimized Backbone Torsion Parameters Maciej Jasi#ski, Michael Feig, and Joanna Trylska J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00291 • Publication Date (Web): 23 May 2018 Downloaded from http://pubs.acs.org on May 24, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Improved Force Fields for Peptide Nucleic Acids with Optimized Backbone Torsion Parameters Maciej Jasiński1,2, Michael Feig1,*, Joanna Trylska2* 1

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

*corresponding authors Michael Feig Department of Biochemistry & Molecular Biology Michigan State University 603 Wilson Road, Room BCH 218 East Lansing, MI 48824, USA [email protected] +1-517-432-7439 Joanna Trylska Centre of New Technologies University of Warsaw Banacha 2c 02-097 Warsaw [email protected] +48 225543683

1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 45

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 parameterized with modern computational capabilities. We present updated CHARMM and Amber force fields for PNA that greatly improve the stability of simulated PNAcontaining 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.


Page 3 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION 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 Fig. 1).1-2 PNAs can form duplexes with DNA, RNA or other PNAs as well as triplexes with DNA, binding via complementary WatsonCrick or Hoogsteen base pair formation and assuming structures similar to regular nucleic acids.3-6 PNA-DNA complexes can be more stable than DNA-DNA and DNA-RNA 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 ribosomal RNA1112

, as inhibitors of viral replication13, or in gene therapies targeting host DNA14-16 or microRNAs17-18.

Figure 1: Peptide nucleic acid (left) vs. DNA and RNA (right).

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 3 ACS Paragon Plus Environment


Page 4 of 45

degradation by cellular enzymes19 but efficient delivery of PNA remains one of the major challenges2023

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 important28. High-resolution structural information about PNA-PNA29-36, PNA-DNA5, 37-39, or PNA-RNA6, 40 complexes is available from crystallography5-6, 29, 32-36, 38-39 and NMR30-31, 37, 40. Crystal structures of PNA-PNA and PNA-DNA duplexes have revealed parallel and antiparallel arrangements5 and PNAPNA duplexes may be right- or left-handed helices29, 32-33, 35-36. Right-handed helices follow mostly canonical A-like and B-like base geometries in PNA-RNA and PNA-DNA complexes, whereas PNAPNA duplexes are seen to prefer P-type forms with a small twist angle, large x-displacement, and a 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 backbone34. 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 PNA41-48, doublestranded PNA-PNA/PNA-DNA/PNA-RNA duplexes31, 45, 48-65, and PNA-containing triplexes66. 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 cobalamin25. 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 surfaces,43 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.


Page 5 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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. re-derived a CHARMM-based force field for cyclopentane-modified PNA46-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 studies25, 44, 52, 54. Other works use Amber-based parameters with re-derived charges and 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

re-derived 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 field50. Generally, all of the PNA force fields borrow bonded parameters including torsion potentials from existing force fields and focus on re-deriving partial charges for the PNA backbone and, in some cases also the bases53. 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 calculations. Less effort has been made to reproduce the relative energetics of different conformations. Only the cyclopentane-modified 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 parameterization remains desirable. The PNA force fields were primarily validated based on 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 resolution37,

40

. Another means for validating force field parameters is via comparison to thermal 5 ACS Paragon Plus Environment


Page 6 of 45

melting data of PNA duplexes48, 63, but the simulation of melting and re-annealing 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 sub-µs regime, with most in the 1-10 ns time regime and only a few recent studies reaching 100-ns time scales25, 43, 48, 65. Two studies have employed enhanced sampling techniques via accelerated molecular dynamics (AMD)63 and metadynamics48 to enable melting and re-annealing of PNA duplexes but with only partial success48. Recently, new structures of PNA-RNA complexes were reported from crystallography6 and new high-resolution structures of PNA-PNA duplexes are available for a few years35, 38. The new structures and the ability to extend simulations to microsecond time scales creates an opportunity to re-assess 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 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.


Page 7 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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/5EMF6, 3MBS35, 1PNN5) (see Fig. 2). The 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. Table 1. Studied PNA systems System

Type

Sequence

Handedness

PNA-PNA

duplex

right

PNA-PNA

duplex

PNA-DNA

duplex

PNA-RNA

duplex

PNA-RNA

duplex

PNADNA-PNA

triplex

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)

left right

PDB ID resolution 3MBS 1.27 Å

Box size

Na+ ions

72.4 Å

0

3MBS 1.27 Å 1PDT NMR

72.4 Å

0

66.8 Å

7

right

5EME 1.14 Å

69.9 Å

7

right

5EMF 1.14 Å

68.7 Å

7

right

1PNN 2.5 Å

79.6 Å

8



Page 8 of 45

Figure 2: Experimental structures for the PNA-PNA (3MBS), PNA-DNA (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 grey, and amino acids in magenta (there is a peptide linker in the 1PNN triplex used in crystallization).

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 300K in 10K 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. 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 were performed to adjust the box size with a 8 ACS Paragon Plus Environment

Page 9 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 RNA74. The CHARMM version of the TIP3P75 model was used to describe water and Na+ parameters used recent NBFIX modifications76. 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 PNA 44 in addition to double-stranded PNA71. The Amber ff99OL3 parameters were used for RNA (including parmbsc0 α/γ77 and χOL378 corrections to ff9979) and Amber OL15 (with parmbsc0 α/γ77 and χ/ζOL180, χOL481, 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 shortrange 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 CHARMM86 (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 = 1/(4 ) adjusted slightly to match results obtained with sander program from AmberTools1688. One simulation each was run



Page 10 of 45

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 re-fitting backbone dihedral torsion parameters to match energies from quantum mechanics (QM) torsion scans. The charges, LennardJones parameters, and bond and angle parameters were not modified in order to maintain compatibility with the Amber and CHARMM force field families.

Figure 3: Model compound used for optimizing backbone torsion angles. The torsion angles optimized here are indicated with Greek letters. Their definitions are given in Table 2. The atom names are given according to the naming in 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). A dotted magenta line is shown to indicate the bond separating two adjacent PNA residues. The model compound shown in Fig. 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 15o increments. Two sets of conformations were 10 ACS Paragon Plus Environment

Page 11 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 6-31+G(d) level. Energies were then evaluated at the MP2 aug-cc-pVQZ level. QM calculations were carried out with Psi4 (version 1.1).89

Table 2. Standard values for PNA backbone torsion angles fixed during dihedral scans for two sets of conformations. For atom names see Figure 3. 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*

2D torsion maps were also calculated from the MM force fields for the model compound shown in Fig. 3 via umbrella sampling in vacuum. For each map, two torsion degrees were varied simultaneously in 5 degree 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: ( ) = ∑ 1 + cos − ,

(1)

with multiplicities n, offsets , and force constants . To compare the MM and QM maps, the QM maps were interpolated at 5 degree 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 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 11 ACS Paragon Plus Environment


Page 12 of 45

function in a Monte Carlo optimization procedure where the force constants in Eq. 1 were varied to improve the agreement between the MM and 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-N-C6’-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 Set91, the MDAnalysis library92-93, and the VMD visualization package94. 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 acceptor atoms and an acceptor-proton-donor angle between 150o and 180o. 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. 12 ACS Paragon Plus Environment

Page 13 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Snapshots from the 1PDT trajectories were clustered to separate the dominant helical state from partially distorted structures. 10,000 frames were 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 AmberTools1688. 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.

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 Fig. 4 (see also final structures in Fig. 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 PNA-containing structures on microsecond time scales. 13 ACS Paragon Plus Environment


Page 14 of 45

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 *) is with respect to a refined structure obtained via simulation (see text).

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 non-covalently 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 proteins101-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 Fig. 3) were considered in two conformations (see Table 2). Only one set of torsions applies for the α/ε pair with not strictly adjacent torsions (see Fig. 3). The pair γ/δ was not scanned and optimized because the planar amide around N2’ results either in values around 75o/60o for right-handed PNA or -75o/-90o for left-handed PNA.35


Page 15 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 Tables 3 and 4. There are significant differences for most torsion angles between the original and 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.



Page 16 of 45

Table 3. Original and optimized PNA backbone torsion parameters for CHARMM and Amber force fields Torsion

n

φ0

CHARMM kopta

korig

1.8 0.5 0.5 0.45 0.3 0.5 1.6 -

0.96 0.17 0.41 0.27 1.59 0.72 0.03 0.73 1.52 0.48 0.31 0.62 0.68 0.71 0.33 -

2.0 2.0 0.4 0.1556 -

korig α

β

γ

δ

ε

a

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

Amber

a

a

kopta 0.29 0.12 0.10 1.54 0.46 0.18 0.21 0.61 0.71 0.16 0.24 0.70 0.42 0.05 0.24 0.12

force constants in kcal/mol

Figs. 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 180o while Amber has the major minimum at 75o when β is near 50o. The QM surface, however, has two minima, at 75o and 16 ACS Paragon Plus Environment

Page 17 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


200o. Another significant differences are too favorable β torsion angles near 180o 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 precisely102,

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 limitations in the fixed-charge, classical nature of the MM force field and, 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.



Page 18 of 45

Figure 5: Comparison of torsional maps between QM and MM before and after optimization for 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 Fig. 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. 18 ACS Paragon Plus Environment

Page 19 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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



Page 20 of 45

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 Figs. 7-8 and Table 4) and the large structural distortions seen with the original force field were mostly avoided (see Fig. S2). Canonical base pairing was largely maintained (see Fig. 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 QM reference interaction energies106. Van der Waals-based stacking energies were similar with both force fields and generally indicate stable duplexes with favorable base stacking (Fig. 7). The observed variations in base stacking energies between different bases reflect weaker stacking between pyrimidine vs. purine bases. Furthermore, the results agree with empirical nearestneighbor parameters for natural nucleic acids. Guanines involved in GC/CG and GG/CC base-pairs

display very stable stacking energies as found previously in RNA-RNA107 and DNA-DNA108 complexes, whereas the thymine bases with the low stacking energy in 3MBS and 1PDT with CHARMM and the adenine bases in the 5EMF with Amber are part of TA/AT and CT/GA base pairs that are the weakest base pairs in DNA duplexes108.


Page 21 of 45 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 Fig. S3. With the new force fields, the average RMSD values for non-terminal residues were close to 1 Å for the PNA-PNA and PNA-RNA duplexes. These duplexes were very stable over the entire length of the simulations with CHARMM, similar to the performance with the original force field (see Figs. 4 and 8). With the new Amber force field, the PNA-PNA systems were much more stable than before. However, in two out of three replicates for the right-handed PNA-PNA system, the helical structures became perturbed towards the end of the simulations after 800 ns and 900 ns, respectively (see Fig. 8). In these cases, canonical base-pairing was lost at one side of the duplex and single bases subsequently stacked in an alternating fashion (see Fig. S2). This suggests that the Amber force field, although much improved from the original version, still has some difficulties with maintaining stable PNA-PNA duplexes. The PNA-DNA-PNA triple helix was significantly more stable with the new CHARMM force field (2.2 Å) as none of the three replicates showed the structural deterioration seen with the original force 21 ACS Paragon Plus Environment


Page 22 of 45

field. With Amber, two replicates maintained the triple helix well, but in one simulation one of the PNA strands at the end opposite from the peptide linker melted partially after around 200 ns (see Figs. 8 and S2). This suggests again that base-pairing in the Amber force field may be slightly too weak.

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 *) is with respect to a refined structure obtained via simulation (see text). 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 fields System

PDB

PNA-PNA (r) PNA-PNA (l) PNA-DNA

3MBS 3MBS 1PDT 1PDT* 5EME 5EMF 1PNN

PNA-RNA PNA-RNA PNA-DNA-PNA

RMSD old ff, CHARMM [Å] 1.34 (

Improved Force Fields for Peptide Nucleic Acids with Optimized

Recommend Documents